Examining AGN UV/optical Variability Beyond the Simple Damped Random Walk

We present damped harmonic oscillator (DHO) light-curve modeling for a sample of 12,714 spectroscopically confirmed quasars in the Sloan Digital Sky Survey Stripe 82 region. DHO is a second-order continuous-time autoregressive moving-average (CARMA) process, which can be fully described using four independent parameters: a natural oscillation frequency ($\omega_{0}$), a damping ratio ($\xi$), a characteristic perturbation timescale ($\tau_{\mathrm{perturb}}$), and an amplitude for the perturbing white noise ($\sigma_{\mathrm{\epsilon}}$). The asymptotic variability amplitude of a DHO process is quantified by $\sigma_{\mathrm{DHO}}$ -- a function of $\omega_{0}$, $\xi$, $\tau_{\mathrm{perturb}}$, and $\sigma_{\mathrm{\epsilon}}$. We find that both $\tau_{\mathrm{perturb}}$ and $\sigma_{\mathrm{\epsilon}}$ follow different dependencies with rest-frame wavelength ($\lambda_{\mathrm{RF}}$) on either side of 2500 \AA, whereas $\sigma_{\mathrm{DHO}}$ follows a single power-law relation with $\lambda_{\mathrm{RF}}$. After correcting for wavelength dependence, $\sigma_{\mathrm{DHO}}$ exhibits anti-correlations with both the Eddington ratio and the black hole mass, while $\tau_{\mathrm{perturb}}$ -- with a typical value of days in the rest-frame -- shows an anti-correlation with the bolometric luminosity. Modeling AGN variability as a DHO offers more insight into the workings of accretion disks close to the supermassive black holes (SMBHs) at the center of AGN. The newly discovered short-term variability (characterized by $\tau_{\mathrm{perturb}}$ and $\sigma_{\mathrm{\epsilon}}$) and its correlation with bolometric luminosity pave the way for new algorithms that will derive fundamental properties (e.g., Eddington ratio) of AGN using photometric data alone.

The UV/optical luminosity of active galactic nuclei (AGN) 1 , or quasars at the bright end, is known to vary at the 10% flux level on average from weeks to years (Vanden Berk et al. 2004;Sesar et al. 2007). The time variability of AGN luminosity has been known for decades (Matthews & Sandage 1963), however, the physical mechanisms driving such variability are still unclear. Nevertheless, the success of reverberation map-ping (Peterson et al. 2004) has shown that the broad emission lines respond to and lag behind the continuum fluctuations, suggesting an accretion disk origin of the UV/optical continuum variability.
Under the assumption that the optical variability originates from the accretion disk close to the SMBH, various models have been proposed to explain the observed variability. Based on the standard α disk model (Shakura & Sunyaev 1973), some have shown that the observed optical variability could be driven by variations in the mass accretion rate (Pereyra et al. 2006;Li & Cao 2008;Liu et al. 2008). Meanwhile, others suggested that it is also possible for the accretion disk to passively reprocess the radiation from the X-ray corona given the observed short time lags between UV/optical continuum light curves (Wanders et al. 1997;Sergeev et al. 2005).
Observationally, significant efforts/progress have also been made to investigate the physical origin(s) of the UV/optical variability by virtue of exploring the correlations of variability signatures with the physical properties of AGN. Notable results include: an anti-correlation of variability amplitude with rest-frame wavelength and an anti-correlation of variability with luminosity and/or Eddington ratio (e.g., Vanden Berk et al. 2004;Wilhite et al. 2007;Bauer et al. 2009;MacLeod et al. 2010;Simm et al. 2016;Caplar et al. 2017). Correlations between variability and black hole mass have also been reported in numerous studies (Wold et al. 2007;Wilhite et al. 2007;Bauer et al. 2009;Kelly et al. 2009;MacLeod et al. 2010;Caplar et al. 2017).
Among the different techniques utilized to characterize AGN variability, Kelly et al. (2009) started a new era of directly modeling (inherently non-periodic) AGN light curves using stochastic diffusion processes, in particular, as a damped random walk (DRW) model. The DRW model features a fixed-slope power spectrum density (PSD) at high frequencies (short timescales) and a flat PSD at timescales longer than a characteristic timescale (τ DRW ). Modeling AGN variability as a DRW has provided a lot of insight into how AGN luminosity varies in UV/optical and what might be driving it (MacLeod et al. 2010). However, better sampled light curves from the Kepler Mission (Borucki et al. 2010) cast doubt on the DRW description of AGN variability because of the steeper slopes observed in the PSDs at high frequencies (Mushotzky et al. 2011); investigations carried out by other groups also arrived at similar conclusions (Kasliwal et al. 2015;Simm et al. 2016;Smith et al. 2018). This discrepancy motivates the search for new models (and methods) to analyze AGN light curves. Given that the DRW model is the simplest case of a more general class of stochastic diffusion processes, namely the continuous-time autoregressive moving-average processes (CARMA ;Brockwell 2001;Roux 2002), Kelly et al. (2014) set up a more flexible framework to model astronomical time series as CARMA processes, where the PSDs of higher-order CARMA processes can take more flexible shapes, for example, a wide range of PSD slopes can be achieved at high frequencies. Later, Kasliwal et al. (2017) demonstrated that the CARMA(2,1) model is a better fit than all other models of CARMA for a well-monitored object Zw229-15. Inspired by the aforementioned discrepancy and the pilot investigations carried out by Kelly et al. (2014) and Kasliwal et al. (2017), Moreno et al. (2019, hereafter M19) conducted an in-depth exploration of the CARMA(2,1) model, oth-erwise known as the (perturbation-driven) damped harmonic oscillator (DHO) model, and established guidelines for modeling AGN light curves as DHOs.
Here, we build upon the work performed by Kelly et al. (2014), Kasliwal et al. (2017), and Moreno et al. (2019) to model a large statistical sample of AGN as DHOs, examine the variability signatures extracted by the DHO model, and explore the potential correlations between DHO parameters and the physical properties of AGN. We acknowledge that the DHO model (alike the DRW model) is a statistical model rather than a physical model; however, stochastic diffusion processes such as CARMA are a natural choice for parameterizing AGN light curves, and it can reveal interesting variability features embedded in the light curves that could not be unfolded otherwise (see Vio & Andreani (2018) for a discussion about the limitations of CARMA modeling).
In Section 2, we introduce the data set utilized in this work and outline our initial light curve construction procedures. In Section 3, we provide an overview of the DHO model, its key features, and how to use Gaussian process (GP) to extract DHO parameters from light curves. In Section 4, we present the results of fitting DHO to our quasar light curves, layout and test our bad-fit identification algorithm, and explore the effects of light-curve sampling and photometric accuracy on the best-fit DHO parameters. In Section 5, we describe the observed correlations of DHO features with the physical properties of the quasars in our sample, and discuss the associated implications in Section 6. Finally, we summarize our results and provide an outlook for future work in direct modeling of AGN light curves using stochastic diffusion processes in Section 7.

THE DATA SET
We compiled a sample of 12,714 spectroscopically confirmed quasars discovered in the Sloan Digital Sky Survey (SDSS) Stripe 82 region (York et al. 2000;Annis et al. 2014), a 120 • long and 2.5 • wide stripe centered along the celestial equator, from the quasar catalog of SDSS Data Release 16 (DR16Q; Lyke et al. 2020). We refer to this initial sample of quasars as the main sample; quasars in this main sample either have their fundamental physical properties (e.g., black hole mass) estimated by Shen et al. (2011) or C IV emission line properties (i.e., equivalent width and blueshift) determined by Rankine et al. (2020). Figure 1 shows the distribution of the main sample quasars in the luminosity and redshift space, where the luminosity (i-band absolute magnitude) has been k-corrected to z = 2 (Richards et al. 2006). Distribution of the main sample quasars in the luminosity (absolute i-band magnitude k-corrected to a redshift of 2) and redshift space. This sample is constructed from the SDSS DR16 quasar catalog-the latest release of SDSS quasar spectra (Lyke et al. 2020); quasars in this sample either have their physical properties estimated by Shen et al. (2011) or C IV emission line properties measured by Rankine et al. (2020).

The Sloan Digital Sky Survey (2000-2008)
The Sloan Digital Sky Survey obtained images for more than 10, 000 deg 2 of the northern hemisphere down to limiting magnitudes of 22.5, 23.2, 22.6, 21.9, 20.8 at the 50% completeness level in the u, g, r, i, z bands, respectively. The SDSS Stripe 82 (S82) region was observed repeatedly over a 8-year long baseline providing up to 90 single-epoch observations (Frieman et al. 2008;Sako et al. 2008;Abazajian et al. 2009). The SDSS light curves used in this investigation span two phases of the Sloan Digital Sky Survey, namely, the SDSS Legacy Survey (York et al. 2000) and the SDSS-II Supernova Survey (Frieman et al. 2008;Sako et al. 2008), where the latter was performed under less photometric conditions. A so-called "ubercalibration" (Ivezić et al. 2007;Padmanabhan et al. 2008;Bramich et al. 2008), which takes advantage of the overlap between adjacent imaging runs to arrive at a uniform internal calibration, was utilized to achieve a ∼1% photometric accuracy in griz and 2% in u for both photometric and less photometric observations. Such a calibration method has become the default since the seventh data release of SDSS (DR7) (Abazajian et al. 2009), from which our light curves were constructed.

SDSS Light Curves
The SDSS light curves for our quasars were generated by cross matching against the photoobj table under the Stripe82 context on CasJobs 2 using a 1" matching radius. We imposed an initial quality cut on the matched detections-that is, their photometry must be "clean" and their photometric uncertainties (psfmagerr) must be smaller than 1 mag. More on the "clean" flag (and other photometry flags) can be found on SDSS-IV's website 3 . The raw light curves are then post-processed following the recipe described below: 1. We first require light curves to have at least 30 epochs.
2. We then remove data points that deviate more than 3σ away from the 3-point running median. This process removes any abnormally largeamplitude variability in the photometry (Graham et al. 2014).
3. Lastly, we remove data points with photometric uncertainties that are more than 5 times larger than the median uncertainty of all photometry in the corresponding light curves.
Note that the post-processing described above is on a per-band basis, that is, failing to have a good light curve in a particular passband does not exclude an object from our sample. The final collection contains ≈12,400 light curves in each u, g, r, i, z band. Those light curves are then fitted with the DHO model. The distributions of four basic statistics of those light curves are shown in Figure 2.
3. METHODOLOGY We model the time variability of AGN UV/optical luminosity as a DHO, which is formally defined as the solution to the following stochastic differential equation, where (t) is Gaussian white noise with an amplitude of unity 4 , α 1 and α 2 are called the autoregressive (AR) coefficients 5 , and β 0 and β 1 are called the moving-average (MA) coefficients. The differentiation is with respect to time.
Here, x represents the brightness (magnitude 6 in this work) of the modeled quasars. For a comparison, the stochastic differential equation for the DRW model-a CARMA(1,0) process-has the form 7 ,

Damped Harmonic Oscillator (DHO)
The utility of modeling quasar light curves as DHOs has been explored and discussed extensively in M19; here we provide a brief introduction to the DHO model under the framework of an impulse-response dynamical system.
In short, we can think of time series (or light curves in this context) following the stochastic differential equation (SDE) shown in Equation 1 as impulse-response dynamical systems. The left-hand-side (LHS) of the SDE describes how the systems respond to impulse perturbations (in the differential form) and the right-handside (RHS) describes how the systems are being perturbed/excited. This interpretation connects the DHO process naturally to the classical damped harmonic oscillator, where the classical counterpart has a deterministic driving force on the RHS rather than a stochastic one. Given this analogy, we can rewrite Equation 1 as, where 2ξω 0 = α 1 and ω 2 0 = α 2 . ξ is the damping ratio of the damped oscillator and ω 0 is the natural oscilla-tion frequency (i.e., when there is no damping). We can classify DHO processes into underdamped (ξ < 1) and overdamped (ξ > 1) DHOs, each corresponds to a different class of dynamical systems that can be interpreted using the impulse-response framework (Moreno et al. 2019). On the RHS of Equation 3 (compared to Equation 1), we renamed β 0 to σ and defined τ perturb as β 1 /β 0 . σ can be treated as the amplitude of the short-term perturbing white noise ( (t)) with the unit 8 of magnitude/time 3/2 . For unit consistency, τ perturb obtains a unit of time and manifests as a characteristic timescale of the perturbation process (Kasliwal et al. 2017;Moreno et al. 2019).
Given Equation 3, we can fully define a DHO process using ξ, ω 0 , σ , and τ perturb , where the first two are from the LHS of Equation 3 and the last two are from the RHS of Equation 3. This new set of independent parameters and the four (α 1 , α 2 , β 0 , and β 1 ) from Equation 1 will be used side by side throughout this manuscript, where the original parameters (from Equation 1) will be used mostly in the technical sections (e.g., Section 3) and the newly defined parameters will be used primarily in the discussion of scientific implications (e.g., Section 5).
Various intrinsic timescales can be extracted from the LHS of Equation 1 (or Equation 3) based on the roots (r 1 , r 2 ) of its characteristic equation, When ξ < 1 (underdamped DHOs) the two roots are complex conjugates, a decay timescale (τ decay ) and a damped oscillation period (T dQPO ) can be obtained, The unit of σ can be derived by matching the units of the LHS and the RHS of Equation 3 (Kasliwal et al. 2017) given that (t) (≡ dW/dt) has the unit of 1/ √ dt (Roux 2002).
where 2π/ω 0 is the natural oscillation period (T QPO ) associated with ω 0 . T dQPO is effectively the oscillation period in the presence of resistance/damping. For no resistance/damping (ξ ≈ 0), T dQPO and T QPO become equivalent. When ξ > 1 (overdamped DHOs), both roots are real leaving us a rising timescale (τ rise ) and a decay timescale (τ decay ), These four timescales (two for each class: underdamped and overdamped) derived from the LHS of Equation 3 set the foundation for describing the response of the modeled dynamical system to a delta function impulse perturbation, whereas τ perturb from the RHS of Equation 3 characterizes how the dynamical system is being driven/excited. In Figure  On top of the intrinsic timescales, M19 also defines a decorrelation timescale, which characterizes the timescale at which the system becomes de-correlated from an earlier excitation (forgets about its past self), Lastly, the asymptotic root-mean-square (RMS) amplitude of a DHO process (σ DHO ), which is jointly determined by the parameters in the LHS and RHS of Equation 1, can be computed using (Brockwell 2001),

DHOs as GPs
For a given time series (light curve), CARMA (DHO in this work) parameters are commonly extracted by 9 We note that the driving component PSD presented here is not necessarily representative of that for the true underlying physical process, rather the features embedded (e.g., τ perturb ) in this PSD will inform us about important characteristics of the true physical model. We refer interested readers to Jones & Ackerson (1990), Brockwell (2001), and Kasliwal et al. (2017) for further details on the analytic form of DHO's PSD. maximum likelihood, where the likelihood function can be calculated through Kalman recursion in the "statespace" of CARMA (Jones & Ackerson 1990;Brockwell 2001;Kelly et al. 2014;Kasliwal et al. 2017). Recently, Foreman-Mackey et al. (2017) introduced a new algorithm for performing fast Gaussian process (GP) modeling and suggested that the likelihood function of CARMA processes can be computed in O(N J 2 ) based on this new algorithm (N is the number of data points in a time series and J is the autoregressive order (p) of a CARMA model); this new algorithm is also up to 10 times faster than the Kalman recursion method. In this work, we adopt the algorithm introduced by Foreman- Mackey et al. (2017) and express DHO processes as a special class of GPs.
Although rarely recognized, a CARMA process driven by Gaussian noise (i.e., (t) is Gaussian) is a Gaussian process, which makes it viable to calculate the likelihood function of CARMA models using GPs. The speedup demonstrated in Foreman-Mackey et al. (2017) originates from the fact that the CARMA auto-covariance function can be formulated as a sum of complex exponentials allowing the covariance matrix to be semiseparable, thus enabling faster computation of the likelihood function (Ambikasaran 2015;Foreman-Mackey et al. 2017). Below we demonstrate how to represent a DHO's auto-covariance function (or entries in the autocovariance matrix for discretely sampled data) in terms of celerite kernels-the actual implementation of the algorithm presented in Foreman-Mackey et al. (2017). The full derivation for CARMA processes of all orders is beyond the scope of this work and can be found in (Yu et al. 2022, in preparation).
From Equation 4 in Kelly et al. (2014), we can write out the auto-covariance function of a DHO process, where r 1 , r 2 are the two roots of the characteristic polynomial associated with the LHS of Equation 1 and τ is the positive time lag between any two timestamps. Note that Kelly et al. (2014) factored out β 0 and called it σ-equivalent to σ defined in the previous section. When r 1 , r 2 are real (overdamped DHOs: ξ > 1), A 1 , A 2 are also real, and Equation 9 becomes a sum of two real exponential GP kernels (celerite real term). In the celerite framework, Equation 9 for overdamped  DHOs can be written as: where τ nm is the positive time lag between the n th and m th data point in a time series, σ n is the measurement uncertainty on the n th data point; a 1 = A 1 , a 2 = A 2 , c 1 = −r 1 , and c 2 = −r 2 . Note that for a discretely sampled DHO process, k (τ nm ) gives the entries in the corresponding auto-covariance matrix.

Fitting DHO to Quasar Light Curves
In the previous section, we demonstrated that the celerite framework can be utilized to compute the likelihood function of CARMA (used DHO as an example), however, the native celerite software does not come with the functionality (the actual code) to fit an arbitrary CARMA model that is more complex than the DRW to light curves. To facilitate general-purpose CARMA modeling taking advantage of celerite, we implemented the appropriate mapping from the CARMA parameterization to the celerite parameterization in a new Python package-EzTao, which was used to fit DHO to our quasar light curves.
The likelihood landscape of a complex GP kernel (such as that of a DHO) is usually non-convex (e.g., having multiple local optima). Thus, during the fitting process we randomly initialized 100 optimizers across the DHO parameter space that can be probed by the temporal sampling of the input light curves and selected the maximum a posterior (MAP) estimation as the best-fit solution, where very broad flat priors were used to prevent the potential numerical overflow/underflow caused by catastrophic runaways of the optimizers 10 . We performed the fitting process 5 times for each light curve to make sure that a robust fit was obtained. Our experience showed that neither increasing the number of optimizers nor repetitions will change the final distribution of the best-fit DHO parameters for our quasar sample.
After we obtained the best fit for each object, we used emcee (Foreman-Mackey et al. 2013), a python implementation of Goodman & Weare's Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler (Goodman & Weare 2010), to sample the posterior distribution with the MCMC walkers initiated at the MAP position. One additional prior was used to restrict the MCMC walkers from potential catastrophic runaways: 10 −3 days < log (τ perturb ) < 10 5 days. We ran MCMC for 15,000 steps using 32 walkers and discarded the first 3,000 steps as the "burn-in". The largest autocorrelation time for all chains is around 500 steps. The uncertainty of the best-fit DHO parameter is taken as the "1-sigma" range (one half the central 68.3% interval) of the marginalized posterior distribution.

BEST-FIT DHO PARAMETERS
10 The AR coefficients (α 1 and α 2 from the LHS of Equation 1) have a boundary of [-15, 15] in the natural log scale, and the MA coefficients (β 0 and β 1 from the RHS of Equation 1) have a boundary of [-23, 7] in the natural log scale

The Distribution of DHO Parameters and The Identification of Bad Fits
Just because a fit is robust (in that it does not change across multiple independent runs) does not mean that it is accurate and precise. We identified two main "dead zones" in the DHO parameter space that are hosts of bad DHO fits, these are MAP fits that are catastrophic failures potentially as a result of insufficient temporal sampling or large photometric uncertainty. DHO fits that end up in those regions, as listed below, were flagged as bad and removed from our sample: 1. DHO fits corresponding to timescales (see Section 3.1) that are either longer than the span of the light curve (Max ∆t) 11 or shorter than one half the minimum separation between any two observations (Min ∆t).
The first cut above alone removed ≈70% of the objects from our initial sample in every photometric band. Independently from the first cut, the second cut alone removed ≈55% of all objects in u and z, and ≈20% of all objects in gri. After these two cuts, we were left with 2997, 4570, 3530, 2999, 1180 DHO fits in ugriz bands, respectively. Next, we used an isolation forest (an outlier detection algorithm) (Liu et al. 2008) to further identify and remove bad fits. An isolation forest takes a data set and splits it randomly until no more split can be made. If anomalies/outliers are rare and different from the main population, it would take longer (more splits) to isolate a "regular" data point than an outlier. Therefore, those that are first isolated are identified as outliers by an isolation forest. Running an isolation forest on the current sample of DHO fits in each band separately removed an additional 5% of objects, our final sample contains 2847, 4341, 3353, 2849, 1121 good DHO fits in u, g, r, i, z bands, respectively.
The distribution of good DHO fits (using g-band as an example) is shown in the top two panels of Figure  In the top left panel of Figure 3 we can identify two main clusters, one for underdamped DHOs (log(ξ) < 0) and one for overdamped DHOs (log(ξ) > 0), separated by the vertical solid line at log(ξ) = 0. The clustering of underdamped/overdamped populations is not a characteristic of the DHO model rather that of the quasar light curves, that is, any point in the DHO parameter space constitutes a valid model. We can further classify underdamped DHOs into quasi-periodic (QPO) and oscillatory DHOs with a (empirically chosen) dividing T dQPO of 600 days (dashed line). Despite both being underdamped, QPOs generally have larger natural oscillation frequencies (ω 0 ) and smaller damping ratios than oscillatory DHOs. Both underdamped and overdamped DHOs are colored by their τ decorr but using different color maps. The dotted purple lines on top of the overdamped population demonstrate the evolution of τ rise across the response parameter space, where the numbers give the τ rise in the unit of days. In the top right panel (the perturbation parameter space), objects are also colored by τ decorr . Since we do not see an obvious color gradient as τ perturb increases (or decreases) in the overdamped cluster, we argue that τ decorr and τ perturb are largely uncorrelated for overdamped DHOs, which suggests that the perturbation component (RHS of Equation 3) and the response component (LHS of Equation 3) of the modeled quasars might be decoupled. Such decoupling is expected if the perturbation process is external to instead of originating in the accretion disk of AGN (see Section 6 for further discussion).
Three representative points are chosen from the distribution of good DHO fits, one for each identified class: overdamped DHO, QPO, and oscillatory DHO, to investigate/visualize the intrinsic variability signatures embedded in the DHO parameters; they are labeled using a square (overdamped DHO), triangle (QPO), and a cross (oscillatory DHO). The bottom three panels of Figure 3, from left to right, show the perturbation spectrum-the driving process PSD (RHS of Equation 3), the impulseresponse function (LHS of Equation 3), and the final (simulated) DHO light curves, of the dynamical systems described by the corresponding DHO parameters at the selected points.
It can be clearly seen from these three example DHOs that different regions of the parameter space correspond to different unique characteristics. The underdamped population-log(ξ) < 0 in the top left panel of Figure 3 and the top-left corner in the top right panel of Figure 3-has large τ perturb and the light curve is smooth at short timescales but bumpy at long timescales; comparatively, the QPO subclass shows a much stronger peri-odicity than the oscillatory subclass. On the other hand, the overdamped population-log(ξ) > 0 in the top left panel of Figure 3 and the bottom-right portion in the top right panel of Figure 3-has smaller τ perturb and the light curve is bumpy at short timescales but smooth over large timescales.
Next, to confirm that our bad-fit identification procedure works as expected, we simulated ∼12,000 DHO light curves sampled at the exact temporal cadence and photometric accuracy of the quasars in our main sample, and then fitted them with DHO. The input model parameters for the simulated light curves were drawn from the distribution of g-band good fits (see top two panels Figure 3). Figure 4 shows the distribution of the MAP best fits obtained on these simulated data in the ξ, ω 0 parameter space. Based on the quality cuts described at the beginning of this section, DHO fits that passed the cuts are color-coded by the difference in log(ω 0 ) between the MAP value and the input; the simulated overdamped DHOs that failed to pass the cuts are shown using gray dots whereas the failed underdamped DHOs are shown using magenta crosses. We can see that the distribution of the identified good fits largely overlaps that of their input outlined by the red contour. We can also see that the bad fits that were simulated as overdamped DHOs spread over the whole ξ, ω 0 parameter space and the bad fits that were simulated as underdamped DHOs are more concentrated around the distribution of their input. The orange dashed line at the top left corner corresponds to a T QPO of one day suggesting that best fits accumulated in this region might be fitting the daily observing cadence rather than the intrinsic timescales in the quasar light curves. A cut of log(ξ) -log(ω 0 ) > 1 (the cyan dash-dotted line in Figure 4) effectively removes those suspicious fits. Given the resulting distribution of good/bad fits obtained on simulated light curves, we can confirm that the good fits identified through the criteria established at the start of this section are effective. However, we do find two flavors of overestimation/underestimation from the distribution of good fits: a systematic offset (between the colored distribution and the red contour) and a parameter-dependent trend (the color gradient within the colored distribution); we will discuss both scenarios with more detail in the next section.

The Effects of Time Sampling and Photometric Accuracy on Best-fit Parameters and Possible Corrections
It has been previously reported that the accuracy/precision of the maximum likelihood estimate of DRW parameters is sensitive to the "quality" of the  Figure 3. The DHO fits that are identified as good using the algorithm described at the beginning of Section 4.1 are color-coded by the difference between their best-fit log(ω0) and the input. Identified bad fits that were simulated as overdamped DHOs are shown using gray dots and those simulated as underdamped DHOs are shown using magenta crosses. The orange dashed line marks a TQPO of one day. The cyan dash-dotted line shows the cut used to remove the bad fits associated the daily observing cadence (surrounding the orange dashed line). The red contour outlines the distribution of the input DHO parameters.
light curves, more specifically, the ratio between the DRW decorrelation timescale (τ DRW ) and the light curve length (Max ∆t) (Koz lowski 2017, 2021). DHO and DRW are the same class of stochastic diffusion processes-CARMA, thus we expect similar trends to appear for best-fit DHO parameters determined by maximum likelihood (or maximum a posterior with wide flat priors). The light curve "quality" measurements used in (Koz lowski 2017) depend on both the true parameters of the underlying process (e.g., τ DRW ) and some basic properties of the light curve data (e.g., Max ∆t), where the former are hardly known but the latter are easily measurable. Thus, rather than trying to fully characterize the accuracy/precision of our best-fit DHO parameters with respect to those "quality" metrics, we attempted to calibrate the best-fit DHO parameters using simulated data. The same simulation process as de-scribed at the end of Section 4.1 was carried out for all five bands, and the best-fit parameters from different bands were joint together to determine the correction. We first modeled the offset of the best-fit DHO parameters on simulated light curves relative to the input parameters (log(X Output )−log(X Input )) as a multivariate linear function of three basic properties of the light curves: the total length of the light curve (Max ∆t), the minimum separation between any two observations (Min ∆t), and the median photometric uncertainty. The coefficients of the best-fit linear regression suggest that the offset is most correlated with the median photometric uncertainty and Max ∆t-consistent with the results from previous work for DRW (Koz lowski 2017).
By applying a correction derived from the best-fit multivariate regression, we were able to remove the systematic offset between the input and the output distribution of DHO parameters. The top six panels of Figure 5 show the comparison between the input parameters and the output (MAP best-fit) parameters before (orange) and after (blue) the applied correction; the middle six panels show the distributions of the offsets (log(X Output )−log(X Input )). The corrected fits have better correspondence with the input parameters than the original ones (i.e., the blue histograms in the second row are more centered at zero than the orange ones). The panels on the third row of Figure 5 show the offsets of the corrected best fits as a function of the input parameters; any parameter-dependent trends as shown were not accounted for in our correction and are left for future investigations. From simulations, we can see that given the light-curve cadence and photometric accuracy of S82 quasars, σ DHO , τ perturb , and σ are best constrained among all other DHO features-in terms of the size of the dispersion and the level of the parameterdependent trend (modulo the extremes) of the offsets.
In addition to the overestimation/underestimation of DHO parameters, it would be interesting to investigate how reliable is the DHO subclass classification (underdamped DHO vs. overdamped DHO) given the lightcurve cadence and photometric accuracy of S82 quasars. Using the simulated data introduced above, we computed the classification precision (and recall) for both populations. Given a clean sample, for a particular classification, the precision is defined by the percentage of best-fit DHOs that are also simulated with the same classification (e.g., a best-fit overdamped DHO is also simulated as an overdamped DHO) and the recall is defined by the percentage of simulated DHOs that are correctly classified. The final result is shown in Table 1. Since the classification of underdamped DHOs into QPOs and oscillatory DHOs is based on an empir- ically chosen diving T dQPO , it is not as informative to show the statistics for those two classes. In summary, given the temporal sampling and photometric accuracy of SDSS Stripe 82 quasar light curves, more than 99.9% of the overdamped DHOs have the correct classification across all five SDSS bands. Among those classified as underdamped DHOs, the precision is between ≈9% and ≈29% with the highest in the g and r bands and the lowest in the z band. The recall for the overdamped population is between ≈88% and ≈98% with the highest in the g and r bands and the lowest in the z band, and for the underdamped population all bands have a recall of 100% except for the r band (92.68%). Since failed overdamped DHO fits can end up in the region of the underdamped population and failed underdamped DHO fits tend to stay close to their input (see grey dots and magenta crosses in Figure 4), we suspect that the low classification precision of the underdamped population is partially caused by the overabundance of the overdamped population. That said, given that our input parameters were drawn from the distribution of MAP best fits obtained on real SDSS Stripe 82 quasar light curves, the true relative abundance and classification precision of underdamped DHOs should be smaller than what we are showing in Figure 3 and Table 1.

CORRELATIONS OF DHO PARAMETERS WITH WAVELENGTH AND PHYSICAL PROPERTIES OF AGN
We investigated the correlations between the best-fit DHO parameters (corrected using the linear coefficients determined in Section 4.2) and the estimated physical properties of the quasars in our sample (Shen et al. 2011;Rankine et al. 2020). Our investigation was focused on the overdamped population provided that the classification for underdamped DHOs is highly unreli-able (see Table 1). We also limited our analysis to three DHO features: σ DHO , σ and τ perturb , because they are best constrained given the light-curve cadence and photometric accuracy of SDSS Stripe 82 quasars (see Figure 5) and that they are least affected by our initial timescale-based cut for removing bad DHO fits (the first criterion listed at the start of Section 4.1). More specifically, the distribution of τ perturb lies between Min ∆t and Max ∆t, thus, is not affected by the cuts associated with those two timescales. On the other hand, τ rise is on the scale of Min ∆t (≈ 1 day) and τ decay is on the scale of 1 2 Max ∆t (≈ 1000 days), therefore the distribution of τ rise and τ decay are subject to selection bias, and similarly for τ decorr ≈ π 2 (τ rise + τ decay ). Lastly, through examining the MCMC samples of our good DHO fits we found that some of them have bi/multi-modal and/or very broad posterior distributions (see Appendix A for examples). We suspect this to be caused by the irregular/sparse sampling of S82 light curves that leaves certain timescales less well probed and/or the (lack of) variability of strong emission line(s) in the particular photometric bands. Further investigations are needed to verify the origin(s) of those bi/multi-modal posterior distributions. Nevertheless, we removed those objects from the following analysis using a well-defined metric (see Appendix A).

Wavelength Dependence of DHO Parameters
Given the large redshift range that our quasar sample spans and the five different photometric bands that we are utilizing in this investigation, best-fit DHO parameters should first be evaluated against and corrected for any wavelength dependence before being used to correlate with the physical properties of quasars. We explored the wavelength dependence of σ DHO , τ perturb , and σ by first correcting them for redshift dependence and then plotting them as a function of the effective wavelength of their photometric bands in rest frame (Schneider et al. 1983). σ DHO , τ perturb , and σ scale with redshift as (1+z) 0 , (1+z) −1 , and (1+z) 3/2 , respectively; the derivation for the redshift dependence of σ can be found in Appendix B. The effective wavelength, computed based on a power-law continuum with a spectral index α v of −0.5, of the SDSS ugriz bands in the observer's frame are 3541, 4653, 6147, 7461, and 8904 angstroms, respectively (Fukugita et al. 1996;Richards et al. 2001;Kaczmarczik et al. 2009). Figure 6 shows σ DHO , τ perturb , and σ as a function of the effective wavelength of their photometric bands in rest frame (λ RF ), the color indicates the specific passband from which the best-fit parameters are obtained. λ RF for a given band is computed as: λ RF = λ eff /(1+z),  Figure 6. Distribution of σDHO, τ perturb , and σ as a function of rest wavelength (λRF). The color gives the photometric band from which the best-fit parameters are obtained. The black points in the middle and bottom panels mark the median values in bins centered at the corresponding λRF; the error bars have been made twice as large for better visibility. Top: The contours show the 85% mass for each band. σDHO follows an overall monotonic relation with λRF; the solid line shows the best linear fit. Middle: τ perturb increases with λRF longward of ≈2500Å; the solid line shows the best linear fit to the median values with λRF > 2500Å. Bottom: σ decreases with λRF shortward of ≈2500Å; the solid line shows the best linear fit to the median values with λRF < 2500Å.
where λ eff is the effective wavelength in the observer's frame and z is the redshift of an object. From the top panel we can see that σ DHO follows an overall monotonically decreasing trend with λ RF -in agreement with previous findings (e.g., Vanden Berk et al. 2004;MacLeod et al. 2010). We performed a bisector linear regression to derive the best-fit relation between σ DHO and λ RF , log(σ DHO ) = (−0.73 ± 0.070) log(λ RF ) +1.69 ± 0.003.
Note that σ DHO follows an increasing rather than a decreasing trend with λ RF within each individual band. We suspect this "misleading" increasing trend to be caused by the anti-correlation of σ DHO with L bol (see Section 5.2) and the selection bias intrinsic to fluxlimited samples (i.e., selected objects from higher redshift are systematically more luminous than those from lower redshift).
The middle panel of Figure 6 shows that τ perturb (the characteristic timescale of the perturbation component) increases with rest wavelength at λ RF > 2500Å with a best-fit relation of log(τ perturb ) = 0.99 * log(λ RF ) − 2.89, but is nearly independent of wavelength at λ RF < 2500Å.
In the bottom panel of Figure 6, we see that σ (the amplitude of the driving white noise) decreases with an increasing λ RF until ≈2500Å and then becomes nearly independent (or only weakly dependent) of wavelength. At λ RF shorter than 2500Å, the best-fit relation between σ and λ RF is log(σ ) = −1.39 * log(λ RF ) + 3.21.
We note that the trends of σ DHO , τ perturb , and σ with λ RF shown in Figure 6 persist when we replace the corrected DHO parameters with their un-corrected version, however, the v-shaped correlation shown in the bottom panel for σ will appear less obvious.
Both τ perturb and σ exhibit clearly different dependencies with wavelength on either side of 2500Å while σ DHO only shows a largely monotonic dependency. We will further evaluate this interesting feature in Section 6.3 and discuss what could be implied regarding accretion disk models.

DHO Amplitude
The variability amplitude has been studied most extensively among all other variability metrics in terms of searching for correlations with the fundamental properties of AGN (e.g., luminosity) (Vanden Berk et al. 2004;Wilhite et al. 2007;MacLeod et al. 2010;Zuo et al. 2012). It is worthwhile to check if σ DHO exhibits similar correlations with those fundamental properties of the quasars in our sample. In this part of our analysis, σ DHO is corrected to a rest-frame wavelength of 2500Å using Equation 13, and for each quasar only the DHO fit from the photometric band with the smallest uncertainty in σ DHO (determined from MCMC) is used. We also limit our quasars to 0.7 < z < 1.9 for the sake of enforcing a consistent determination of L bol , M BH , and L/L edd (i.e., L 3000 was used to estimate L bol and the Mg II emission line was used to estimate M BH ; Shen et al. 2011). Indeed, we confirm that σ DHO is anti-correlated with the recorded values of L/L edd and M BH for our quasars.   tude. In addition to the expected anti-correlation between EQW [O III] and R Fe II , we see a trend of increasing σ DHO toward the top left corner of this plot-consistent with the result from an earlier investigation conducted using a smaller sample (Ai et al. 2010). The EV1 sequence has long been argued to be driven by the diversity in L/L edd (Boroson 2002;Shen & Ho 2014), that is, data points at the top left corner should have smaller L/L edd than those at the bottom right corner. Indeed, this argument is consistent with that implied by the anticorrelation between σ DHO and L/L edd . At high redshift, C IV EQW and C IV blueshift alone are proposed to be indicators of L/L edd (Shemmer & Lieber 2015;Rankine et al. 2020). The right panel of Figure 8 shows the distribution of our sample in the C IV parameter space; we used the same technique as utilized in the EV1 analysis to bin our data. We found that our quasars tend to have larger variability amplitude than average when C IV EQW is large and C IV blueshift is small (the top left corner); on the other hand, our quasars appear to have smaller amplitude than average when C IV EQW is small and C IV blueshift is large (the bottom right corner). This trend is consistent with the anti-correlation between σ DHO and L/L edd and the suggestions that C IV EQW and C IV blueshift are L/L edd indicators (Shemmer & Lieber 2015;Rankine et al. 2020). A similar trend was also found by Rivera et al. (2020) using multi-epoch spectroscopy, where a hybrid metric combining C IV EQW and C IV blueshift was defined to locate quasars along this trend. This hybrid metric was later referred to as the "C IV distance"  in Richards et al. (2021) and Rivera et al. (2021)-the C IV distance increases its value going from the top left corner to the bottom right corner along the distribution of quasars in the C IV parameter space. Thus, an anti-correlation between σ DHO and C IV distance should be expected.

The Perturbation Parameters: τ perturb and σ
In the DHO framework, τ perturb and σ characterize the driving perturbation to the modeled dynamical system, where σ gives the amplitude of the driving white noise and τ perturb describes a characteristic timescale beyond which the perturbation process loses power (see the bottom left panel of Figure 3 for a reference of the perturbation PSD). In the context of AGN variability modeling, we could expect τ perturb and σ to capture the characteristics of the physical mechanisms that drive the observed UV/optical variability, which might correlate with the fundamental properties of AGN. To that end, we explored the evolution of L/L edd , L bol , and M BH across the τ perturb , σ parameter space. Since τ perturb and σ do not increase/decrease monotonically with λ RF , our analysis only used objects with 3.35 < log(λ RF ) < 3.45, where both τ perturb and σ show clear linear dependence on λ RF so that we can calibrate them to the rest wavelength of 2500Å. As with Section 5.2, we restrict the sample to 0.7 < z < 1.9. The top row of Figure 9 shows the selected objects (0.7 < z < 1.9 and 3.35 < log(λ RF ) < 3.45) in the τ perturb , σ parameter space, and each panel, from left to right, colors objects by their L bol , M BH , and L/L edd , respectively. In the top left panel, we can see a trend of decreasing L bol toward the top right corner (large τ perturb and large σ ); no apparent evolution of M BH or L/L edd across the parameter space can be found in the top middle or right panels.
Since the rest-frame τ perturb is on the scale of days-comparable to the light-crossing time associated with the size of the accretion disks of our quasarsand since the characteristic radius for emission at fixed wavelengths scale linearly with M BH in the log space (Shakura & Sunyaev 1973), it is logical to convert τ perturb into a distance scale (r perturb ) expressed in terms of gravitational radius.
Here, c is the speed of light, G is the gravitational constant, and M BH is the mass of the SMBH associated with a quasar. Thus, r perturb is essentially the mass-weighted version of τ perturb . The bottom row of Figure 9 shows the distribution of the selected quasars in the r perturb , σ parameter space. As with the top three panels of Figure 9, from left to right, quasars in each panel are colored by their L bol , M BH , and L/L edd , respectively. As a guide, we plotted the best-fit regression (dashed lines) in the form of: log(r perturb ) = A * log(σ ) + B * log(M BH ) + C for M BH of 10 8 M , 10 9 M , and 10 10 M . Note that the scatter seen in the top middle panel is caused by objects with different M BH having the same τ perturb , weighting τ perturb by M BH (bottom middle panel) enables us to better see the anti-correlation of τ perturb with σ and facilitate the comparison of different objects on the same physical scale. A new trend of increasing r perturb with increasing L/L edd is also revealed in the bottom right panel of Figure 9. This new correlation has a nonparametric Spearman rank-order correlation coefficient of 0.35 with a two-tailed p-value of 10 −20 .
To further elucidate the contributions of L bol , M BH , and L/L edd to the diversity of perturbation parameters shown in Figure 9, we selected a subsample of quasars with −1.1 < log(L/L edd ) < −0.9 and plotted their distribution in Figure 10 and colored them by their L bol . From Figure 10 we can see that this subsample spans almost the full range of the distribution shown in the bottom three panels of Figure 9 and exhibits a much clearer anti-correlation (smoother color gradient) between the perturbation parameters and L bol (compared to the bottom left panel in Figure 9). The large span of this subsample in the perturbation parameter space and the cleaner/smoother color gradient of L bol suggest that most of the diversity in r perturb (and τ perturb ) is driven by L bol (or M BH ) and that L/L edd works independently of L bol (or M BH ) in terms of determining the observed distribution and plays only a minor role. 6. DISCUSSION

Underdamped DHOs: QPOs and Oscillatory DHOs
In Section 5, we omitted QPOs and oscillatory DHOs from comparing DHO features with the derived physical properties of our quasars because of their low classification precision (see Table 1). However, both subclasses exhibit interesting signatures that can be expected from real physical systems.
As the name suggests, QPOs vary quasi-periodically (see the simulated light curves in Figure 3). Such signals in quasars might be expected from super-massive black hole binaries orbiting each other closely and appearing as a single point source in the image (Begelman et al. 1980). SMBH binaries are expected to emit low-frequency gravitational waves when they merge. QPOs are most likely candidates for those systems and therefore will provide a large pool of potential sources for current and future low-frequency gravitational wave projects (Hobbs 2013;McLaughlin 2013;Amaro-Seoane et al. 2017). At the same time, persistent quasi-periodic oscillations can also be expected from single SMBHs where the inner accretion flow is geometrically thick and undergoes Lense-Thirring precession (Ingram et al. 2009;Graham et al. 2015). Oscillatory DHOs feature larger τ perturb and smaller σ compared to overdamped DHOs (see the top right panel of Figure 3). If most of the variability in overdamped DHOs can be attributed to X-ray reprocessing (which will be discussed further in Section 6.2), then the variability revealed in oscillatory DHOs might originate in the local accretion disk. More specifically, the larger τ perturb and smaller σ could be interpreted as characteristics of a perturbation mechanism different than X-ray illumination, e.g., changes in the mass accretion rate, which might features a longer characteristic perturbation timescale (τ perturb ) and a smaller short-term variability amplitude (σ ) (Pereyra et al. 2006;Li & Cao 2008;Arévalo et al. 2008). Such variability signatures can be expected from quasars with extremely high L/L edd where the radiation from the X-ray corona is weak relative to the intrinsic radiation from the disk (Kubota & Done 2018) or blocked by a puffed up inner disk (see Figure 15 from Leighly (2004) or Figure 18 from Luo et al. (2015) for a reference), both would lead to little or no X-ray reprocessing.
Besides these extreme objects, we might expect to see a mixture of signatures from both an oscillatory DHO and an overdamped DHO in real light curves, that is, assuming that intrinsic disk variability due to changing mass accretion rate and X-ray reprocessing are contributing comparably to the observed variability (Arévalo et al. 2008). Such scenarios could also lead to bi-modal distributions in the MCMC samples (see Appendix A). In those cases, we will need more advanced modeling tools to decouple the light from different processes. Nonetheless, we have removed objects with bi/multi-modal posterior distributions from our analysis and the analysis presented in Section 5 is concentrated on the overdamped population only, therefore, it is logical to suspect that the variability features revealed in this investigation are likely dominated by one single mechanism (e.g., X-ray reprocessing).

Long-term Amplitude of Overdamped DHOs:
Primarily Determined by L/L edd ?
We find that the long-term asymptotic variability amplitude of AGN in the overdamped subclass, as characterized by σ DHO , is anti-correlated with Eddington ratio and black hole mass. This finding is consistent with that resulted from previous investigations utilizing other methods (Vanden Berk et al. 2004;Wilhite et al. 2007;MacLeod et al. 2010;Ai et al. 2010;Simm et al. 2016). The anti-correlation of σ DHO with L/L edd is expected in a model where the size of the hot X-ray corona relative to the accretion disk anti-correlates with L/L edd and the optical variability is largely due to reprocessing of X-ray photons (e.g., Kubota & Done 2018;Giustini & Proga 2019). More specifically, a high L/L edd corresponds to a small X-ray corona and a large/strong disk, therefore, less reprocessing of X-ray photons occurs in the disk-leading to a smaller long-term variability amplitude; on the other hand, a low L/L edd indicates a large X-ray corona relative to the disk, thus, more reprocessing of X-ray photons occurs in the disk and a larger long-term variability amplitude can be expected. We note that the correlations discussed in this section can also be produced by changing mass accretion rate in the disk (Li & Cao 2008), moreover, recent work on intensive multi-band AGN reverberation mapping for a handful of objects have provided potential evidence for such scenario (Edelson et al. 2017;Starkey et al. 2017;Edelson et al. 2019;Cackett et al. 2020), see Cackett et al. (2021) for more discussion on this topic.
If L/L edd sets the basic level of AGN variability, then the additional anti-correlation with M BH can be explained by the recognition that the part of the accretion disk that emits at a fixed effective temperature changes radius with increasing/decreasing M BH . The flux emitted per unit area at a radius R >> 6R g on the disk is defined as where T is the effective temperature at a radius R,ṁ is the Eddington ratio, M BH is the black hole mass, and R g is the gravitational radius (Novikov & Thorne 1973). According to this equation, the radius (R) on the disk that emits at a fixed temperature (a fixed wavelength), given a known Eddington ratio, moves outwards with increasing M BH resulting in a larger distance to the illuminating X-ray corona. In a 'lamppost' geometry, the larger the distance to the X-ray corona, the smaller the intensity of the illuminating X-ray photons and the smaller the variability amplitude due to X-ray reprocessing. Therefore, we argue that the variability amplitude could be largely determined by L/L edd with M BH acting as a secondary parameter.

Perturbation Parameters of Overdamped DHOs:
Indicating a Two-component Accretion Disk?
In Figure 6, we saw that τ perturb and σ do not strictly follow a power-law relation with λ RF with a break point at ≈2500Å; the change in dependency with λ RF could signify that two different physical processes are involved in shaping the observed variability. Our speculation can be explained by a physical model of the accretion flow consisting of: a hot X-ray corona extending from the innermost stable circular orbit (ISCO) to the inner edge of a truncated disk (r hot = R hot /R g ), a warm Comptonization region (for producing the soft X-ray excess) spanning from r hot to an intermediate radius (r warm = R warm /R g ), and a standard cold thin disk going from r warm to the outer edge of the disk (RóżaŃska & Czerny 2000;Czerny et al. 2003;Sobolewska et al. 2004;Done et al. 2012;Kubota & Done 2018). In this model, the warm Comptonization region features two slabs of warm electrons sandwiching the standard disk; it is also assumed that the hot X-ray corona is the main driver of the UV/optical variability. A schematic view of said geometry is shown in Figure 11, where the blue annulus represents the hot X-ray corona, the green slab represents the warm Comptonization region, and the purple slab corresponds to the outer standard disk. Note that the disk does not extend to the ISCO and truncates at r hot .
More concretely, given this two-component accretion disk model (warm Comptonization region + standard thin disk region) and the fact that the effective temperature of a disk annulus scales with its radius (see Equation 18), we might suspect the observed break point (at ≈2500Å) in the wavelength dependence of τ perturb and  Figure 11. Diagrams providing schematic views of the accretion disk geometry proposed in Kubota & Done (2018) for a variety of L/L edd and MBH. Note that the diagrams are drawn in the units of gravitational radius and Equation 18 is used as a reference when discussing the changing disk geometry. The blue annulus represents the hot X-ray corona, the green slab shows the warm Comptonization region of the disk, and the purple slab corresponds to the standard disk region. The effective emitting wavelength of the disk increases going from r hot to rout. The r perturb arrow connects the X-ray corona to the disk annulus emitting at an effective wavelength of 2500Å. Left column: MBH is held fixed and L/L edd increases from top to bottom. The size of the X-ray corona decreases and the 2500Å mark shifts outwards with increasing L/L edd . The effective temperature/wavelength of the disk at the boundary between the warm Comptonization region (green slab) and the standard disk region (purple slab) also increases/decreases with L/L edd . Right column: L/L edd is held fixed and MBH increases from top to bottom. When L/L edd is fixed, the disk temperature (wavelength) at the boundary between the green and purple slab should stay relatively fixed, thus, increasing MBH shifts both rwarm and the 2500Å disk annulus inwards. σ to correlate with the transition point from the warm Comptonization region to the standard disk region. At r > r warm (or λ RF > 2500Å given our data), the hot X-ray corona directly illuminates the disk, and τ perturb probes the light-crossing time between the emitting disk annulus and the corona. At r < r warm (or λ RF < 2500Å given our data), the hot X-ray photons are scattered by the warm electrons in the upper atmosphere of the disk, therefore, the perturbation process exhibits an increased τ perturb from that expected for direct illumination (see the best-fit line and the median values in the middle panel of Figure 6). At the same time, the down-scattered photons (with reduced frequencies) are better absorbed at the disk surface compared to the source X-ray photons because the absorption coefficients for both freefree and bound-free adsorptions are inversely proportional to frequency (Rybicki & Lightman 1986), which leads to an intensified perturbation amplitude (σ ). We acknowledge that the empirically identified break wavelength of 2500Å is an ensemble average over the true break wavelengths for our selected quasars and that the break wavelength could shift as a function of L/L edd (see the left column of Figure 11), as suggested by Kubota & Done (2018). However, the size of the error bars in our inferred parameters prohibit us from revealing it. Future investigations utilizing better light curves are needed to characterize that dependency.
The trend revealed by the bottom panel of Figure 6 for σ and the change in behaviors at ≈2500Å also agrees with the result presented in Wilhite et al. (2005), which used a composite difference spectrum constructed from ≈300 SDSS quasars to demonstrate a similar λ RF dependency of σ . However, such a trend is not apparent in the correlation of the long-term variability amplitude with λ RF -as shown with σ DHO and other metrics (e.g., σ DRW ; MacLeod et al. 2010). We suspect the missing imprint of the short-term variability amplitude (σ ) on the long-term variability amplitude (σ DHO ) could be related to the thermalization of the illuminating photons with the disk over long timescales. It is also possible that the correlation of the long-term variability amplitude with λ RF is itself a function of L/L edd and/or L bol and that the quasar samples used in previous work and this one for deriving such correlation span a large range of L/L edd and L bol . Indeed, the quasars used in this work cover a slightly larger range, as characterized by the dispersion of the distribution, of L/L edd than that used in Wilhite et al. (2005). The L/L edd distribution of our quasars have a median absolute deviation (MAD) and a IQR of 0.54 dex and 0.27 dex, respectively, whereas the same statistics for quasars used in Wilhite et al. (2005) have a value of 0.43 dex and 0.22 dex, respectively. Moreover, we found σ DHO to follow a similar v-shaped trend with λ RF as for σ when binned with L bol or L/L edd . A larger sample of quasars with trusted physical properties and high-quality light curves is needed to investigate and characterize how the correlation of σ DHO with λ RF depends on the value of L bol and L/L edd . Based on the model proposed in Kubota & Done (2018) and the fact that we can approximate the effect of the illuminating X-ray corona at r >> 6 using a point source on the spin axis at a height of H = r hot (Gardner & Done 2017), r perturb at the break wavelength is effectively r 2 hot + r 2 warm , assuming it probes the distance from the X-ray corona to the disk at r warm . According to the best-fit parameters for Mrk 509 and PG1115+407 shown in Table 2 of Kubota & Done (2018), r 2 hot + r 2 warm should exhibit a negligible correlation with L/L edd given a fixed M BH . Indeed, we only see a weak trend of increasing r perturb with increasing L/L edd from the bottom right panel of Figure 9. We also suspect that trend as a result of incorrect break wavelength being adopted for quasars having a large range of L/L edd , that is, the break wavelength should decrease with increasing L/L edd and calibrating all τ perturb to the λ RF of 2500Å produces a false-positive trend of increasing r perturb with increasing L/L edd (see the diagrams in the left columns of Figure 11). Nonetheless, the distance scale reflected by τ perturb (r perturb : 30-150 R g ) is comparable to the numbers quoted in Table 2 of Kubota & Done (2018). In addition, for a given L/L edd , r perturb should decrease with increasing M BH (and L bol ), because r warm (if corresponding to a constant temperature) decreases when M BH increases according to Equation 18 (see the right column of Figure 11) and H = r hot is rel-atively fixed for a fixed spin (Kubota & Done 2018); the selected subsample of quasars with −1.1 < log(L/L edd ) < −0.9 as shown in Figure 10 demonstrates this correlation.
We currently do not understand the origin of the anticorrelation between r perturb and σ (at a fixed M BH ), but it is consistent with a picture where the closer the disk is to the X-ray corona the larger the perturbation amplitude. The most straightforward explanation is that the range of spins of the central black holes has produced the diversity in r perturb /τ perturb given all other parameters of the SMBH held fixed; however, further investigation is required to verify this hypothesis.
Our results lead to conclusions: 1) the warm Comptonization region might be responsible for the observed break in the wavelength dependence of τ perturb and σ ; 2) τ perturb at λ RF > 2500Å might be associated with the light crossing time from the X-ray corona to the accretion disk; 3) the short-term variability amplitude decreases with wavelength at λ RF < 2500Å and stays roughy constant (or drops at a much slower rate) with rest-frame wavelength-consistent with that shown in Wilhite et al. (2005).

SUMMARY
In this work, we have investigated the UV/optical variability of ≈12,000 SDSS S82 quasars by modeling their light curves in the ugriz bands as DHO processes. A DHO process can be fully described by four basic parameters: a natural oscillation frequency (ω 0 ), a damping ratio (ξ), a characteristic perturbation timescale (τ perturb ), and an amplitude for the perturbing white noise (σ ). The asymptotic (long-term) variability amplitude is characterized by σ DHO (a function of the above four parameters). We explored the correlations of the best-fit DHO parameters and the derived features with the physical properties of our quasars estimated by Shen et al. (2011). The main results are summarized below: 1. The distribution of the best-fit DHO parameters splits naturally into two main clusters: an overdamped DHO population and an underdamped DHO population (see Figure 3). The overdamped population exhibits similar variability signatures as those characterized by a DRW model. The underdamped population can be further classified (empirically) into the QPO subclass and the oscillatory DHO subclass. Both QPOs and oscillatory DHOs have light curves that are smoother at short timescales (smaller ξ) but bumpier at long timescales than those of overdamped DHOs. The QPO subclass features observable quasiperiodicity in its light curve.
2. σ DHO , τ perturb , and σ of overdamped DHOs evolve with rest-frame wavelength (λ RF ). σ DHO follows a largely monotonic trend with λ RF with a powerlaw index of −0.73 ± 0.070. τ perturb also follows a power-law relation with λ RF but only at λ RF > 2500Å. For σ , a clear power-law relation is only observed at λ RF < 2500Å.
3. After correcting for the wavelength dependence, σ DHO (of overdamped DHOs) exhibits anticorrelations with both L/L edd and M BH -in agreement with the results of previous work (e.g., Wilhite et al. 2007;MacLeod et al. 2010). However, our best-fit regression (Equation 14) suggests steeper anti-correlations of σ DHO with L/L edd and M BH than that reported in MacLeod et al. (2010) for a DRW model. 4. We found that the short-term variability of overdamped DHOs as parameterized by τ perturb and σ is anti-correlated with L bol .
5. We argue that the different characteristics of AGN UV/optical variability revealed by overdamped DHOs as listed above can be connected together using a physical picture: 1) L/L edd determines the size (extension) of the X-ray corona relative to the truncated accretion disk; 2) the majority of the observed variability can be attributed to the reprocessing of X-ray photons by a relatively passive accretion disk; 3) from inside out, the accretion disk is divided into two regions-a warm Comptonization region and a standard thin disk region, where the warm Comptonization region is sandwiched by layers of warm electrons that produce the soft X-ray excess (Kubota & Done 2018).
Limited by the cadence and photometric accuracy of the dataset used in this work, we are only able to investigate the correlations of three features of the overdamped DHO population with the physical properties of AGN. With light curves of better temporal sampling and/or higher photometric accuracy such as those from current/future time-domain surveys (Chambers et al. 2016;Bellm et al. 2019;Ivezić et al. 2019), we can further exploit DHO modeling in AGN variability study. For example, the QPO subclass can provide a pool of SMBH binary candidates, when utilized jointly with other SMBH binary discovery methods, both the selection efficiency and completeness can be improved (Graham et al. 2015;Liu et al. 2016;Charisi et al. 2016). In addition, oscillatory DHOs are likely candidates for AGN that accrete at extremely high L/L edd , where the inner disk has puffed up to prevent the X-ray corona from directly illuminating the disk thus resulting in a smoother light curve at short timescales (Leighly 2004;Luo et al. 2015). Moreover, the discovered anticorrelation of τ perturb and σ with L bol , together with the known anti-correlation of σ DHO with L/L edd and M BH (or L bol ), can be utilized to develop new algorithms that will derive L/L edd for millions of AGN using photometric data alone. Lastly, our sample of quasars span only a limited range in M BH (10 8 M to 10 10 M ), applying the same analysis technique used here to AGN of much smaller M BH (e.g., 10 6 M ) will help further elucidate the correlations between variability and fundamental properties of AGN.
Despite the better sampling and high S/N of light curves coming from current and future time-domain surveys (Bellm et al. 2019;Ivezić et al. 2019), we stress the need to develop reliable methods that can effectively merge light curves from multiple surveys for the sake of extending the baseline. Given the current likelihoodbased inference technique, the long-term variability of AGN (on the scale of years) can only be best constrained when the light curves are much longer than the intrinsic timescales (Koz lowski 2017, 2021). However, the stochastic nature of AGN variability and its correlation with wavelength make merging light curves by median/mean magnitude a sub-optimal solution. In addition to constructing longer light curves, better inference algorithms that can specifically tackle the effects of sparse sampling and short baseline for light curves will be very helpful. We also emphasize that an ultimate algorithm that can efficiently fit light curves from different bands simultaneously will be essential to perform a 'true' CARMA modeling of AGN light curves given that the light curves from different passbands (and the information encoded therein) should be inter-correlated (see Hu & Tak (2020) for a heuristic example of such approach).
Last but not least, stochastic diffusion processes like DRW and DHO are statistical models rather than physical models. Therefore, care should be taken when interpreting the timescales extracted from stochastic modeling. One interesting discovery that we made while comparing the DRW features and the DHO features derived from the dataset used in this work is that the decorrelation/decay timescale of DRW (τ DRW ) is not correlated with the decorrelation/decay timescale of DHO (τ decorr or τ decay ) as we would have expected, but instead exhibits a tight correlation with the ratio between τ decay and τ perturb of DHO (see Figure 14). We discuss plausible origin(s) of this "mis-match" in Appendix C. By visually examining the posterior distribution of our MCMC samples, we found that some of them appear to be bi/multi-modal (see the left panel of Figure 12). We suspect this bi/multi-modality to have a mixture of different origins: intrinsic degeneracy (observed variability coming from different physical processes with comparable contributions), the poor sampling of our light curves, and strong emission lines landing in the range of the specific photometric bands. To automatically remove those fits from our analysis, we designed a bi-modality index S BM , S BM = σ τ perturb + y∈P [log(γ y ) + g(τ y ) * log(τ y /150)], P = [α 1 , α 2 , β 0 , β 1 ] (A1) γ = (p 75 − p 25 )/(p 99 − p 1 ) (A2) where σ τ perturb is the "1-sigma" range of the posterior distribution of DHO's τ perturb in the log scale, τ is the autocorrelation time/step for the MCMC (see Foreman-Mackey et al. (2013)), and p n is the nth percentile of the MCMC samples. Here, γ measures the tailness of the marginalized distribution. S BM penalizes for large γ (small tails) and large τ (MCMC takes more steps to converge), both of which are indicators of potential bi/multi-modal distributions; at the same time, no incentives are given for small τ (implemented through g(τ ) as an "activation function"). For a Gaussian distribution, γ has a nominal value of ≈0.3. We also found that the MCMC chains for objects with single-mode posterior distributions (see the right panel of Figure 12) usually converge within 150 steps. σ τ perturb plays the role of identifying unconstrained posteriors. The distribution of σ τ perturb shows a natural break at 1.4, which is reasonable given that the cadence of our light curves is only sensible on timescales from 1 day to 2000 days-covering approximately 3 dex. That said, DHOs with σ τ perturb greater than 1.4 (a range of 2.8 dex) should be considered unconstrained. These realizations lead to a selection cut of S BM > 2.6, where such DHO fits were not used for the analysis presented in Section 5. If the bi/multi-modal distribution is solely caused by emission lines contaminating the continuum variability, we should expect to see a higher percentage of objects being removed by our cut when the rest-frame effective wavelength is around an emission line or the Balmer Continuum (BC), which is also expected  Figure 13. The percentage of DHO fits removed by a cut of SBM > 2.6 as a function of the rest-frame effective wavelength of each photometric band. The shaded areas show the rest-frame coverage of the photometric bands centered at the C IV and Mg II emission lines. Some local peaks can be spotted at the wavelengths of emission lines (or within the shaded regions), e.g., at the C IV line for g and r bands, and at the Mg II line for r and z bands.
to be external to the disk. We investigated this expected effect and show the result in Figure 13. From Figure 13, we can see that in each band there are peaks occurring around the wavelengths of selected emission lines (or BC).
Although not every emission line in every band has a matching local peak and the peaks are not as prominent, this is still a supportive evidence for our speculation that emission lines might contribute to the bi/multi-modal posterior distribution shown in Figure 12, at least partially. Emission line variability being a likely source of contamination for our analysis would be a good indicator of the promise of photometric reverberation mapping with the upcoming Rubin C. Observatory Legacy Survey of Space and Time (Chelouche & Daniel 2012;Chelouche et al. 2014;Ivezić et al. 2019 ).

C. OVERDAMPED DHO VS. DRW
The overdamped DHO model is similar to a DRW, but is more flexible. It would be interesting to compare the variability information extracted by the two models side by side using a same dataset. DRW can be generally characterized by an asymptotic amplitude (σ DRW ) and a long-term decay timescale (τ DRW ). The features extracted by DHO that are comparable to those two DRW parameters are σ DHO and τ decay . We fitted both models to the g-band light curves of our quasars. We compared σ DHO with σ DRW , σ (the β 0 coefficient in Equation 1 and Equation 2) of DHO with that of DRW, and τ decay of DHO with τ DRW ; the results are shown in the top left, top right, and bottom left panel of Figure 14, respectively. From those three comparisons, we only found a good one-to-one correlation between σ DHO and σ DRW . One would have expected τ decay of DHO to correlate with τ DRW , since both characterize how long it takes for the modeled system to forget about its past self, however, we did not see that. Instead, we found a strong correlation between the ratio of τ decay to τ perturb and τ DRW (see the bottom right panel of Figure 14). The slight "overestimation" of τ decay /τ perturb at large τ DRW can be explained by the well known underestimation of τ DRW when the light curve span is shorter than 10 times the intrinsic τ DRW (Koz lowski 2017), which is exactly what is being shown. We are not certain about the origin of this correlation, but if τ perturb is physical then we suspect that the DRW model might be overlooking this information. More specifically, since DRW requires the perturbation process to be a white-noise with a flat power spectrum, then τ DRW might be in units of τ perturb because the DHO perturbation process only behaves like a white-noise at a timescale longer than τ perturb (see the bottom left panel of Figure 3 for the power spectrum density of the DHO perturbation process). In addition, studies utilizing the power spectrum density technique have also suggested a potential second short-term characteristic timescale (relative to the known long-term decorrelation timescale) in AGN light curves, where the reported short-term timescales are comparable to the τ perturb revealed in this work (Zu et al. 2013;Stone et al. 2022). We stress that the correlation shown between τ decay /τ perturb and τ DRW does not necessarily invalidate a DRW (or DHO) description of AGN variability, but argues for more careful examinations of the timescales returned by such modeling.