Negative Lags on the Viscous Timescale in Quasar Photometry and Prospects for Detecting More with LSST

The variability of quasar light curves can be used to study the structure of quasar accretion disks. For example, continuum reverberation mapping uses delays between variability in short and long wavelength bands ("short"lags) to measure the radial extent and temperature profile of the disk. Recently, a potential reverse lag, where variations in shorter wavelength bands lag the longer wavelength bands at the much longer viscous timescale, was detected for Fairall 9. Inspired by this detection, we derive a timescale for these"long"negative lags from fluctuation propagation models and recent simulations. We use this timescale to forecast our ability to detect long lags using the Vera Rubin Legacy Survey of Space and Time (LSST). After exploring several methods, including the interpolated cross-correlation function, a Von-Neumann estimator, javelin, and a maximum-likelihood Fourier method, we find that our two main methods, javelin and the maximum-likelihood method, can together detect long lags of up to several hundred days in mock LSST light curves. Our methods work best on proposed LSST cadences with long season lengths, but can also work for the current baseline LSST cadence, especially if we add observations from other optical telescopes during seasonal gaps. We find that LSST has the potential to detect dozens to hundreds of additional long lags. Detecting these long lags can teach us about the vertical structure of quasar disks and how it scales with different quasar properties.


INTRODUCTION
Quasars are among the most luminous objects in the Universe.As a result, they are important to a variety of fields within astrophysics, from the study of accretion disks to galaxy evolution (Fabian 2012;Kormendy & Ho 2013).With the exception of the extraordinary efforts by the Event Horizon Telescope Collaboration et al. (2019) and Gravity Collaboration (2018), it is generally not possible to resolve the microarcsecond scales of the accretion disk, broad line region (BLR), and event horizon.Since spatial resolution is so difficult to achieve and quasars vary over a range of frequencies, temporal resolution is instead often used to study the accretion disks of quasars.
In particular, Blandford & McKee (1982) proposed a technique known as reverberation mapping to measure distances to the BLR, by measuring time lags (of order 100 days) between the continuum emission and broad emission lines (see also Kaspi et al. 2000;Peterson et al. 2004;Bentz & Katz 2015;Grier et al. 2017).A second type of reverberation mapping, first introduced by Collier et al. (1999) and known as disk continuum reverberation mapping, uses time lags between variability in different wavebands to measure the radial extent of the disk, using the lamppost model (e.g., Sergeev et al. 2005;Cackett et al. 2007Cackett et al. , 2018a;;De Rosa et al. 2015;Edelson et al. 2015Edelson et al. , 2017Edelson et al. , 2019a;;Jiang et al. 2017;Fausnaugh et al. 2016a;Starkey et al. 2017;Homayouni et al. 2022).The lamp-post model proposes that high-frequency (Xray or hard-UV) emission from the corona will be reprocessed by the disk and re-emitted as UV/optical light.Because it takes light time to travel from the inner hotter regions of the disk, where it is re-emitted in shorter wavebands, to the outer cooler regions of the disk, where it is re-emitted in longer wavebands, there will be a time lag of order days between variability in short wavelength bands and long wavelength bands.Because this lag is on the order of days, we refer to it here as the "short" lag.
Using the Shakura & Sunyaev (1973) accretion disk model of an optically thick, geometrically thin disk to determine the temperature dependence of the disk, this time lag should scale as τ c ∝ λ 4/3 , and also depend on the mass of the supermassive black hole (SMBH) and accretion rate (see also, Lynden-Bell 1969;Pringle & Rees 1972;Novikov & Thorne 1973;Friedjung 1985).However, recent observations have suggested that the radial extent of the disk is larger than predicted by the standard thin disk model (Jiang et al. 2017;Cackett et al. 2018b;Mudd et al. 2018;Guo et al. 2022), although this discrepancy may be due to nuclear extinction (Gaskell 2017;Gaskell et al. 2023) or an underestimation of SMBH masses due to unknown BLR geometry (Pozo Nuñez et al. 2019).Other models suggest that AGN disks may appear more extended due to large azimuthal temperature fluctuations (Dexter & Quataert 2012), disk winds (Laor & Davis 2014), or magnetically elevated accretion (Dexter & Begelman 2019), none of which are included in the standard thin disk model.
There is also growing evidence that the lamp-post model may be insufficient to fully explain the optical/UV variability in quasar light curves (e.g., Arévalo et al. 2009;Shappee et al. 2014;Edelson et al. 2019b;Yu et al. 2022) Alternative sources of optical/UV light curve variability, which may occur on the thermal timescale, are magneto-rotational instability turbulence (Balbus & Hawley 1991;Burke et al. 2021) and convection due to enhanced opacity in the optical/UV region of the disk (Jiang & Blaes 2020).As these turbulent perturbations are accreted inwards on the viscous timescale, they can also modify the variability of light curves in inner regions of the disk (Lyubarskii 1997).Because these perturbations are accreted inwards the lag will be "negative" compared to the short lag.In other words, light curve variability from disk turbulence in short wavelength bands will lag variability in long wavelength bands.The viscous timescale is (Davis & Tchekhovskoy 2020), where α is the Shakura & Sunyaev (1973) stress parameter, Ω = (GM/r 3 ) 1/2 is the orbital frequency, M is the mass of the SMBH, G is the gravitational constant, h is the scale height of the disk, and r is the radial position in the disk.
In the standard thin disk model τ visc can be of order hundreds of years, because in the radiation-dominated region of the disk h is constant as a function of r, and h/r ≲ 0.01.However, a recent simulation by Jiang & Blaes (2020) suggests that due to variations in the scale height of the disk, the viscous timescale could be on the order of 100 days.Furthermore, two recent independent works have suggested the detection of negative lags in the nearby quasar Fairall 9 (Hernández Santisteban et al. 2020;Yao et al. 2022).In particular, using the SWIFT and Las Cumbres Observatory (LCO) light curves with a baseline of ∼ 270 days and sub-daily cadence from Hernández Santisteban et al. (2020), Yao et al. (2022) found a long negative lag of ∼ −70 days between the W2 -and z -band.This sub-100 day negative lag supports the conclusions in the simulation in Jiang & Blaes (2020) and suggests that it may be possible to detect long lags in many more quasars.
In this paper, we derive a new empirical model for the "long" negative lag based on the model fit to the long lags detected for Fairall 9 by Yao et al. (2022).We then use this model to examine our ability to detect long lags in the future.Detecting more long lags will help to constrain our model and provide information on the vertical structure of disks.Additional lag measurements would also allow us to model how the vertical structure of disks change as a function of quasar properties such as SMBH mass and accretion rate.
Because long lags can last tens to hundreds of days, a long baseline of high cadence observations is necessary to detect them.The Vera Rubin Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) will carry out a comprehensive survey of the southern sky and will contain over 100 million quasars.LSST will provide photometric light curves in ugrizy bands, allowing us to study the wavelength dependence of quasar lags.Kovacevic et al. (2022), Pozo Nuñez et al. (2023), and others found that LSST should be able to detect short lags in quasar light curves.Here, we generate mock LSST quasar light curves reprocessed with several different short and long lags.We then test our ability to detect long lags using several lag detection methods, primarily focusing on two: javelin (Zu et al. 2011) and a maximum-likelihood Fourier method developed by Miller et al. (2010) and adapted for use on quasars by Zoghbi et al. (2013) and Cackett et al. (2022).
We derive our model in Section 2, describe in detail how we generate our mock LSST light curves in Section 3, and describe our lag detection methods in Section 4. In Section 5 we report the accuracy of our lag detection methods for light curves reprocessed with only long lags (Section 5.1), light curves reprocessed with long and short lags (Section 5.2), light curves observed with different proposed LSST cadences (Section 5.3), light curves from the simulation in Jiang & Blaes (2020) (Section 5.4), and light curves for dimmer magnitude quasars (Section 5.5) We provide a brief summary of the accuracy of our lag detection methods for all light curves in Section 5.6.In Section 6 we make a prediction for the number of long lags detectable by LSST over its ten year baseline.We summarize our results in Section 7 and discuss the importance of improved light curve models (Section 7.1) and detection methods (Section 7.2).

WHAT IS THE LONG LAG TIMESCALE?
The long lag timescale recovered by Yao et al. (2022) is much shorter (70 days versus 100 years) than the viscous timescale predicted by the standard Shakura & Sunyaev (1973) thin-disk model.However, this timescale is in good agreement with the recent simulation by Jiang & Blaes (2020) that shows h/r > 0.01 in the radiation pressure dominated region.Disks with larger aspect ratios are also often invoked to explain line driving broad absorption line quasars (e.g.Murray et al. 1995;Leighly 2004).If this shorter timescale for the long lag is accurate, it should be possible to detect long lags with long baseline campaigns, like LSST, depending on the properties of the quasar and waveband coverage.We therefore use the results from Yao et al. (2022) to derive a new empirical model for the long lag timescale.
As in Yao et al. (2022), we start by converting the viscous timescale to a velocity, v visc = r/τ visc , and integrate from r 0 , the radius corresponding to the reference band, to r x , the radius corresponding to band x, . (2) Here we have written h ≡ h 0 (r/r 0 ) β in order to capture the r-dependence of h.In the standard thin disk model, β = 0, as h is independent of radius.Following Fausnaugh et al. (2016b) and Jiang et al. (2017) the radius at which the disk radiates at a wavelength λ is, where k B is Boltzmann's constant, ĥ is Planck's constant, c is the speed of light, σ is the Stefan-Boltzmann constant, and Ṁ is the accretion rate of the SMBH.X = 2.49 is a factor to account for emission at this wavelength from a range of temperatures.f i is a factor to account for the irradiation from an X-ray/UV source near the SMBH by changing the normalization of the tem-perature profile.We set f i = 1, because three dimensional radiation magneto-hydrodynamic (MHD) simulations show that the UV/Optical region of the disk is very optically thick, making it difficult for X-ray photons to convert to UV photons, because the free-free absorption opacity for X-rays is very small and the Compton process is not important (Jiang et al. 2016(Jiang et al. , 2019;;Jiang & Blaes 2020).As a result, the UV/Optical emission should be dominated by photons generated by the local disk (f i = 1).f i = 1 implies that X-ray light will not be directly reprocessed by absorption and re-emission in the disk.Instead in this model the X-ray short timescale variability is reprocessed by perturbing the local photosphere via scattering on the X-ray variability time scale.While Equation ( 3) is derived for the standard thin disk model, r ∝ λ 4/3 does not rely on the assumption that h/r ≪ 1, only the assumption that the emission is local.Since we care most about the wavelength dependence of the long lag timescale we use Equation (3) here.However, the mass and accretion rate dependence of the long lag may need to be corrected when a more accurate model for the long lag can be derived either through simulations or additional observations.In particular, we note that h/r and thereby β may depend on the accretion rate and mass of the quasar.For simplicity, here we will define h/r and β independently of any quasar property.
Substituting Equation (3) into Equation (2) gives a viscous lag timescale of the form, where, λ W = 1928 Å, and h W /r W is the aspect ratio at the disk radius corresponding to the W2 -band.The parameterizations and bands are based on the recent fitting that Yao et al. (2022) performed on the wavelengthdependent long lag they detected for Fairall 9 to this function for a reference band λ W .They found τ 0 = 66.8 ± 2.3 and β = 3.59 ± 1.22.For Fairall 9, M = 2.55 ± 0.56 × 10 8 M ⊙ (Peterson et al. 2004) and the Eddington ratio l/l edd = 0.02 ± 0.004 (Hernández Santisteban et al. 2020), which gives τ 0 = 66.8 days for α = 0.1 and h W /r W = 0.2.Viscous Timescale [days] Figure 1.The leftmost panel shows the rest-frame short lag between the u-and y-band in days calculated using the standard thin disk model and the remaining panels show the long lag timescale between the u-and y-band in days for disks with aspect ratios from left to right of 0.03, 0.1, and 0.3 as a function of Eddington ratio and SMBH mass in M⊙.The white point shows the approximate SMBH mass and Eddington ratio of Fairall 9 (Peterson et al. 2004;Hernández Santisteban et al. 2020).; and in the bottom right panel we use the ratios between the long and short lag for different wavebands in Yao et al. (2022).
The wavebands for LSST range from the u-band to the y-band.For λ u = 3600 Å we get, for our fiducial values.
We show the lag between the u-and y-band at zero redshift for a range of quasar masses and Eddington ratios in Figure 1.For comparison we show the short lag timescale for the standard thin disk model in the left panel using Fausnaugh et al. (2016a), where r(λ 0 ) is calculated from Equation (3).The remaining three panels are for different values of h W /r W : h W /r W = 0.03, which is similar to its proposed value in the standard thin disk model, and h W /r W = 0.1, 0.3, which are slightly smaller and larger, respectively, than the value found from fitting the lag found in Yao et al. (2022).Depending on h W /r W , at low redshift the long lag can range from a fraction of a day to dozens of years, for realistic quasar properties.
In addition to the mass and Eddington ratio of the quasar, the long lag also depends strongly on the redshift, z, of the quasar.First, time dilation due to increasing redshift should increase lag times by a factor of (1+z).Second, the rest frame waveband of a quasar will change by a factor of 1/(1+z), i.e. λ 0 will decrease by a factor of 1/(1+z).Since τ 0 scales as λ 4/3(7/2−2β) 0 for β ≈ 3.6 overall τ 0 scales as, The top left panel of Figure 2 shows the normalized distribution of τ long from the u-to y-band for quasars in the Sloan Digital Sky Survey Data Release 7 (SDSS DR7, Shen et al. 2011) for three different redshift ranges, using the β ≈ 3.6 found in Yao et al. (2022).z = 1 is the highest redshift shown for two reasons.First, at z > 1 the observed-frame u-band corresponds to rest frame bands shorter than the shortest UV band in Yao et al. (2022) and we do not want to extrapolate our models beyond the data available.Second, for z ≳ 1 the long lag is over three years long for a majority of quasars if β is indeed as large as predicted by the model fit in Yao et al. (2022).Lags longer than three years would be difficult to detect even with the ten year baseline of LSST.
In the top right panel of Figure 2, instead of using the exact fit from Yao et al. (2022) we assume β = 2.For β = 2 the lag timescale would still flatten at longer wavelengths, as in Yao et al. (2022), but the aspect ratio would only increase linearly as a function of radius.If the average value of β is closer to this value, lags in quasars out to z ≲ 2 should be measurable.
We now note a few things about the value of β.First, due to the denominator of τ 0 in Equation ( 5), we restrict β ̸ = 7/4.Second, the dependence of the long lag on reference wavelength would be very different for β < 7/4.The term in the brackets in Equation (4) would become a positive factor and τ 0 would become negative.The increase in h/r becomes less steep than the other radial dependencies of the long lag.Therefore the long lag would increase with reference wavelength as is the case with short lags (see Equation (7)).
For example, for β = 0, as in the standard thin disk model, τ long ∝ λ 14/3 0 , which is a stronger dependence on λ 0 than the τ short ∝ λ 4/3 0 for the short lag.As a result, the observed long lag from the u-to y-band would actually decrease as a function of redshift, because emission from the inner disk where the lag timescale is shorter gets redshifted into longer wavebands.The bottom left panel of Figure 2 shows the distribution of τ long for three different redshift ranges for β = 0 and h/r = 0.01, which are typical assumed values for the standard thin disk model.These distributions peak around 10 5 −10 6 years.
As an alternate way of scaling from Fairall 9 to other targets, we can assume that the ratio between short and long lags is preserved between Fairall 9 and other targets.We use this alternative way of scaling to examine how removing our various model assumptions changes the long lag distributions for our different redshift ranges.To calculate these lags we take the redshift of each quasar in SDSS DR7 with z < 1 and use that to determine the rest-frame wavebands that correspond to the observed frame u-and y-band for that quasar.We then take the ratio between the long and short lag detected in Yao et al. (2022) in the rest-frame and multiply it by the short lag (calculated using Equation ( 7)) for that quasar.In practice we use the error-weighted mean of this ratio in the rest-frame waveband corresponding to each quasar and both adjacent wavebands.
Calculated in this way the lags are shifted towards slightly shorter timescales, relative to our model with β = 3.6.This shift is especially the case for z < 0.5 quasars and 0.75 < z < 1 quasars for which the distributions peak around 20 days instead of 40 days and 250 days instead of 400 days, respectively.This shift towards shorter lags would enable the detection of long lags for light curves with shorter baselines than LSST and make it possible to detect long lags out to z > 1 with LSST.The main reason for the shift to shorter timescales is that this method removes the mass and luminosity dependence of the lag timescale, by scaling all lags to Fairall 9.In the bottom right panel of Figure 2 we show the long lag timescale calculated using this approach.

GENERATING MOCK LIGHT CURVES
With a prescription for calculating lags in hand, we turn now to creating mock multi-band LSST-like light curves incorporating this negative lag.We outline our process for generating mock LSST light curves that have been reprocessed with lags.First, we either generate a driving light curve using a DRW model or use light curves from the radiation MHD simulation in Jiang & Blaes (2020).Then, we reprocess these light curves with Because there is no short lag and it is at frequencies higher than the wrapping frequency (vertical grey line), the highest frequency bin is consistent with zero.
a Gaussian response function to add lags.Finally, we sub-sample these light curves with proposed LSST cadences.

Damped Random Walk Light Curves
For light curves analyzed in Sections 5.1, 5.2, and 5.3 we start by generating a highly sampled (∆t = 0.001 days) 2000 year-long light curve, using a damped random walk (DRW) model as in Kelly et al. (2009).A DRW is a stochastic process that goes from a ν −2 powerlaw, where ν is frequency, at high-frequency to white noise at low-frequency at some characteristic damping timescale τ damp .The DRW model has been shown to be a decent representation of quasar light curves for timescales of days to years (e.g. Kelly et al. 2009;Koz lowski et al. 2010;MacLeod et al. 2012;Zu et al. 2013).The physical interpretation of this damping timescale is still debated, although Kelly et al. (2009) and Burke et al. (2021) found that it may be related to the thermal timescale of the disk.Burke et al. (2021) found that τ damp scales with the SMBH mass as, To determine the properties of our mock light curves we take quasars from SDSS DR7, and put them in three redshift bins: z < 0.5, 0.5 < z < 0.75, and 0.75 < z < 1 (see Table 1).In this paper we only investigate our ability to detect lags in quasars at z < 1, because, as can be seen in Figure 2, a majority of quasars at z ≈ 1 have long lags over 500 days long, which would be difficult to detect with the LSST baseline of 10 years.We use the median quasar mass in each bin to determine τ damp .We assume the structure function at infinity, SF ∞ = 2.5 in flux units, which is a typical value for quasar light curves from the Dark Energy Survey (Yu et al. 2020) and describes the long term variability of the quasar.
We split our 2000 year light curve into 100 individual light curves which we use as underlying driving light curves, d(t).We generate and subsample a 2000 year light curve instead of generating 100, 20 year-long driving light curves because it is important when generating a DRW light curve that the baseline is significantly longer than the damping timescale (Koz lowski 2017; Stone et al. 2022).Done this way our baseline DRW light curve is over 3000 times longer than the largest damping timescale we use.

Simulated Light Curves
While the DRW model does a good job describing quasar variability over a range of timescales, it is not based on a physical model of the disk.The DRW model may also not model the PSD properly at timescales less than 1 month (Mushotzky et al. 2011a;Caplar et al. 2017;Smith et al. 2018;Sánchez-Sáez et al. 2018), and be unreliable on the longest timescales no matter how long the modeled DRW is (Stone et al. 2022).Therefore, in Section 5.4 we test the robustness of our results by using light curves from a simulation performed by Jiang & Blaes (2020).
Jiang & Blaes (2020) used athena ++ to perform a global radiation MHD simulation of a quasar disk between 30 r g and 100 r g , where r g is gravitational radii.This region of the disk, which roughly corresponds to the region that emits in the optical, is dominated by the Rosseland mean opacity, which increases around 1.8 × 10 5 K due to the iron opacity bump.Jiang & Blaes (2020) found that the iron opacity bump causes the disk to become convectively unstable leading to turbulence which puffs up and heats the disk.Eventually the puffing up of the disk decreases the opacity enough for the disk to become convectively stable and cool, until it becomes optically thick again.In Jiang & Blaes (2020) this cycle repeats itself every few years and leads to variability in the luminosity of the disk by factors of around 3-6.
Figure 14 in Jiang & Blaes (2020) shows the luminosity scaled by the Eddington luminosity emitted from r ≤ 60r g and r ≤ 80r g over roughly 57 years for their simulated 5 × 10 8 M ⊙ quasar.Over the 57 years the Eddington ratio can vary from l/l edd ≈ 0.01 − 1.We split both of these light curves into eight chunks of 20 years each, and linearly interpolate them down to subdaily cadence from an original sampling of roughly every three days.This interpolation should not have a large impact on our results because the LSST cadence we down-sample the light curve with has a cadence of roughly every 10-30 days (see Section 3.3).

Reprocessed Light Curves
Next we reprocess our driving light curves with a Gaussian response of the form, where C is the center of the Gaussian (i.e. the input lag) and S is the width of the Gaussian, set to S ≡ C/5, for simplicity.For a negative lag (C < 0) we replace t with −t.We choose S ≡ C/5 such that there is some width to the response that scales with the lag duration.However, the exact shape of the response function is unknown, especially for the long lag.
In cases where we are simulating both a long and a short lag we add two response functions together.We scale the relative strength of the two response functions by fixing the amplitude of the long lag response function to some fraction of the amplitude of the short lag response function.As a result, the integrated Gaussian for the long lag will no longer sum to one, which means the amplitude of the reprocessed light curve will differ from the amplitude of the driving light curve.
Because we want to examine a broad range of long lag timescales, we use long lags of −50, −130, and −400 days.From Figure 2 we can see this range of timescales does a good job representing quasars from SDSS DR7 at z ≲ 1.To match long lag timescales to other quasar light curve properties we estimate that each quasar bin z < 0.5, 0.5 < z < 0.75, and 0.75 < z < 1 corresponds to a long lag of −50, −130, and −400 days, respectively.We also calculate the short lag corresponding to each bin using Equation ( 7).The long and short lags for each redshift bin are in Table 1.
To get the reprocessed light curve we Fourier transform the driving light curve (d(t) → D(ν)) and response function (ψ(t) → Ψ(ν)).In a linear reprocessing model the reprocessed light curve is, which in Fourier space is equivalent to, We can then inverse Fourier transform R(ν) → r(t) to acquire our reprocessed light curve.

LSST Cadence
Now that we have a driving and reprocessed light curve, we down-sample their ∆t = 0.001 cadence to a variety of LSST cadence possibilities.The LSST cadence is not final.Over the past decade the LSST community has conducted a thorough investigation of prospective cadences (Connolly et al. 2014).Different proposed cadences are available through the LSST Operations Simulator (OpSim, Reuter et al. 2016) 1 .OpSim simulates the field selection and image acquisition for LSST observations over its proposed 10 year baseline.It takes into account sky brightness, bad weather predicted from localized weather models, seeing, telescope maintenance times, slew times, and filter changes.
The LSST observing strategy includes two main components, the Wide Fast Deep Survey and the Deep Drilling Fields (DDFs).The Wide Fast Deep Survey will be at least 18000 deg 2 where each 9.6 deg 2 field has a median of 825 visits summed over all six filters.The DDFs are five relatively small (9.6 deg 2 ) fields which will have deeper coverage and a higher cadence than the Wide Fast Deep Survey.In this paper, we focus on objects observed in the DDFs.Specifically, here we use the proposed cadences for the XMM-LSS DDF.Other DDFs have similar cadences.
In this work, we use the baseline v2.0 DDF cadence, and the DDF cadence with the longest available season length (ddf season length slf0.10 v2), which uses 80% of the available season, or a total of roughly 210 days.We primary use the u-and y-band cadences, to take advantage of the largest wavelength range available for LSST.Over ten years the baseline v2.0 cadence has 1156 and 4515 samples for the u-and y-bands, respectively, while the long season cadence has 842 and 4434 samples for the u-and y-bands, respectively.We down-sample this number to eliminate visits that occur within 1/100 of a day of each other and end up with 186 (150) and 443 (443) samples for the baseline (long season) cadence, for the u-and y-bands, respectively.This down-sampling still leaves a large fraction (15-35% depending on the band and cadence) of visits that occur within roughly 15 minutes of each other.For the baseline cadence, most of the remaining samples are separated by roughly 1 to 3 days.For the long season cadence the separation is typically around 10 to 30 days, with more visit separations on the longer side for the u-band.
We find a longer season length, and as a result shorter season gaps, to be beneficial for detecting long lags in our mock light curves, even though the longer season has lower cadence during the season (see Figure 6 and discussion in Sections 4.2 and 5.3).Therefore, in Sections 5.1, 5.2, 5.4, and 5.5 we use the long season cadence.In Section 5.3, we discuss the relative performance of our lag detection methods on light curves with the baseline cadence and ways to improve our results if the baseline cadence is selected.
After down-sampling our driving light curve, we add band-dependent dust extinction and measurement error by calculating the signal to noise ratio (SNR) for each observation using the sky brightness and instrument noise from OpSim.In this calculation we use the average magnitudes in each SDSS quasar bin for quasars with i < 20 mag (see Table 1), to first test our ability to detect the long lag in the most luminous quasars.The median SNR for our mock quasars over all bins is ∼ 150.We show that we should be able to detect long lags for fainter quasars as well in Section 5.5.
The top panel of Figure 3 shows an example DRW light curve.The driving light curve is shown in blue, and the same light curve reprocessed with a long lag of −50 days is shown in orange.The middle panel shows the same light curves sub-sampled with the long season LSST cadence for the u-band and y-band, for the driving and reprocessed light curves, respectively.

LAG DETECTION METHODS
In this section, we outline several methods developed for detecting short lags and BLR lags that we will apply to the detection of long lags.First, we briefly discuss different cross-correlation methods (Gaskell & Peterson 1987;White & Peterson 1994;Peterson et al. 2004), and the Von-Neumann estimator presented by Chelouche et al. (2017).Next, we give a more thorough description of the two methods we explore in detail in this paper, javelin (Zu et al. 2011) and a Fourier method (Zoghbi et al. 2013).We test a variety of methods on our mock light curves because each one has different strengths and drawbacks.These drawbacks mostly arise because these methods need to work for unevenly sampled UV/optical light curves, like those LSST will produce, which will not have an even cadence because of day-time gaps, bad weather, seasonal gaps, shared telescope time, etc.

Cross-Correlation Methods and the Von-Neumann Estimator
With uniformly-sampled light curves, time lags can be estimated by determining the delay time, τ , with the largest cross-correlation coefficient, Multiple lags can in principle be detected by finding multiple peaks in the correlation coefficients.We show the cross-correlation coefficients as a function of lag for mock DRW light curves with an evenlysampled, daily cadence and ten-year baseline in Figure 4.The light curves have been reprocessed with long (short) lags of −50 (7), −130 (9), and −400 (12) days, in the top, center, and bottom panel, respectively.As the ratio between the amplitude of the long lag and the amplitude of the short lag decreases the peak around the long lag decreases.When the amplitude ratio is 0.01, the peak around the long lag nearly entirely disappears for input long lags of −50 and −400 days, although there remains a small peak around the long lag for the light curve with a −130 day input long lag.
On the other hand, as the amplitude ratio decreases, the peak around the short lag becomes more pronounced.For an input long lag of −50 days it becomes the largest peak when the amplitude ratio is 0.1.For the other two input long lags, the peak around the short lag becomes dominant when the amplitude ratio is 0.01.We show the positive and negative lag values with the highest correlation coefficients for these evenly-sampled light curves in the top four rows of Table 2.
To use this cross-correlation technique on unevenly sampled data, such as LSST light curves, we use the interpolated cross-correlation method (ICCF, Gaskell & cross-correlation on evenly sampled light curves.That is, the uneven cadence does not prevent us from being able to find an accurate lag, but as the amplitude ratio between the long and short lag decreases we are less (more) able to recover the input long (short) lag.We do not try other cross-correlation methods on our light curves, such as the z-transformed discrete correlation function (ZDCF, Alexander 2013), but we expect these methods to perform similarly to the ICCF.We do however try the Von-Neumann estimator presented by Chelouche et al. (2017).This Von-Neumann method shifts two light curves by some time τ and measures the mean-square successive-difference, which is a measure of the randomness, between the two light curves as a function of τ .The lag is therefore the τ with the smallest measure of randomness.This method is unique because it involves no model fitting, interpolations, or binning, and is therefore ideal for unevenly sampled light curves.
Figure 5 shows the mean of the mean-square successive differences for 1000 Monte Carlo Markov Chains (MCMCs) for an example mock light curve reprocessed with a long (short) lag of −50 (7) days in the top panel, −130 (9) days in the center panel, and −400 (12) days in the bottom panel.The ratio of the amplitude of the long lag to the amplitude of the short lag is 0.2, 0.1, and 0.01, in the left, center, and right panels, respectively.Because the Von-Neumann method assumes a lag that is constant with frequency it has trouble identifying two lags simultaneously.For the shortest duration, −50 day, lag this assumption means the Von-Neumann method combines the two lags into a single wider trough when the amplitude ratio is 0.2.That trough gets skewed towards the short lag when the amplitude ratio is 0.1.For the longer duration lags the assumption of only one lag results in the Von-Neumann method finding only one of either the long or the short lag for most amplitude ratios.Only for the light curve with an input long lag of −130 days and an amplitude ratio of 0.1 are both lags detected at all, with significant dips in the distribution near each lag.As with the cross-correlation method, the single-lag assumption is a drawback of using the Von-Neumann method to detect a multi-lag signal, unless one signal can be filtered out in some way.However, for an amplitude ratio of 0.2, we find that the Von-Neumann method could be a helpful tool to detect long lags in LSST light curves with long seasonal gaps.We discuss this further in Section 5.3.

Javelin
Our first main lag detection method is javelin2 , previously known as spear, which is a stochastic process estimator used to look for time lags in spectroscopic (Zu et al. 2011) or photometric (Zu et al. 2016) light curves.Li et al. (2019) performed a comprehensive comparison of the performance of javelin versus cross-correlation methods over a broad range of light-curve properties for mock SDSS light curves and found that javelin performs best at accurately detecting BLR lags.As our lags are of a similar timescale we use javelin as one of two main lag detection methods here.
javelin works on unevenly sampled light curves by fitting a model to the available data and running an MCMC to find the posterior distributions for the DRW parameters (τ damp and SF ∞ ) of each light curve.javelin then uses these posteriors to statistically interpolate each light curve as it shifts them temporally and uses an MCMC to simultaneously derive the posterior distributions for the DRW parameters and the response function between different light curves.For its default settings, which we use here, javelin models the response function as a top hat with a finite width centered at time lag τ .
The DRW statistical interpolation works well for shorter gaps in unevenly sampled light curves but is unable to fit the large seasonal gaps in our data, which becomes a problem for the LSST baseline cadence.Figure 6 shows a javelin distribution for an example mock baseline cadence light curve.The largest peaks are around ±180 days, which is roughly the length of the gaps in between each observing season.We discuss solutions to this problem in case the baseline cadence is chosen in Section 5.3.
The bottom left panel of Figure 3 shows the javelin distribution for the example light curves in the panels above, which have been sub-sampled with the long season cadence and reprocessed with an input long lag of −50 days.The probability distribution peaks very sharply around −50 days and the median of the negative values of the distribution is −50 days.
Because the long lag is negative and the short lag is positive, we take the long (short) lag as the median of the negative (positive) values of the javelin distribution.This method could bias our values for the short lag towards slightly larger values, because the short lags are near zero and so the peak related to the short lag could extend below zero.However, we find that peaks in javelin distributions around the input short lags rarely extend past zero.
While javelin is a useful method for detecting lags in quasar light curves it does make assumptions about the light curve that are model dependent.Also, because the second light curve is modeled as a smoothed and shifted version of the first, the lag can depend on which light curve you model first.In addition, because javelin uses a DRW model, which is what we use to generate our mock light curves, we could be biasing our comparison to favor javelin.We address this potential bias in Section 5.4, where we test our methods on light curves from radiation MHD simulations instead of DRW mock light curves.

Maximum-Likelihood Fourier Method
Our second main lag detection method is a Fourier method.Fourier methods are useful because they allow us to separate lags as a function of frequency.They are straightforward as well because calculating the cross-spectrum is a simple multiplication.We can use Equation ( 12) and the power-spectral densities (PSD) of the reprocessed and driving light curve, |D(ν)| 2 = D * (ν)D(ν) and |R(ν)| 2 = R * (ν)R(ν), respectively, to derive the cross-spectrum of the two light curves to be,  The phase lag, ϕ, is then just the phase of the complex cross-spectrum, which is related to the time lag as (Uttley et al. 2014), We show the lags as a function of frequency recovered by this Fourier method between two 10 year baseline, evenly-sampled, daily cadence light curves in Figure 7.The input long (short) lags are −50 (7) and −130 (9) days in the top and bottom panel, respectively.We do not use this method to detect long lags of −400 days because the LSST baseline is not long enough to recover enough data points at frequencies < 10 −3 days −1 .The vertical grey lines show the phase wrapping frequencies, which are equivalent to ν w = 1/2C for a symmetric response function (Uttley et al. 2014).Phase wrapping, the sudden flip in the sign of the lag, occurs because phase lags are limited to between ±π, so a positive or negative shift of half a wavelength cannot be distinguished.
For light curves with only an input long lag, at frequencies lower than the first phase wrapping the measured lag is −51 and −131 days for input long lags of −50 and −130 days, respectively.At frequencies above the long lag's phase wrapping, the lag values oscillate around and then approach zero.When a short lag is added and the ratio between the amplitude of the long lag and the amplitude of the short lag is 0.2, the lags measured at low frequencies become roughly −27 and −95 days for input long lags of −50 and −130 days, respectively.For an amplitude ratio of 0.1 the duration of these lags decrease again to −17 and −74 days for input long lags of −50 and −130 days, respectively.For an amplitude ratio of 0.01 the lag measured at the lowest frequencies becomes positive for an input long lag of −50 days and is only −8 days for an input long lag of −130 days.
The decrease in the duration of the recovered long lags as the amplitude ratio decreases is due to the growing influence of the short lag on lower frequencies.While the short lag has the strongest impact at frequencies around the inverse of the short lag timescale, it is also applied to all frequencies below that.Therefore it reduces the value of the long lag.
On the other hand, because the long lag only acts at lower frequencies, the short lag is easier to isolate in Figure 7.For light curves with two input lags, at frequencies higher than the long lag phase wrapping frequency the lag values oscillate before quickly converging on the short lag value for all amplitude ratios.The smaller the amplitude ratio the sooner they converge.The lag values stay at the input short lag, until the short lag wrapping frequency is reached, when the lag values again oscillate and then in this case approach zero.We summarize the long and short lag values in Figure 7 in the top four rows of Table 2.This Fourier method is very useful for recovering the shape of a response function for evenly sampled X-ray data (e.g.Papadakis et al. 2001;McHardy et al. 2007;Arévalo et al. 2008;Fabian et al. 2009;Zoghbi et al. 2010;De Marco et al. 2013;Cackett et al. 2013;Uttley et al. 2014), but requires a work-around for unevenly sampled data, which can not be fast Fourier transformed.
Here we use a maximum-likelihood method described in Zoghbi et al. (2013) and first presented by Miller et al. (2010).
First, we create N = 4 bins in frequency space ranging from f min = 1/T , where T = 10 years is the length of the light curve, to f max = 0.055 day −1 .As in Zoghbi et al. (2013) we also add a bin on either end to be discarded later to reduce bias in the minimum and maximum bins.f max is chosen based on the maximum frequency for which there is significant power for an LSST cadence.Figure 8 shows the distribution of the inverse of the time between each observation for the LSST long season cadence.Because the length of the light curves is 10 years, frequencies less than 10 −2 days −1 are common, but because the cadence is roughly every ten days, frequencies greater than a few times 10 −2 days −1 are rarer.
It is easier to obtain accurate results with fewer bins.Previous testing shows that errors for each bin are smaller when the width of bins are at least a few times 1/T .This requirement makes it difficult to detect −400 day lags, and so we do not attempt to detect that lag with this method.Otherwise, four bins are sufficient for our purpose here: detecting lags from two Gaussian response functions representing the long and short lag, which occur at order of magnitude different frequencies.However, additional bins would be useful for detecting the shape of more complicated response functions.
Using these bins we fit a piece-wise PSD to each light curve and calculate the likelihood using the autocovariance function, where |D(ν)| 2 = i D i is the PSD as a piecewise sum over a range of frequency bins.We then maximize this likelihood to find the PSD for each light curve.We use these PSD as input parameters for an MCMC, which gives us a distribution for the power in each bin.We take the median of the distribution for each bin as the power for that bin and combine them to form a piecewise PSD over the full range of frequencies.Next, we maximize the log-likelihood of the crosscovariance, to obtain the best fit values of Ψ(ν) = i Ψ i and ϕ(ν) = i ϕ i as a function of frequency using the PSD from the MCMC distributions.We again use the best fit parameters as input for an MCMC, which gives us the distribution of ϕ(ν).For each bin we use the median of the distribution as that bin's ϕ(ν) and the standard deviation as the error of ϕ(ν).Finally, we can use ϕ(ν) to determine the lag for each frequency bin using Equation ( 15).
We show an example of the lag as a function of frequency detected using this maximum-likelihood method in the bottom right panel of Figure 3, for the same light curve with an input long lag of −50 days shown above.The vertical grey line shows the phase wrapping frequency for −50 days.In Figure 3 the lowest frequency bin has the largest error-bar and is consistent with zero.In general for our mock light curves this bin tends to have the largest error and to be the least consistent with the input long lag.For these reasons we neglect the lag found in this bin, and then calculate the long lag as the error-weighted mean of all recovered lags < −10 days.We apply a cut-off of −10 days, because lags will approach zero at the wrapping frequency, and we only care about negative lag values when searching for the long lag.In Figure 3 the two bins below the wrapping frequency have an error-weighted mean value of −49 days, which is similar to the input lag.The bin above the wrapping frequency is consistent with zero for this light curve, because there is no input short lag.For light curves with input short lags we also use the maximum-likelihood method to detect the short lag.To avoid interference from a long lag signal, we follow Yao et al. (2022) and calculate a short lag by taking the lag value from the one bin centered between 0.01 < ν < 0.1.

Detecting Long Lags Alone
First we examine the performance of javelin and the maximum-likelihood method on light curves with the long season LSST cadence (see Section 3.3) that have been reprocessed with a long lag only.Figure 9 shows the distributions of the medians of the negative values  Table 2. Above divider: The peaks in the cross-correlation coefficient and the long and short lags detected by the Fourier method for the even-cadence light curves with different input lags and amplitude ratios (between the long and short lag) shown in Figures 4 and 7. Below divider: The median long and short lags detected by javelin and the maximum-likelihood method for 100 mock light curves with each cadence, amplitude ratio, and input long and short lag used throughout the paper.The error ranges represent one standard deviation.to detect with the maximum-likelihood method for light curves with the LSST baseline of ten years.For input lags of −50 and −130 days the maximum-likelihood method is less accurate than javelin.However, the maximum-likelihood method is able to detect the long lag to within 30% of the input lag for around 70% of light curves.We show the medians and standard deviations of the javelin and maximum-likelihood method distributions for these mock light curves in row (5) of Table 2.

Detecting Short and Long Lags Together
We now add a short lag to our mock light curves with the long season LSST cadence to examine how multiple lag signals can affect our ability to detect either lag.Our fiducial value for the ratio of the amplitude of the long lag to the amplitude of the short lag is 0.2, as we expect the short lag signal to be stronger.We also examine our ability to detect lags when this ratio is 0.1 and 0.01.For readability, in this section we only show figures for light curves with an input long lag of −50 days and an input short lag of 7 days.Unless otherwise noted, the input lag duration does not have an impact on the accuracy of our lag detection methods.We repeat the analysis in this section for light curves reprocessed with long (short) lags of −130 (9) and −400 (12) days in Appendix A.1.
In the top panel of Figure 11 we compare the javelin results for light curves reprocessed only with a long lag (in orange) to javelin results for light curves reprocessed with a long and short lag (in blue).For light curves with both lags we take the median of the negative (positive) values in the javelin distribution as the long (short) lag.When the amplitude ratio between the long and short lag is 0.2 or 0.1 the distributions of long lags detected by javelin are consistent with the distribution of long lags detected for light curves without a short lag.Adding in a short lag only makes the distributions somewhat broader.Comparing these results to the correlation coefficients for evenly-sampled daily cadence light curves, shown in Figure 4, the long lag tends to be more prominent than we might expect, in particular for an amplitude ratio of 0.1.It may be that the roughly every 10-30 day cadence of the long season LSST cadence, helps us to detect the long lag, because it hides some of the signal from the short lag.
On the other hand, when the ratio is decreased to 0.01, we fail to detect the correct long lag for nearly all mock light curves with an input long lag of −50 days.In Figure 23 in Appendix A.1, for an input long lag of −130 days there is still a peak in the distribution around −130 days.It is unclear whether this peak is tied to the cadence of the light curve or the injected long lag.For one, there is also a peak around −130 days in the lag distributions for input long lags of −50 and −400 days when the amplitude ratio is 0.01.There is also some evidence for a peak around −130 days in our null tests using light curves that have not been reprocessed with a lag (see top left panel of Figure 33  In the top panel we zoom in on the portion of the distribution between 0 and 20 days to better show the accuracy of the short lag when detected.The input short lag, shown as the vertical dashed line, is 7 days.
cross-correlation coefficient for the evenly-sampled light curve with a −130 day input long lag and an amplitude ratio of 0.01, while there is no corresponding peak in the correlation coefficient for other input lags when the amplitude ratios are 0.01.Therefore it is unclear if the peak in the distribution in the right middle panel in Figure 23 is due to a recovered lag or not.
We show the medians of the positive values of the javelin distribution for 100 mock light curves with input lags of −50 and 7 days in the top panel of Figure 12.We zoom in to the portion of the distribution falling between 0 and 20 days to better show the accuracy of the short lag when it is detected.When the ratio of the amplitude of the long lag to the amplitude of the short lag is 0.01 javelin is able to detect the short lag, although the distribution is skewed towards larger values than the input short lag.In previous studies, javelin has detected short lags that appear to be skewed towards longer frequencies due to diffuse BLR emission (e.g., Cackett et al. 2018a).Something similar may be happening here due to the presence of both a long and short lag.
For an amplitude ratio of 0.1, javelin detects a short lag of around 13-16 days for roughly 60% of light curves with input lags of −50 and 7 days.The null tests do not show a peak around 13 days for this cadence (see top right panel of Figure 33 in Appendix B.1).In addition, for the even cadence light curve with an input long lag of −50 days shown in the top panel of Figure 4, the peak in the correlation coefficient around the short lag is larger than the peak around the long lag when the amplitude ratio is 0.1.The fact that this peak around the short lag is already larger at this amplitude ratio suggests that it might be easier to detect the short lag with javelin in this particular case.Therefore it is possible that the light blue peak in the top panel of Figure 12 is due to the input short lag, but heavily skewed towards longer lags because of the presence of the long lag.For other input long lags javelin is unable to detect a short lag for most light curves unless the amplitude ratio is 0.01 (see Figure 25 in Appendix A.1.) The bottom panel of Figure 11 is the same as the top panel, but instead shows the long lag detected by the maximum-likelihood method.As with the even cadence light curves in Figure 7, when a short lag is added, the long lag is shifted towards shorter timescales.The median of the distribution for light curves with an input long lag of −50 days goes from −49 days for light curves with no short lag, to −25 days and −20 days for light curves with amplitude ratios of 0.2 and 0.1, respectively.When the amplitude ratio is 0.01, for light curves with input long lags of −50 days there are very few negative lags detected at all and for light curves with an input long lag of −130 days there are only about 20 negative lags detected, mostly between −10 and −65 days (see Figure 24 in Appendix A).It is clear that like javelin the maximum-likelihood method is unable to isolate the long lag when the amplitude ratio is 0.01.
The bottom panel of Figure 12 shows the distribution of short lags detected by the maximum-likelihood method.
Unlike javelin, the maximum-likelihood method is able to detect short lags that are accurate to within a few days of the input lag for over 60% of light curves with amplitude ratios of 0.2 and 0.1.Detecting short lags more accurately for a wider range of amplitude ratios is therefore one benefit to using the maximum-likelihood method in addition to javelin.
The accuracy of the maximum-likelihood method on light curves where we include two input lags appears to be consistent with Figure 7, where we use the Fourier method on uniform 1-day cadence, 10 year-long, DRW light curves.It is promising that the maximumlikelihood method is able to perform nearly as well on light curves with the long season LSST cadence as the Fourier method performs on light curves with an even cadence.However, it will be difficult with this method to disentangle the length of the long lag versus the amplitude ratio between the long and short lag.We show the medians and standard deviations of the distributions of lags recovered by javelin and the maximum-likelihood method for our 100 mock light curves in rows ( 6), (7), and (8) in Table 2 for light curves with amplitude ratios between their input long and short lags of 0.2, 0.1, and 0.01, respectively.In Figure 13 we show the fractional error of the medians of the distributions of lags recovered by javelin and the maximum-likelihood method for the 100 mock light curves as a function of these amplitude ratios for all input long lags.The error-bars show the fractional error of one standard deviation for each distribution.The accuracy of the recovered lags is similar regardless of the durations of the input lags.

Improving the Baseline Cadence
The LSST collaboration is still evaluating different possible cadences.We have evaluated and compared the ability of our lag detection methods to accurately pick out long lags for several cadences.Above we asserted that the most significant factor for improving our ability to detect lags is increased season length.In the previous sections we use the long season LSST cadence.In this section we compare the detected lags for light curves with this long season cadence to light curves with the baseline v2.0 cadence.We again only show results here for light curves with input lags of −50 and 7 days.We show results for all three lag bins in Table 1 in Appendix A.2. Unless otherwise noted, our results for light curves reprocessed with a long lag of −50 days are consistent with results for light curves reprocessed with other input long lags.The top panel of Figure 14 shows the long lags detected by javelin for 100 mock light curves sub-sampled with the long season and baseline cadences, in orange and green, respectively.Unfortunately, when using the baseline cadence javelin fails to accurately recover any long lags.Instead, javelin primarily recovers long lags of around −180 days.The second panel from the top in Figure 33  panel, respectively.The positive and negative distributions peak strongly around ±180 days, like the distribution shown in green in top panel of Figure 14, suggesting that the peak around ±180 days is due to the cadence, not a lag.The maximum-likelihood method, on the other hand, still performs reasonably well on the light curves subsampled with the baseline cadence.
Because the maximum-likelihood method separates the signal out in frequency space it is easier to separate out the lag signal from the cadence signal.The bottom panel of Figure 14 shows the distribution of lags recovered by the maximum-likelihood method for light curves with the long season or baseline cadence in orange and green, respectively.The distribution of long lags recovered by the maximum-likelihood method is slightly less skewed towards shorter timescales, with a median of −31 days.In Figure 28 in Appendix A the distribution of long lags for light curves with an input long lag of −130 days is more skewed towards shorter timescales, with a median of −58 days.It could be that the presence of a gap at longer timescales actually helps to skew the long lag back towards longer timescales for light curves with input lags significantly shorter than the observing gaps (∼ −50 days), while longer lags (∼ −130 days) closer to the gap length are affected differently.
In order to improve the performance of javelin on light curves with large seasonal gaps, we can add observations from other telescopes to fill in the 180 day gap.Here, we explore the effectiveness of filling these gaps with observations from LCO. LCO has been used for numerous RM campaigns and has wavebands ranging from the u-y-band.Here we experiment with adding observations in the g-and z -bands to each seasonal gap on roughly a ten day cadence.We use the g-and z -bands instead of the u-and y-bands to increase the SNR.We determine the additional LCO cadence in the gaps using the LCO g-and z -band cadences from the Hernández Santisteban et al. ( 2020) Fairall 9 observing campaign, which are roughly daily.We randomize these Fairall 9 observing times over a uniform distribution with a range of half the distance to the following observation, and then down-sample to every tenth observation.
For each lag bin in Table 1, we make 100 mock light curves sub-sampled at times corresponding to the LSST baseline cadence for the g-and z -band and these additional LCO observations.The LCO observations add roughly 18 observations per year or 180 total.We add measurement errors corresponding to the SNR for LSST and LCO.The SNR for our mock LCO observations were calculated using the LCO Exposure Time Calculator 3 and range from 19 to 77, depending on the band and redshift bin.
We show the long lags recovered by javelin for 100 of these mock light curves with input lags of −50 and 7 days in pink in the top panel of Figure 15.For comparison we show the lags recovered by javelin for 100 mock light curves with the same input lags and the u-and y-band LSST baseline cadences (in green) and g-and z -band LSST baseline cadences (in light blue).Light curves with the g-and z -band LSST baseline cadences experience the same cadence-related issues as light curves with the u-and y-band LSST baseline cadences (see also Appendix B.1).However, adding LCO observations eliminates most of the false lag detections from the large gaps in the baseline cadence.javelin recovers a long lag within 10% of the input lag for 74% (100% are within 30%) of light curves with input long lags of −50 days.The decline in false lag detections 3 https://exposure-time-calculator.lco.global/ with this LSST+LCO cadence can also be seen in our null tests in Appendix B.1.In the bottom panel of Figure 33, the distributions of lags detected by javelin for light curves with this LSST+LCO cadence and no input lag peak around zero, as expected.
The bottom panel of Figure 15 is the same as the top panel except for the maximum-likelihood method.The distribution of long lags detected for light curves with the LSST+LCO cadence is similar to the distribution for light curves with the long season cadence.The median long lag detected for light curves with the LSST+LCO cadence is −22 days for an input long lag of −50 days.This median is similar to, although skewed slightly shorter than, the median found for light curves with the long season cadence (−25 days).
We show the distributions of short lags detected by the maximum-likelihood method for the same light curves as above in Figure 16.For all cadences the distributions peak around the input short lag.However, for the baseline cadences only around 50% of short lags detected are within 2 days of the input short lag, whereas when the LCO observations are added around 80% of the short lags detected are within 2 days of the input short lag, which is slightly more accurate than with the long season cadence.For both javelin and the maximum-likelihood method it is clear that adding just 18 LCO observations a year in the gaps of the baseline cadence should be sufficient to erase the gap issue and detect lags at rates consistent with those for the long season cadence.We show the medians and standard deviations of the distributions of lags recovered by javelin and the maximumlikelihood method for light curves with the baseline uand y-band cadences, baseline g-and z -band cadences, and baseline+LCO cadence in rows ( 9), (10), and (11) of Table 2, respectively.
Finally, we find that the Von-Neumann method (see Section 4.1) is more effective at detecting long lags for mock light curves with the baseline LSST cadence than javelin or the Fourier method.The Von-Neumann method is more accurate than these other methods in this case because it does not involve any interpolation or binning.
Figure 17 shows the distributions of long lags recovered by the Von-Neumann method for 100 mock light curves with input lags of −50 and 7 days and either the long season cadence (in orange) or the baseline v2.0 cadence for the u-and y-band (in green).Unlike with javelin, the Von-Neumann method is still able to detect an accurate long lag for light curves with the baseline cadence.The median long lag for both distributions are −40 +8 respectively.Therefore, the Von-Neumann method appears to be a very useful tool for recovering the long lag if the baseline cadence is selected.

Beyond DRW Light Curves
Although the DRW model has been shown to be a good representation of quasar light curves (e.g., MacLeod et al. 2012), the DRW model is not physically motivated and may be unreliable at the smallest and largest scales.Furthermore, since javelin uses a DRW model to interpolate unevenly-sampled light curves, it is prudent to investigate its performance on light curves that are not generated using a DRW.Therefore, in this Section we use the light curves from the radiation MHD simulation in Jiang & Blaes (2020) (see Section 3.1.2).These mock light curves have the long season LSST cadence and the fiducial amplitude ratio of 0.2.The simulation in Jiang & Blaes ( 2020) is for a 5 × 10 8 M ⊙ SMBH and the light curves vary in Eddington ratio from l/l edd ≈ 0.01 − 1.A DRW can be fit to the light curves with τ damp = 200, although this fit is not unique.
Figure 18 shows the distribution of the median of the negative values of the javelin distribution for the 16 mock light curves made using the simulated light curves from Jiang & Blaes (2020).For an input long lag of −50 and −130 days all medians fall within ∼ 10% of the input lag.14 out of 16 lags detected for light curves with input long lags of −400 days are within 20% of the input long lag, with only two outliers.Both of these outliers are from light curves made using only the inner radius light curve (see Section 3.1.2),and may be from a portion of the light curve with lower variability, making it more difficult to detect the lag.Nonetheless, despite using a DRW model to interpolate light curves, javelin still performs well on light curves that are not generated as a DRW.
We show the lags detected by the maximum-likelihood method for the 16 mock light curves with input long lags of −50 and −130 days in Figure 19.We again see that the distribution of long lags is skewed towards shorter timescales as was the case for light curves reprocessed with a long and short lag in Section 5.2.The shifts in these distributions are mildly larger.The maximumlikelihood method finds long lags of −20 days (40% of the input lag) versus −25 days (50% of the input lag) and −95 days (73% of the input lag) versus −98 days (75% of the input lag) for light curves with input long lags of −50 and −130, respectively.
As with DRW mock light curves where the ratio of the amplitude of the long lag to the amplitude of the short lag is 0.2, javelin is unable to accurately detect the short lag.The maximum-likelihood method, however, is able to detect the short lag to within 50% accuracy for all light curves and to within a day for a majority of light curves.
We show the medians and standard deviations of the distributions of lags recovered by javelin and the maximum-likelihood method for these non-DRW modeled light curves in row ( 12) of Table 2. Using light curves from the radiation MHD simulation in Jiang & Blaes (2020) instead of DRW light curves appears to have little affect on our ability to detect long and short lags.Because we derive similar success rates for more physically motivated light curves, we are hopeful that our results for light curves generated with a DRW model will hold for future observed LSST light curves.

Dimmer Magnitude Quasars
In the previous sections we test our lag recovery methods on light curves with i < 19 mag, which have SNR of around 150.To test if we can detect lags at lower SNR, we generate mock long season cadence LSST light curves with the same input long and short lags as in the previous sections and a ratio of 0.2 between the long and short lag amplitudes, but set the band magnitudes to i = 21 for all input lags.For LSST this band magnitude gives an SNR of around 30.
The colored lines in Figures 21 and 22 show the distributions of lags recovered by javelin and the maximumlikelihood method, respectively, for these mock light curves.These distributions are very similar to the distributions, shown in gray, of lags recovered by both of these methods for light curves with the band magnitudes in Table 1.The median long lags recovered by javelin are −50, −130, and −410 days for input long lags of −50, −130, and −400 days, respectively.
The median long lags recovered by the maximumlikelihood method for these light curves, −26 and −98 days for input long lags of −50 and −130 days, respectively, are also very similar to the median long lags recovered for light curves with brighter magnitudes.We show these median values, their standard deviations, and these values for the short lags in the bottom row of Table 2.Because the results for these dimmer magnitude quasars are nearly identical to the results in Section 5.2, we can expect to be able to detect long lags for quasars brighter than i = 21 mag, and suggest it might be possible to look for lags for even fainter quasars.The colored lines are the distributions for light curves where we set i = 21 to explore how a lower SNR changes our ability to detect lags.The distributions of these lags are very similar to the distributions shown in gray which are for light curves with the i-band magnitudes in Table 1.

Summary of Results
We now provide a summary of Section 5. First, in Section 5.1 we find that javelin and the maximumlikelihood method both detect long lags with high accuracy in mock light curves with the LSST long season to explore how a lower SNR changes our ability to detect lags.The distributions of these lags are very similar to the distributions shown in gray which are for light curves with the i-band magnitudes in Table 1.
cadence that only have an input long lag.When we also reprocess these mock light curves with a short lag in Section 5.2, javelin is still able to accurately recover the long lag as long as the ratio of the amplitude of the long lag to the amplitude of the short lag is 0.2 or 0.1.The long lags recovered by the maximum-likelihood method get skewed to shorter and shorter durations as this amplitude ratio decreases, a trend we also see in our tests of the Fourier method on even cadence light curves (see Figure 7).This trend could potentially help us to determine the amplitude ratio between the long and short lag in cases where we can make a robust long lag prediction from an alternative method.Both detection methods fail to recover the long lag for most light curves when the amplitude ratio is 0.01, although javelin may in some cases be able to detect a long lag for this ratio when the input long lag is −130 days.javelin is unable to detect a short lag for an amplitude ratio of 0.2 for most light curves.For an amplitude ratio of 0.1, javelin is only able to detect the short lag for a significant amount of light curves for light curves with input lags of −50 and 7 days.When the amplitude ratio is 0.01 javelin is able to detect a short lag for a majority of light curves, although the short lags detected are often skewed towards longer durations.The maximum-likelihood method, on the other hand, is always able to detect the short lag, which is a real strength of the method.
In Section 5.3 we find that javelin performs poorly on light curves with the LSST baseline v2.0 cadence, identifying the season gaps instead of the long or short lag.The maximum-likelihood method and the Von-Neumann method perform better on light curves with the baseline cadence, seeming to be only minorly affected by the change in cadence.We show that adding just 18 observations a year from LCO in the gaps in the baseline cadence fixes the cadence-issues for javelin, resulting in highly accurate lag detections.Adding these additional LCO observations also improves the accuracy of the short lags recovered by the maximum-likelihood method.
In Section 5.4 we find that both javelin and the maximum-likelihood method detect lags with roughly the same level of accuracy for the light curves from the radiation MHD simulation in Jiang & Blaes (2020) as for our DRW light curves.That our methods perform this consistently on more physically motivated light curves suggests that the results in this section should be applicable to real LSST light curves in the future.Finally, in Section 5.5 we find that when we generate dimmer mock light curves with a magnitude of i = 21, which lowers their SNR, the distributions of lags detected by javelin and the maximum-likelihood method are nearly identical to the distributions for the brighter magnitude mock light curves (i < 19) in Section 5.2.

NUMBER OF DETECTABLE LONG LAGS
We now provide an estimate for the number of long lags we may expect to measure with LSST.To do this estimate we will use quasar properties from SDSS and metrics developed for the robustness of BLR and short lag detections.
First, we select only quasars at z < 1, because for our models based on Fairall 9 Figure 2 suggests that we should be able to detect the long lags for a majority of SDSS quasars at z < 1.Note that if β = 0 and h/r = 0.01, as is often assumed for the standard thin disk model, the long lag will be at least on the order of 100 years for bright SDSS quasars, and therefore not detectable using LSST.Next, we select only quasars with an i -mag brighter than 21.In Section 5.5 we show that our lag detection methods are accurate down to this magnitude.
Whether a long lag is detectable or not also depends on the variability of the quasar.When looking for BLR lags in SDSS quasar light curves, Grier et al. (2017) defined a criterion SNR2 based on the continuum and line light curve rms variability.The rms variability is the intrinsic variability of the light curve about a fitted linear trend, divided by the uncertainty of the estimated intrinsic variability.
To detect C IV lags, Grier et al. (2019) require an SNR2 > 20.They found that 348 out of 492 (∼ 71%) quasars in their sample from the SDSS reverberation mapping project that fall within 1.35 < z < 4.32 met that criterion.Homayouni et al. (2020) look for Mg II lags for 453 quasar light curves with redshifts more similar to our target redshifts here, 0.35 < z < 1.7, and find that 198, or ∼ 44%, have an SNR2 > 20.The number of quasars with SNR2 > 20 is fewer in this case because Mg II tends to be less variable.
A second criterion for determining lag significance for both BLR and short continuum lags is setting a minimum threshold for r max , the maximum ICCF correlation coefficient (see Section 4.1).Setting a minimum threshold for r max ensures that the light curves are well correlated.Guo et al. (2022) require r max > 0.8 for their sample of z < 0.8 quasar light curves from the Zwicky Transient Facility.Of 4255 quasars, 559 (13%) met that criterion.
SDSS DR14 has ∼ 5 × 10 5 quasars over 9376 deg 2 (Pâris et al. 2018).Limiting ourselves to z < 1 and i < 21 mag quasars we are left with ∼ 9.5 × 10 4 quasars.Accounting for the five planned 9.6 deg 2 DDF, we can expect to detect around 500 of these quasars within the DDFs.If we take the more conservative fraction of quasars with SNR2 > 20 from Homayouni et al. (2022) we find that around 200 quasars should have detectable lags.If we assume only 13% of light curves will have r max > 0.8 as in Guo et al. (2022) we would be able to detect long lags for 65 quasars within the DDFs.
Here, we only look at z < 1 quasars because the long lag gets very long at z > 1 and LSST has a ten year baseline.However, it is possible LSST will operate for 20 years.In addition, at higher redshifts we could potentially use a narrower range of wavebands (for example, g-i ) to detect lags with lengths similar to those tested here for z < 1.Finally, we could also expand on the 10 year LSST baseline by using lower cadence archival data.Future work should examine if such a cadence enables lag detection.If we use the same i -band magnitude cut-off for SDSS DR14 quasars out to z < 2 we would roughly double the number of detectable long lags.
The Wide-field Infrared Survey Explorer and eROSITA provide other catalogues of quasars not included in SDSS.In particular these catalogues could be useful for studying lower luminosity, and thus lower mass, or redder, more dust-obscured quasars.Since τ long ∝ M 1/2 , the lag timescales for these quasars will be shorter.Shorter lag timescales could help us detect a higher yield of longs lags for lower mass quasars, although this exact yield is difficult to estimate.
Even with our conservative estimates, we should be able to detect long lags for dozens to hundreds of new quasars over the ten year lifetime of LSST.In addition, Yao et al. (2022) were able to potentially detect a lag of ∼ 70 days using light curves with a less than 1 year baseline, so we speculate that we should be able to detect long lags in some fraction of nearby quasars with durations ≲ 50 days within the first 1-3 years of LSST.Roughly a third of SDSS quasars with i < 21 and z < 1 are within z < 0.5, where we expect lags to have durations of ≲ 50 days.Therefore it could be possible to detect long lags for roughly 20-100 quasars within the first couple of years of LSST.

DISCUSSION AND SUMMARY
In this Section we discuss caveats related to our simulated light curves (Section 7.1) and lag detection methods (Section 7.2) that we plan to address in future work.We then provide a summary of the paper in Section 7.3.

Improved Simulated Light Curves
In this paper we create mock observations of our driving and reprocessed light curves with both the proposed LSST baseline v2.0 and longest season cadence from OpSim.While OpSim models SNR and observing gaps by taking into account sky brightness, bad weather predicted from localized weather models, seeing, telescope maintenance times, slew times, and filter changes, the cadences and SNR measurements may still be optimistic.However, small changes in cadence should not have a large impact on our ability to detect lags on the order of months.In addition, our tests in Section 5.5 using SNR corresponding to i = 21 mag light curves show that even if the SNR corresponding to the magnitudes in Table 1 is too high by an order of magnitude, we should still be able to recover lags with similar accuracy.
We use two types of light curves, self-generated DRW light curves and light curves from the simulation in Jiang & Blaes (2020).While often used to model quasar light curves, the DRW models are not physically motivated and may fail to model the quasar PSD properly at short and long timescales.In particular, PSD from highly sampled light curves from the Kepler mission appear to have a steeper slope at short timescales (Mushotzky et al. 2011b;Kasliwal et al. 2015;Smith et al. 2018) and the length of the light curve may have an effect on the PSD at long timescales (Stone et al. 2022).On the other hand, the simulated light curves from Jiang & Blaes (2020) are not intended to capture quasar variability for a wide range of quasar parameters.Therefore, future work to simulate light curves for a wider range of quasar parameters is needed to better understand our ability to detect long lags in their variability.
We reprocess our light curves with a Gaussian response function.However, the shape of the response function for a long lag is unknown.For this paper we set the width of the Gaussian to S ≡ C/5 (see Equation ( 10)), but preliminary tests suggest that javelin has more difficulty accurately recovering lags for response functions with larger widths (S ≳ C/2).In addition, using other functional forms for the response function, such as the log normal distribution, may make it more difficult to detect a long lag using javelin.Therefore the width and shape of the response function are important parameters.Here, we use the default javelin response function model, a top hat with no priors on the width.It is possible that adjusting the model or priors for the width of the response function could improve the accuracy of javelin.Still, javelin may not be the best tool if detecting long lags relies on finding lags in variability that is reduced or smeared as it moves inwards.The maximum-likelihood method is better designed to capture the shape of the response function and should therefore be able to handle broader or more complex response functions.
Our response function is centered on frequencies corresponding to the timescale of the long lag.However, if the perturbations that are accreted inward leading to the long lag are on the thermal time scale it may be that, unlike the short lag timescale, the frequencies most affected by the long lag do not correspond exactly to the inverse of the long lag timescale.The thermal timescale is similar to the viscous timescale we use to determine the long lag timescale, for the large aspect ratios in our model.However, a future improvement would be to experiment with the recoverability of long lags depending on the frequencies most affected by the long lag.
Another simplification is that we model the radial dependence and value of the aspect ratio as independent of SMBH mass or Eddington ratio, when that is very likely not the case.Most likely, as the Eddington ratio of a quasar increases, its aspect ratio will increase, which will shorten its long lag timescale.
Ultimately, simulations of the reprocessing of high frequency light in the UV/optical emitting regions of turbulent quasar disks for a wide range of quasar parameters are needed to better model a more accurate response function for the long lag.To perform these simulations it will be necessary to use a multi-band radiation MHD code, such as the code developed by Jiang (2022) for athena ++.Radiation MHD is necessary to capture the Maxwell stress, magnetic pressure, radiative viscosity, and radiation pressure, all of which may play a role in locally generated turbulence that can be accreted inwards on the viscous timescale (Jiang & Blaes 2020).The multi-band capabilities will be useful because it will allow us to inject high-frequency radiation and compare it to the radiation reprocessed and emitted by the local disk.

Improved Lag Detection Methods
Here we use two main lag detection methods.The first, javelin, very accurately detects long lags when the ratio between the amplitude of the long lag and the amplitude of the short lag is ≳ 0.1.On the other hand, it performs very poorly when this amplitude ratio is 0.01.javelin has been shown previously to be biased in the presence of multiple lag signals (e.g., Cackett et al. 2018a).
We can potentially improve the performance of javelin by simultaneously fitting several reprocessed light curves, instead of just fitting one driving and one reprocessed light curve.Fitting several reprocessed light curves will be possible for LSST, which will have six different wavebands.With multiple bands, javelin will be less likely to misidentify seasonal gaps as lags when the lags of the different wavebands are significantly different from each other (Zu et al. 2016), which would be especially helpful if the baseline v2.0 cadence is selected.However, in our model for the long lag, the lag timescales between longer wavelength bands becomes very short.As a result, fitting multiple wavebands simultaneously may not be as helpful as it is when fitting BLR lags, which differ significantly from band to band.Another potential improvement would be to use a custom model for the response function in javelin that includes a positive and negative lag, instead of using the default model which is just a single top hat function.
The maximum-likelihood method is less accurate overall than javelin at detecting the long lag.The long lags also become increasingly skewed towards zero as the amplitude of the short lag becomes stronger relative to the long lag. Figure 7 shows that this skewing towards zero is due to the shape of the combined long and short lag response function.However, because of this bias, it will be difficult to disentangle the exact value of the long lag versus the ratio of the long lag amplitude to the short lag amplitude.Utilizing the full frequency dependent information instead of summarizing the long lag with a single average value may help disentangle the two and reveal more details about the shape of the response function.More details could also be obtained by increasing the number of bins.In this paper, we only use four frequency bins for the Fourier method, because of the relatively low sampling rate of LSST, and because it is easier to detect a single lag with fewer bins.
There are several ways that using all of the frequency dependent lag information can provide additional information on the shape of the response function.For one, regardless of the ratios between the long and short lags, the wrapping frequency seems to consistently occur very near 1/2τ long in Figure 7. Therefore, we may be able to more effectively exploit the location of the wrapping frequency to derive the lag.In addition, as can be seen in Figure 7, for larger amplitude ratios even if the value of the long lag at lower frequencies is skewed towards zero, right before the lag wraps it drops to roughly the value of the input long lag.If, on the other hand, the lags increase towards zero or positive values as the lag approaches the wrapping frequency, as is the case for smaller amplitude ratios, that would provide information on the possible presence of a long lag with a small amplitude relative to the short lag.
Overall, there is room for improvement in the lag recovery techniques used here, especially if we want to detect long lags with amplitudes less than 10% of the amplitude of the short lag.We have mentioned several ways our employment of javelin and the maximumlikelihood method could be improved.It is possible that alternative methods used for detecting lags in quasar light curves may be more accurate.The Von-Neumann method, which we also test on our mock light curves, seems like a particularly promising alternative, especially if the baseline cadence is chosen for LSST (see Section 5.3).Ultimately, the greatest improvement may be to create a new method, directly tailored to detecting the long lag.

Summary
Inspired by the recent potential detection of a long lag in Fairall 9, we have laid out a model for the timescale of negative long lags, which occur due to turbulence in the outer disk and propagate inward on the viscous timescale.Because the viscous timescale depends on the aspect ratio of the disk, measuring long lags can provide information on the vertical structure of the disk.Based on the measured long lags in Yao et al. (2022) for Fairall 9, we find that the scale height depends on the radial disk position as h ∝ r 3.59 .However, additional long lags need to be detected to properly constrain our initial model for the radial dependence of the scale height.We therefore hope to use LSST to detect these additional long lags.
To study whether LSST will be able to detect long lags, we simulate mock LSST light curves using both DRW modeled light curves and light curves from the simulation in Jiang & Blaes (2020), which we then mock observe with several LSST cadences from OpSim.We reprocess the lagging light curve with a Gaussian response function.We then test several lag detection methods, including the ICCF (Gaskell & Peterson 1987;White & Peterson 1994;Peterson et al. 2004), a Von-Neumann estimator presented by Chelouche et al. (2017), javelin (Zu et al. 2011), and a maximum-likelihood Fourier method (Zoghbi et al. 2013), on our mock light curves.We primarily focus on two main lag detection methods.The first is the commonly used javelin code.We find that when only the long lag is included (Section 5.1), javelin is accurate to within 10% of the input lag for all light curves with input long lags of −50 and −130 days, and accurate to within 10% of the input lag for 79% of light curves with input lags of −400 days in our DRW modeled light curves.
When we include a short lag as well (Section 5.2) javelin remains very accurate, unless the ratio of the amplitude of the long lag to the amplitude of the short lag is around 0.01, at which point the accuracy decreases dramatically.Conversely, the short lag is not detected by javelin unless this ratio ≤ 0.1 (for an input long lag of −50 days) or ≤ 0.01 (for input long lags of −130 and −400 days).javelin detects long lags with a similar level of accuracy even when using a DRW model to interpolate the non-DRW modeled light curves from the simulation in Jiang & Blaes (2020) in Section 5.4.
Because javelin is not aimed at detecting two lags simultaneously we also use a maximum-likelihood Fourier method, which separates lags in frequency space.Separating lags in frequency space should make it easier to detect multiple lags.However, as the amplitude of the long lag decreases relative to the amplitude of the short lag, the influence of the short lag on lower frequencies increases, skewing the recovered long lag towards shorter durations.When the amplitude ratio between the long and short lags is 0.01 the maximum-likelihood method is unable to identify the long lag at all for most light curves, just like javelin.On the other hand, the maximum-likelihood method provides a significant advantage when it comes to detecting the short lag, which it is able to accurately detect for all amplitude ratios.Finally, the performance of this method on the light curves from Jiang & Blaes (2020) is comparable to its performance on DRW light curves.
Based on the results in Section 5.3 and other cadence tests we performed, we believe a long season cadence 4 is helpful for detecting long lags.A long season cadence is particularly beneficial when using javelin to detect long lags, whereas the Von-Neumann method and the maximum-likelihood method for shorter duration long lags (∼ −50 days) are able to detect long lags with similar accuracy regardless of the cadence.However, adding just 18 LCO observations annually to mock light curves sub-sampled with the baseline v2.0 cadence, leads to a comparable or higher lag detection rate than for the long season cadence, for both our main methods for all lag durations.Additional LCO observations or the Von-Neumann method may also be effective tools for finding long lags in the Wide Fast Deep Survey, which will cover a much larger area of the sky at lower cadence.
Even with our conservative estimates in Section 6, we should be able to detect long lags for dozens to hundreds of new quasars over ten years.These detections will significantly improve our ability to model the long lag as a function of wavelength, SMBH mass, and Eddington ratio.These models will help us constrain the vertical structure of quasar disks for a range of quasar parameters.

APPENDIX
A. FULL TESTS In Sections 5.2 and 5.3 we show the distributions of lags detected by our various methods only for light curves with input lags of −50 and 7 days, for the sake of brevity.Here we present the same results, but for all lag bins in Table 1.

A.1. Long and Short Lags
We show the distributions of the long lags detected by javelin for 100 mock long season cadence LSST light curves with input long lags of −50 days (top panel), −130 days (middle panel), and −400 days (bottom panel) in Figure 23.Distributions shown in orange have only long lags, while distributions shown in blue have long and short lags with varying ratios between the amplitude of the long and short lag.The accuracy of the distributions for different input long lags is mostly consistent, with the exception of the distribution for light curves with input long lags of −130 days and an amplitude ratio of 0.01, which appears to be more accurate than the other distributions (see Section 5.2 for more discussion).
Figure 24 is the same as Figure 23 but for the maximum-likelihood method.The results for light curves with an input long lag of −130 days are mostly consistent with the results for an input long lag of −50 days.The maximumlikelihood does do slightly better at detecting some negative lags for an amplitude ratio of 0.01 for −130 day lags.This improvement is expected based on our tests of the Fourier method with even-cadence light curves shown in Figure 7.It appears that because the long lag is greater for light curves with −130 day input lags, it is not skewed all the way positive by the short lag.
We show the short lags detected by javelin and the maximum-likelihood method in Figures 25 and 26, respectively, for different input lags and amplitude ratios.As discussed in Section 5.2, javelin does a better job detecting the short lag for light curves with an input long lag of −50 days and an amplitude ratio of 0.1 than for light curves with other input lags.Otherwise, the accuracy of these distributions for different input lags is similar.

A.2. Improving the Baseline Cadence
We compare the accuracy of javelin and the maximum-likelihood method for the long season LSST cadence and baseline v2.0 cadence in Figures 27 and 28, respectively.javelin identifies the seasonal gaps of −180 days for light curves with the baseline cadence for all input lags.As mentioned in Section 5.3, the maximum-likelihood method does better at still detecting the input lags for light curves with the baseline cadence than javelin, although it does better on light curves with −50 day lags than −130 day lags, perhaps because −130 days is closer to the length of the season gaps.
In Figures 29 and 30, we compare light curves that have both the LSST baseline cadence and a roughly ten day cadence of LCO observations in the LSST season gaps (see Section 5.3 for details), to light curves with the LSST baseline cadence for various wavebands.We find that adding LCO observations in the LSST seasonal gaps is very useful for improving the accuracy of javelin for all input long lags.The distributions of lags detected by the maximum-likelihood method for light curves with different cadences are consistent with each other for all input long lags.However, Figure 31 shows that the maximum-likelihood method does a better job detecting short lags with the added LCO cadence for both input lags.
As mentioned in Section 5.3, the Von-Neumann method performs better on light curves with the baseline cadence than javelin.We show the distribution of lags detected by the Von-Neumann method for light curves with the long season (in orange) and baseline (in green) cadences in Figure 32 the steps outlined in Section 3, but skip Section 3.2, i.e. we do not reprocess either light curve with a response function.
As a result, we have two identical light curves that are then sub-sampled with the cadence for two different wave bands (u and y or g and z, depending on the cadence).Below we discuss the signals recovered by javelin (Section B.1) and the maximum-likelihood method (Section B.2). Any signals recovered could be due to the cadence/gaps in the cadence, or timescales inherent to the light curves themselves, such as the damping timescale.

B.1. Javelin
Figure 33 shows the distributions of the medians of the negative (left panel) and positive (right panel) values of the javelin distributions for 100 mock DRW light curves with no re-processing.Light curves in the top panel have been mock observed with the long season LSST u-and y-band cadence used for the light curves in Sections 5.1, 5.2, 5.4, and 5.5.In the remaining panels light curves have been mock observed with the baseline u-and y-band cadence, the baseline g-and z -band cadence, and the baseline g-and z -band cadence with the additional LCO cadence, which are all used to sub-sample the light curves in Section 5.3.The distribution of the negative medians for light curves with the long season cadence peaks at zero.The duration of the long lag detected by javelin for 70% of light curves is < 15 days, suggesting that the lags recovered by javelin in the previous sections are due to the input lag signal.On the other hand, over 50% of positive lags recovered by javelin are over 5 days, suggesting that our short lag results could be contaminated.The individual javelin distributions for light curves with the long season cadence that have not been reprocessed with an input lag, several of which are shown in Figure 34, peak most prominently at zero, but have noise.Most of

Figure 2 .
Figure2.The normalized distribution of τ long (z) for quasars in SDSS DR7(Shen et al. 2011) that fall within different redshift ranges.We assume α = 0.1, and take the masses, Eddington ratios, and redshifts fromShen et al. (2011).In the top left panel we use β = 3.6 as inYao et al. (2022); in the top right panel we use β = 2, which is a more conservative estimate that the aspect ratio increases linearly with radius; in the bottom left panel we use β = 0 and h/r = 0.01, which are typical values in the standard thin disk model and lead to much long lag timescales (see different x-axis range)

Figure 3 .
Figure 3. Top panel: Simulated driving u-band light curve, in blue, and y-band light curve, reprocessed with a long lag of −50 days, in orange.Middle panel: The same light curve as in the top panel sub-sampled using the long season LSST cadence.Bottom left: javelin probability distribution for the long lag between the light curves above.The long lag is shown as a vertical dashed line.The median of the negative values of the distribution is −50 days.Bottom right: An example of the lags as a function of frequency recovered by the maximum-likelihood method for a mock light curve with an input long lag of −50 days (shown as the horizontal dashed blue line).The points are the lag detected in each bin (shown as the dashed vertical red lines).The vertical lines show the error bars.The lowest frequency bin has the largest error bar.The next two bins are roughly consistent with the input lag and have an error-weighted mean of −48 days.Because there is no short lag and it is at frequencies higher than the wrapping frequency (vertical grey line), the highest frequency bin is consistent with zero.

Figure 4 .Figure 5 .
Figure 4.The cross-correlation coefficients of even daily cadence light curves that have been reprocessed with long (short) lags of −50 (7), −130 (9), and −400 (12) days, in the top, center, and bottom panel, respectively.Distributions shown in orange have been reprocessed with only a long lag.The ratios of the amplitude of the long to the amplitude of the short lag for the light curves with distributions shown in green, light blue, and pink, are 0.2, 0.1, and 0.01, respectively.The dashed vertical lines show the input lags.

Figure 6 .
Figure 6.An example javelin distribution for a mock light curve with the baseline cadence.While there is a very small peak around the input lag (shown as a dashed line) of −50 days, the peaks around ±180 days are larger.

Figure 7 .
Figure7.The time lags as a function of frequency between two evenly sampled, daily cadence light curves found using the phase of their complex cross-spectrum and Equation (15).The input long (short) lags are −50 (7) and −130 (9) days in the top and bottom panel, respectively.Distributions shown in orange have been reprocessed with only a long lag.The ratio between the amplitudes of the long and short lag for the light curves with distributions shown in green, light blue, and pink, are 0.2, 0.1, and 0.01, respectively.The dashed horizontal lines show the input lags and zero and the vertical grey lines show the wrapping frequencies.

Figure 8 .
Figure 8.The distribution of the inverse of the time between each observation in the LSST long season cadence.

Figure 9 .Figure 10 .
Figure 9.The distribution of the long lags detected by javelin for 100 mock long season cadence LSST light curves with input long lags of −50, −130, and −400 days (shown as dashed black lines) in the top, center, and bottom panel, respectively.

Figure 11 .
Figure 11.The distributions of the long lags detected by javelin (top panel) and the maximum-likelihood method (bottom panel) for 100 mock long season cadence LSST light curves with an input long lag of −50 days (shown as a dashed black line).Distributions shown in blue are for light curves that also have a short lag of 7 days.The legends show the ratio of the amplitude of the long lag to the amplitude of the short lag for each blue distribution.

Figure 12 .
Figure 12.The distribution of the short lags detected by javelin (top panel) and the maximum-likelihood method (bottom panel) for 100 mock long season cadence LSST light curves with different amplitude ratios (shown in the legend).In the top panel we zoom in on the portion of the distribution between 0 and 20 days to better show the accuracy of the short lag when detected.The input short lag, shown as the vertical dashed line, is 7 days.

Figure 13 .
Figure13.The fractional error of the median of the distribution of long lags recovered by javelin (points) and the Fourier method (diamonds) for 100 mock long season cadence LSST light curves as a function of the ratio of the amplitude of the long lag to the amplitude of the short lag for input long lags of −50, −130, and −400 days in blue, purple, and orange, respectively.The error bars represent the fractional error of one standard deviation for each distribution and the horizontal dashed line shows zero.A slight scatter is added in the x-direction to make the individual points easier to view.

Figure 14 .
Figure 14.The distributions of the long lags recovered by javelin (top panel) and the maximum-likelihood method (bottom panel) for 100 light curves with input lags of −50 and 7 days.The green and orange distributions are for mock light curves with the long season and baseline v2.0 LSST cadences, respectively, for the u-and y-band.The dashed vertical line shows the input long lag.

Figure 15 .
Figure15.Distribution of the long lags recovered by javelin (top panel) and the maximum-likelihood method (bottom panel) for 100 mock light curves with input lags of −50 and 7 days (top panel).The dashed vertical line shows the input long lag.The green and light blue distributions are for mock light curves with the baseline v2.0 LSST cadence for the u-and y-band and g-and z -band, respectively.The distributions shown in pink are for light curves with the baseline g-and z -band cadence plus an LCO observation taken roughly every ten days in each LSST observing gap.These LCO observations amount to around 18 observations in each gap, or per year, for a total of around 180 observations.

Figure 16 .Figure 17 .
Figure16.Distribution of the short lag recovered by the maximum-likelihood method for 100 mock light curves with input lags of −50 and 7 days.The dashed vertical line shows the input short lag.The green and light blue distributions are for mock light curves with the baseline v2.0 LSST cadence for the u-and y-band and g-and z -band, respectively.The distributions shown in pink are for light curves with the baseline g-and z -band cadence plus an LCO observation taken roughly every ten days in each LSST observing gap.These LCO observations amount to around 18 observations in each gap, or per year, for a total of around 180 observations.

Figure 18 .
Figure 18.Distribution of the long lags recovered by javelin for 16 mock long season cadence LSST light curves made using light curves from the simulation in Jiang & Blaes (2020) for input lags of −50 and 7 days (top panel), −130 and 9 days (middle panel), and −400 and 12 days (bottom panel).The long lags are shown as the dashed vertical line.

Figure 19 .
Figure 19.Distribution of the long lags detected by the maximum-likelihood method for 16 mock long season cadence LSST light curves made using light curves from the simulation in Jiang & Blaes (2020) for input lags of −50 and 7 days (top panel) and −130 and 9 days (bottom panel).The long lags are shown as the dashed vertical line.

Figure 20 .
Figure 20.Distribution of the short lags detected by the maximum-likelihood method for 16 mock long season cadence LSST light curves made using light curves from the simulation in Jiang & Blaes (2020) for input lags of −50 and 7 days (top panel) and −130 and 9 days (bottom panel).The short lags are shown as the dashed vertical line.

Figure 21 .
Figure 21.The distributions of the long lags recovered by javelin for 100 mock light curves with the long season LSST cadence and input lags of −50 and 7 days (top panel), −130 and 9 days (middle panel), and −400 and 12 days (bottom panel).The vertical dashed lines show the input long lag.The colored lines are the distributions for light curves where we set i = 21 to explore how a lower SNR changes our ability to detect lags.The distributions of these lags are very similar to the distributions shown in gray which are for light curves with the i-band magnitudes in Table1.

Figure 22 .
Figure22.The distributions of the long lags recovered by the maximum-likelihood method for 100 mock light curves with the long season LSST cadence and input lags of −50 and 7 days (top panel) and −130 and 9 days (bottom panel).The vertical dashed lines show the input long lag.The colored lines are the distributions for light curves where we set i = 21 to explore how a lower SNR changes our ability to detect lags.The distributions of these lags are very similar to the distributions shown in gray which are for light curves with the i-band magnitudes in Table1.

Figure 23 .
Figure 23.The distribution of the long lags detected by javelin for 100 mock long season cadence LSST light curves with input long lags, shown by the dashed vertical lines, of −50 days (top panel), −130 days (middle panel), and −400 days (bottom panel).Distributions shown in blue also include a short lag of 7 days (top panel), 9 days (middle panel), and 14 days (bottom panel).The ratio of the amplitude of the long lag to the amplitude of the short lag is shown in the legends for each blue distribution.

Figure 24 .Figure 25 .
Figure 24.Distribution of the long lags recovered by the maximum-likelihood method for 100 mock long season cadence LSST light curves with input long lags, shown by the dashed vertical lines, of −50 days (top panel) and −130 days (bottom panel).Distributions shown in blue also include a short lag of 7 days (top panel) and 9 days (bottom panel).The ratio of the amplitude of the long lag to the amplitude of the short lag is shown in the legends for each blue distribution.

Figure 26 .Figure 27 .Figure 28 .
Figure 26.Distribution of the short lags recovered by the maximum-likelihood method for 100 mock long season cadence LSST light curves with different amplitude ratios (shown in the legend).The input short lags, shown as the vertical dashed line, are 7 days in the left panel and 9 days in the right panel.

Figure 29 .Figure 30 .
Figure29.Distribution of the long lags recovered by javelin for 100 mock light curves with input lags of −50 and 7 days (left panel), −130 and 9 days (middle panel), and −400 and 12 days (right panel).The dashed vertical line shows the input long lag.The green and light blue distributions are for mock light curves with the baseline v2.0 LSST cadence for the u-and y-band and g-and z -band, respectively.In addition to the baseline LSST cadence, the distributions shown in pink have an LCO observation taken roughly every ten days in each LSST observing gap.These LCO observations amount to around 18 observations in each gap, or per year, for a total of around 180 observations.

Figure 31 .
Figure31.Distribution of the short lag recovered by the maximum-likelihood method for 100 mock light curves with input lags of −50 and 7 days (left panel) and −130 and 9 days (right panel).The dashed vertical line shows the input short lag.The green and light blue distributions are for mock light curves with the baseline v2.0 LSST cadence for the u-and y-band and gand z -band, respectively.In addition to an LSST cadence, the distributions shown in pink have an additional LCO observation taken roughly every ten days in each LSST observing gap.

Figure 33 .
Figure33.The distribution of the negative (left panels) and positive (right panels) medians of javelin distributions for 100 mock light curves with no input lag.From top to bottom, the light curves have been sub-sampled with the long season LSST uand y-band cadence, the baseline u-and y-band cadence, the baseline g-and z -band cadence, and the baseline g-and z -band cadence + LCO cadence.

Figure 34 .Figure 35 .
Figure 34.Eight individual javelin distributions for light curves with no reprocessing and the long season LSST cadence.Each panel is labeled with the median of the positive and negative values of the distribution.

Figure 36 .
Figure 36.The top panel shows the lags recovered by javelin and the bottom panel shows the lags recovered by the maximumlikelihood method for 16 mock light curves made from the simulation in Jiang & Blaes (2020) with no input lag and the long season LSST u-and y-band cadence.