Improving Damped Random Walk parameters for SDSS Stripe 82 Quasars with Pan-STARRS1

Krzysztof Suberlak

0000-0002-9589-1306]Krzysztof L. Suberlak
Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195, USA

0000-0001-5250-2633]Željko Ivezić
Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195, USA

0000-0003-3422-2202]Chelsea MacLeod
Harvard Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

We use the Panoramic Survey Telescope and Rapid Response System 1 Survey (Pan-STARRS1, PS1) data to extend  the Sloan Digital Sky Survey (SDSS) Stripe 82 quasar light curves. Combining PS1 and SDSS light curves provides a 15 yr baseline for 9248 quasars - 5 yr longer than prior studies that used only SDSS. We fit the light curves with the damped random walk (DRW) model model - a statistical description of their variability. We correlate the resulting DRW model parameters: asymptotic variability amplitude SF$_{\infty}$, and characteristic timescale $\tau$, with quasar physical properties - black hole mass, bolometric luminosity, and redshift. Using simulated light curves, we find that a longer baseline allows us to better constrain the DRW parameters. After adding PS1 data, the variability amplitude is a stronger function of the black hole mass, and has a weaker dependence on quasar luminosity. In addition, the characteristic timescale $\tau$ dependence on quasar luminosity is marginally weaker. We also make predictions for the fidelity of DRW model parameter retrieval when light curves will be further extended with Zwicky Transient Facility (ZTF) and the Rubin Observatory Legacy Survey of Space and Time (LSST) data. Finally, we show how updated DRW parameters offer an independent method of discovering changing-look quasar candidates (CLQSOs). The candidates are outliers in terms of differences in magnitude and scatter between SDSS and PS1 segments. We identify 40 objects (35 newly reported) with tenfold increase in variability timescale between SDSS and SDSS–PS1 data, which is due to a large change in brightness (over 0.5 mag) - characteristic for CLQSOs.

§ INTRODUCTION

Quasars are variable at rest-frame optical wavelengths at an asymptotic rms level of about 0.2 mag. These distant galaxies harbor an actively accreting supermassive black hole: an active galactic nucleus (AGN). Although it is agreed upon that the majority of optical light originates from the thermal emission of the accretion disk, the detailed origin of variability has been debated for the past 50 yr (see Sun et al., 2018 and references therein). Some favor a thermal origin of variability [Kelly et al., 2013] related to the propagation of inhomogeneities (`hot spots') in the disk [Dexter & Agol, 2011, Cai et al., 2016]; others suggest magnetically elevated disks [Dexter & Begelman, 2019] or X-ray reprocessing  [Kubota & Done, 2018].  Indeed, it may well be that the answer involves a combination of these; as [Sánchez-Sáez et al., 2018] suggested, perhaps short-term variability (hours to days) is linked to the changes in X-ray flux, while long-term variability (months to years) is more intrinsic to the disk [Edelson et al., 2015, Lira et al., 2015]. Nevertheless, quasar optical light curves have been successfully described using the damped random walk (DRW) model [Kelly et al., 2009, MacLeod et al., 2010, Kozłowski et al., 2010, Zu et al., 2011, Kasliwal et al., 2015], and the DRW parameters have been linked to the physical quasar properties (MacLeod et al., 2010, hereafter M10).

Variability is also a classification tool, allowing one to distinguish quasars from other variable sources that do not exhibit a stochastic variability pattern [MacLeod et al., 2011]. This property is especially useful for selecting quasars in the intermediate-redshift range, which overlaps the stellar locus in color-color diagrams [Sesar et al., 2007, Yang et al., 2017]. Variability has also been used to increase the completeness in measurements of quasar luminosity function (see McGreer et al., 2013, McGreer et al., 2018, Palanque-Delabrouille et al., 2013, Ross et al., 2013, AlSayyad, 2016).

Power spectral density (PSD) informs us about the distribution of variability across the frequency range, from short timescales (high frequencies) to long timescales (low frequencies). Quasar – or more broadly speaking, AGN– variability, exhibits a broken power-law PSD, of the form
$\log{P(f)} \propto \alpha \log{(f)}$, with $\alpha_l$ at low frequencies and $\alpha_h$ at high frequencies. For a pure DRW process,  $\alpha_{h}{=}-2$ and $\alpha_{l}{=}0$, so that

\begin{equation}
P(f) = \frac{4\sigma^{2}\tau}{1+(2 \pi \tau f)^{2}}
\end{equation}

(with $\sigma = \mathrm{SF}_{\infty} / \sqrt{2}$, $\tau$ the characteristic timescale, and $f$ the frequency), where $P(f) \propto f^{-2}$  at high frequencies $f > (2\pi \tau)^{-1}$, and levels to a constant value at lower frequencies  [Kelly et al., 2014].

There is a debate in the literature about the exact shape of the quasar PSD and any possible  departures from the pure DRW model. Studies using quasar data from wide-field photometric surveys (OGLE, SDSS, PS1) benefit from relatively long baselines (several years), which constrain the low-frequency part of the PSD. Overall, there is no evidence of a significant departure from DRW at these long timescales, i.e. $\alpha_{l} \approx 0$ [Zu et al., 2013, Simm et al., 2016, Kozłowski, 2016, Caplar et al., 2017, Guo et al., 2017, Sun et al., 2018]. However, these ground-based surveys suffer from a sparse sampling, which can be remedied by using a space-based telescope that can carry out near-continuous observations, like the Kepler mission  [Borucki et al., 2010]. Studies using Kepler data that focused on a smaller number of well-sampled AGNs with short baselines (<100 days), found a range of power-law slopes at high frequencies – from -1 to -3.2, which includes the DRW $\alpha_{h} \approx -2$– but further study is needed [Mushotzky et al., 2011, Edelson et al., 2014, Aranzana et al., 2018, Smith et al., 2018]. However, in this paper, the timescales probed are larger than several days; thus, we can assume that DRW is the best working description of quasar variability for available optical light curves. Furthermore, in this work, we directly compare the results of Sloan Digital Sky Survey (SDSS) light curves extended with the Panoramic Survey Telescope and Rapid Response System 1 Survey (Pan-STARRS1, PS1) to M10, who used pure DRW description (see discussion therein on a possible departure from DRW). Therefore, to allow a better comparison of our results with M10 we  use the DRW description of quasar PSD.

Due to its stochastic nature, for an unbiased parameter retrieval of the DRW process, the light curve is required to be several times longer than the characteristic timescale (Kozłowski et al., 2010, Kozłowski, Szymon, 2017, hereafter K17). This means that DRW parameters recovered for short light curves (compared to the recovered timescale) may be biased, which in turn affects the correlations with physical parameters (black hole mass, Eddington ratio, absolute luminosity).

For this reason, while some studies have restricted the probed redshift range, limiting the quasar sample to where one would expect only shorter timescales based on previous studies [Kelly et al., 2013, Simm et al., 2016, Guo et al., 2017, Sun et al., 2018], some have elected not to study timescales at all [Sun et al., 2018, Sánchez-Sáez et al., 2018] or to use the timescales recovered from short light curves primarily for classification [Hernitschek et al., 2016].

By extending the available quasar light curves, we are able to better recover DRW timescales. Since almost a decade ago, when M10 published their study based on SDSS Stripe 82 (S82) data, new data sets (PS1, PTF, CRTS) have become available. They can extend the quasar light curves by almost 50%.  For instance, Li et al., 2018 combined SDSS and Dark Energy Camera Legacy Survey  (DECaLS) data to provide a 15 yr baseline. However, using all SDSS quasars, rather than only those confined to a well-observed S82  equatorial region, meant  that the majority of light curves suffered from a poor sampling (only a few epochs). Therefore, rather than directly fitting individual light curves, they had to resort to an ensemble structure function (SF) approach.

Unlike previous studies, in this work, by combining SDSS and PS1 data for the well-observed S82, we afford an extended baseline (15 yr), a large number (9000) of quasars, and  a good cadence (N > 60 epochs), to which we fit the DRW model. The layout of this paper is as follows. We confirm in Section <ref> that extending the quasar baseline is an important improvement in providing unbiased estimates of the DRW model parameters (K17);  in Section <ref> we describe the data sets employed and their combination into a common photometric system;  in Section <ref>, we simulate the improvement in the recovery of DRW parameters with baseline extension and realistic cadence; in Section <ref>, we describe the main results, analyzing correlations between physical parameters and variability; in Section <ref>, we discuss the physical meaning of relevant timescales; and in Section <ref>, we summarize the main conclusions. In this work, we adopt a $\Lambda $CDM cosmology with $h_{0} = 0.7$ and $\Omega_{m} = 0.3 $.

§ METHODS

§.§ DRW as a GP

The DRW (Ornstein–Uhlenbeck process, Rasmussen & Williams, 2006) can be modeled as a member of a class of Gaussian processes (GPs). Each GP is described by a mean and a kernel: a covariance function that contains a measure of correlation between two points, $x_{n}$ and $x_{m}$, separated by $\Delta t_{nm}$ (autocorrelation). For the  DRW process, the covariance  between two observations spaced by  $\Delta t_{nm}$ is

\begin{eqnarray}
k(\Delta t_{nm}) &=& \sigma^{2}\exp{(-\Delta t_{nm} / \tau)}  \nonumber \\
&=& \sigma^{2} \mathrm{ACF}(\Delta t_{nm})
\label{eq:covariance}
\end{eqnarray}

Here $\sigma^{2}$ is an amplitude of correlation decay as a function of $\Delta t_{nm}$,  while $\tau$ is the characteristic timescale over which correlation drops by $1/e$. For a DRW,  the correlation function $k(\Delta t_{nm})$ is also related to the autocorrelation function.

Not explicitly used in this paper, but of direct relevance to the DRW modeling, is the SF. It can be found from the data as the rms scatter of  magnitude differences $\Delta m$  calculated as a function of temporal separation $\Delta t$ (we drop the subscripts $n$ and $m$ for brevity). The SF is directly related to a DRW kernel $k(\Delta t)$:

\begin{equation}
\mathrm{SF}(\Delta t) = \mathrm{SF}_{\infty} (1-\exp{(-|\Delta t|/\tau)})^{1/2}
\end{equation}

For quasars, SF approximately follows  a power law, $\mathrm{SF} \propto \Delta t^{\beta}$,  and for large time separation $\Delta t$, as epochs in the light curve cease to be correlated, it levels out to a constant value $\mathrm{SF}_{\infty}$: the asymptotic SF.  Note that $\sigma$ in Equation <ref> is related to the asymptotic SF as $\mathrm{SF}_\infty = \sqrt{2} \sigma$ (also see Bauer et al., 2009, MacLeod et al., 2012, Graham et al., 2015 for an overview).

§.§ Fitting

We evaluate the likelihood of the DRW model with a particular set of $\tau,\sigma$ given the data with celerite, a fast GP solver [Foreman-Mackey et al., 2017]. The underlying matrix algebra is similar to that used by [Rybicki & Press, 1992], [Kozłowski et al., 2010], and M10. Also, as in previous work, we use a prior on the DRW parameters that is  uniform  in log space:  $1 / (\sigma \tau)$. The main difference in our approach is that rather than adopting the maximum a posteriori (MAP) as the `best-fit' value for the DRW parameters (as in Kozłowski et al., 2010, K17; Kozłowski, 2016, M10; MacLeod et al., 2011),  we find the expectation value of the marginalized posterior. This is advantageous because of the non-Gaussian shape of the posterior; otherwise, if the posterior was a 2D normal distribution, the expectation value would coincide with the maximum of the posterior (MAP solution).

§.§ The Impact of Light Curve Baseline

It was reported by K17 that one cannot trust any results of DRW fitting unless the light curve is at least 10 times longer than the characteristic timescale. In this section, we revisit the relationship between recovered and input timescales as a function of light curve baseline by following the K17 setup. We confirm that the bias in retrieved DRW timescale depends on how many times the light curve is longer than the timescale. However, we find that the baseline does not have to be as many as 10 times longer to provide meaningful, rather than unconstrained, results. Assuming a fixed baseline of $\Delta T = 8$ yr, we simulate 10,000 light curves, exploring 100 values of input timescales, but identical $\mathrm{SF}_{\infty}=0.2$ mag with either SDSS ($N=60$), or OGLE-like ($N=445$) cadence. Defining $\rho$ as the ratio of input timescale to baseline, we probe a range of $\rho$ between 0.01 and 15, uniform in a logarithmic grid.

Recovery of the input DRW timescale with baseline fixed to $\Delta T = 8$ yr. We explore 100 logarithmically spaced values of $\rho \equiv \tau / \Delta T$, simulating 100 realizations of the DRW process at each $\rho$. The impact of photometric uncertainties and cadence is small in this case: the left panel (SDSS; $N=60$ epochs) does not significantly differ from the right panel (OGLE; $N=445$ epochs). The dotted horizontal and solid vertical lines mark $\rho = 0.1$, i.e. the baseline being 10 times longer than the timescale. The solid diagonal line corresponds to a perfect recovery of DRW parameters (where $\rho_{in} =\rho_{out} $). For any quasar, extending its light curve  moves it toward the top right (well-constrained) portion of the diagram, since for a fixed $\tau_{in}$, increasing $\Delta T$ decreases $\rho$. For baselines shorter than $\tau$, best-fit $\tau$ is underestimated and becomes biased to $\Delta T /2 $.

For each light curve, we simulate the underlying DRW signal $s(t)$ by iterating over the array of time steps $t$.  At each step, we draw a point from a Gaussian distribution, for which the mean and standard deviation are recalculated at each time step (see Equations A4 and A5 in Kelly et al., 2009; Sec. 2.2 in M10; K17). Initially, at $t_{0}$, the signal is equal to the mean magnitude, $s_{0} = \langle m \rangle$. After a time step $\Delta t_{i} = t_{i+1} - t_{i}$, the signal $s_{i+1}$ is drawn from  a normal distribution $\mathcal{N}(loc, stdev)$, with

\begin{equation}
loc = s_{i} e ^ { - r  }  + \langle m \rangle \left( 1 - e ^{ - r }\right)
\end{equation}

\begin{equation}
stdev^{2} =  0.5  \, \mathrm{SF}_{\infty}^{2} \left( 1 - e ^{  - 2 r  }  \right)
\end{equation}

where  $r = \Delta t_{i} / \tau$, and $\tau$ is the damping timescale.

Like K17, we add to the true underlying signal with zero mean $s(t)$ and a mean magnitude  ($\mathrm{r_{SDSS}}=17$ mag and $\mathrm{I_{OGLE}}=18$ mag) and calculate a magnitude-dependent estimate of photometric uncertainty:

\begin{eqnarray}
\sigma_{\mathrm{SDSS}}^{2} &=& 0.013^{2} + \exp{[2 (\mathrm{r_{SDSS}}-23.36)]} \\
\sigma_{\mathrm{OGLE}}^{2} &=& 0.004^{2} + \exp{[1.63 (\mathrm{I_{OGLE}} - 22.55)]}
\end{eqnarray}

To simulate observational conditions, we add the Gaussian noise $n(t) = \mathcal{N}(0,\sigma(t))$:

\begin{equation}
y(t) = s(t) + n(t)
\end{equation}

The resulting distribution of fitted timescales  as a function of input timescales scaled by the  8 yr baseline, $\rho_{out}$, versus $\rho_{in}$, is shown in Figure <ref>. We confirm the findings of K17: for short light curves, the best-fit $\tau$ becomes ${\sim}1/5$ of the light curve length (where $\log_{10}{(\rho_{out})} \approx -0.7$, the `unconstrained' region; bottom left in each panel). However, as long as the light curve is several times longer than the timescale ($1/\rho \gtrapprox 3, i.e. \log_{10}{(\rho)} \lessapprox 0.5$), we can recover the timescale without substantial bias (the dashed line approaches the solid diagonal line in both panels). In summary, in this section, we showed that when using a DRW description for quasar light curves, extension of the light curve baseline moves the fit results from the biased region (Figure <ref>, bottom left of each panel) to the unbiased regime (Figure <ref>, top right of each panel). This is the basis for this study,  in which we extend the baselines of quasar light curves from SDSS only (10 yr) to combined SDSS–PS1 (15 yr).

§ DATA

We focus on the data pertaining to a 290 deg$^{2}$ region of the southern sky known as S82, repeatedly observed by the SDSS between 1998 and 2008. Originally aimed at supernova discovery, objects in this area were reobserved 60 times, on average (see MacLeod et al., 2012, Section 2, for overview and Annis et al., 2014 for details). Availability of well-calibrated [Ivezić et al., 2007], long-baseline light curves spurred variability research [Sesar et al., 2007]. The DR9 catalog [Schneider et al., 2008] contains 9258 spectroscopically confirmed quasars within S82. Using a $0.5''$ matching radius, we find corresponding data for 9248 quasars from  the PS1 DR2 [Chambers et al., 2016, Flewelling, 2018, Flewelling et al., 2020], 7737 from the Catalina Real-Time Transient Survey (CRTS; Drake et al., 2009),  6455 from the Palomar Transient Factory (PTF; Rau et al., 2009), and 8001 from the Zwicky Transient Facility (ZTF) DR1 [Bellm et al., 2019, Masci et al., 2019]. In this section, we first consider the possibility of utilizing data from all surveys and the challenges involved in combining the data from varying filter sets. We conclude that, with its shallower photometry (Figure <ref>) and baseline overlapping with SDSS–PS1, the added data from the CRTS would be compromised by photometric uncertainty (Figure <ref>), and the necessity of converting from a broad CRTS $V$-band filter (white light) to SDSS $r$. Similarly, the ZTF and PTF  were deemed too shallow for optimal data combination (median 0.1 mag uncertainty, as opposed to 0.02 mag for SDSS and 0.03 mag for PS1; see Figure <ref>). Finally, utilizing photometric offsets, although considered and calculated (Figure <ref>), would complicate the investigation into quasar variability by adding an additional  layer of uncertainty, and for this reason, we use only SDSS and PS1 $r$-band data, as these bandpasses are sufficiently similar that no offset is required (see Figure <ref>).

Figure <ref> illustrates the improvement in baseline coverage when combining various surveys. The length of each thick dashed line corresponds to the duration of each survey (survey baseline), and the size of each circle corresponds to the area covered by each survey. The vertical location of each dashed line marks the $5\sigma$ depth in the $r$-band (or equivalent). The Rubin Observatory Legacy Survey of Space and Time (LSST) stands out in that it will provide the best extension of SDSS baseline and depth.

Illustration of survey baseline, sky area covered, and depth. The length of each thick dashed line corresponds to the extent of the real or simulated light curves for S82 quasars for each survey. This includes SDSS DR7, CRTS DR2, PS1 DR2, PTF DR2, ZTF DR1, and, for LSST, the full 10 yr survey. The vertical location of each line corresponds to the $5\sigma$ limiting magnitude (SDSS $r$, PS1 $r$, PTF $R$, ZTF $r$, LSST $r$, CRTS $V$). The size of each circle represents the total survey area (for SDSS, up to DR15).  Note how PS1 and PTF extend the baseline of SDSS by approximately $50\%$ and that inclusion of LSST roughly triples the SDSS baseline. For reference, the area covered by LSST is $20,000$ deg${^2}$.

Combining data from different photometric standards requires applying color transformation or photometric offsets. We first seek to combine PS1 $gri$,  PTF $gR$, and CRTS $V$ into a common SDSS $r$ band (best photometry). To this end, we calculate  color terms using the SDSS standard star catalog [Ivezić et al., 2007]. Focusing on a 100,000 randomly chosen stars, we find their CRTS, PS1, and PTF matches. [CRTS from B.Sesar, priv.comm., PS1 from MAST (<http://panstarrs.stsci.edu>), and PTF from IRSA PTF Object Catalog (<https://irsa.ipac.caltech.edu/>)]

The difference  between the target (SDSS) and source (e.g.,PS1) photometry can be written as a function of the mean SDSS $g-i$ color:

\begin{equation}
m_{\mathrm{PS1}} - m_{\mathrm{SDSS}} = f(g-i)
\end{equation}

Some authors (e.g., Li et al., 2018) allow the transformation to be a higher-order polynomial, but as Figure <ref> shows, quasars occupy a relatively narrow region of $g-i$ color space, and we find that the linear fit is sufficient. The derived linear coefficients for photometric transformations between SDSS $r$ and PS1 $gri$, PTF $gR$, and CRTS $V$ as a function of SDSS $g-i$ color are listed in Table <ref>. We illustrate the process, showing, in Figure <ref>, the SDSS–PS1 standard star data used to calculate the offsets. Note that the PS1 $r$ (middle panel) is very close to the SDSS $r$, within  $0.01$ mag (1% level) across the $g-i$ color range. We focus on the SDSS $r$–PS1 $r$ offset as a function of magnitude in Figure <ref>; the near-equivalence of bandpass coverage is valid at the 1% level up to $r < 20.5$.

In selecting the most beneficial datasets to complement SDSS $r$, we also consider the associated photometric uncertainties (aka `errors'). As shown in Figure <ref>,  PTF and CRTS are shallower than SDSS or PS1 ($\sim20.5$ mag vs $~22$ mag). Therefore, for faint objects like quasars (for the S82 sample, the population median is SDSS $r{\sim}20$ mag), PTF and CRTS have larger photometric uncertainties than SDSS or PS1. Indeed, as Figure <ref> shows, the distribution of median errors for PTF, CRTS, and ZTF quasar data (median $\sim0.1-0.15$ mag) is wider than the corresponding SDSS and PS1 data (median ${\sim}0.02-0.04$ mag). As simulations show (Sec. <ref>), although PTF and CRTS data do extend the SDSS baseline, their error properties decrease their utility in complementing the SDSS data set. After all, the SDSS baseline extension afforded with PTF and CRTS is comparable to that achieved with PS1 data alone (Figure <ref>).

Furthermore, to mitigate problems that could arise when applying photometric transformations (such as spurious variability due to incorrect offsets or color-dependent variability), we choose to combine SDSS $r$ with only PS1 $r$, since, as Figures <ref> and  <ref> show, SDSS $r$ and PS1 $r$ are sufficiently similar (at a 1% level up to 20.5 mag) that no photometric transformation is required.

Finally, we clean the combined SDSS $r$–PS1 $r$ quasar light curves using standard procedures of  $\sigma$-clipping in magnitude and error space and error-weighted day averaging to mitigate the impact of bad photometry and average out the intra-night variability (as in Charisi et al., 2016, Suberlak et al., 2017). Of 9248 SDSS–PS1 quasars, 8516 have the PS1 $r$ data with 662,092 epochs.  We remove points that have errors departing from the median SDSS(PS1) light curve segment by more than $7 \sigma$, and we visually inspect all photometry with magnitudes departing by more than  $7\sigma$ from the median magnitude. Of 585 flagged light curves, 253 required removal of individual epochs containing bad photometry. To avoid unphysically small errors,  we add in quadrature $0.01$ mag if the combined nightly error is $<0.02$ mag. In the final sample, there are 580,321 epochs.

Regions of color-color (top left, top right, bottom left) and color-magnitude (bottom right)  space occupied by SDSS S82 quasars (colors) and stars (contours). We use quasar median photometry from Schneider et al., 2010 and the standard star catalog of Ivezić et al., 2007, showing a random subset of 10,000 stars. As seen in the bottom left panel, quasars occupy a particular range of SDSS $g-i$ color. Therefore, in fitting the linear color transformations, we limit the color range to $-0.35<(g-i)<0.75$ (vertical dashed lines in Figure <ref>). Quasars also overlap other variable sources (e.g., RR Lyrae) not shown here [Sesar et al., 2007].

The SDSS–PS1 offsets, derived with the SDSS standard stars [Ivezić et al., 2007]. From a randomly chosen subset of 100,000 SDSS stars, 95,000 have PS1 DR2 data. To minimize scatter due to larger errors, we select 40,000 stars with  $r < 19$ mag. Vertical dashed lines mark the region in the SDSS color space occupied by quasars, indicated by the contours enclosing  90%, 70%, and 30% quasar data (also see Figure <ref>). Stars with $-0.2<g-i<0.8$ were used to fit the stellar locus with a first-order polynomial, marked by the dotted-dashed line. The best-fit slopes of 0.619, -0.04, and -0.283 for PS1 $g$, $r$, and $i$, respectively, are listed as $B_{1}$ in Table <ref>.

Color terms (photometric offsets) between CRTS, PTF, and PS1 Passbands and SDSS Using the SDSS Mean $g-i$ Color to Spread the Stellar Locus.

Band (x)

CRTS V     -0.0464      -0.0128

PTF g      -0.0294       0.6404

PTF R      0.0058       -0.1019

PS1 g      0.0174        0.6194

PS1 r      0.0065       -0.0044

PS1 i      0.0260       -0.2830

The SDSS $r$ synthetic magnitude, $r_{s}$, can be found as $r_{s} = x-B_{0}-B_{1}(g-i)$. This linear trend is illustrated in Figure <ref>, where we plot $(x-r_{\mathrm{SDSS}})$ as a function of $(g-i)_{\mathrm{SDSS}}$ for $x=g_{\mathrm{P1}}, r_{\mathrm{P1}},i_{\mathrm{P1}}$. To derive the color terms, we used a subset of 100,000 stars randomly chosen from the SDSS standard star catalog [Ivezić et al., 2007]. To minimize scatter, we selected bright stars with $r<19$ mag.

The PS1 $r$ versus SDSS $r$ as a function of SDSS $r$ for 100,000 randomly selected standard stars from the Ivezić et al., 2007 catalog. Almost 95% of SDSS stars have PS1 DR2 photometry. The filled circles represent the median offset, a slight slope at the 1% (0.01 mag) level, up to $r<20.5$ mag.

Distribution of median photometric uncertainties (`errors') in $r$-band real light curves. The PTF and ZTF surveys are shallower than SDSS and PS1 (see Figure <ref>), and thus for faint objects, such as quasars, they have much larger median errors. The CRTS errors (not shown) are, on average, 50% larger than PTF.

In conclusion, we use SDSS S82 data spanning 10 yr in the observed frame extended by the PS1 DR2 data that add, on average, 5 yr of data. Although we considered utilizing ZTF, PTF, and CRTS data, their larger photometric uncertainties offset any gain due to the extended baseline. Furthermore, only the SDSS and PS1 $r$ bands are sufficiently similar (at a ${<} 0.01$ mag level; Figure <ref>) that no photometric offsets are required, which would unnecessarily complicate merging light curves from different surveys.

§ SIMULATIONS : LESSONS LEARNED

We simulate the theoretical improvement of the DRW parameter retrieval in extended light curves. We generate long and well-sampled  `master' light curves, all with input $\tau = 575 $ days, SF$_{\infty} = 0.2$ mag (the median of S82 quasar distribution in M10), with zero mean.  We subsample at real observed epochs for SDSS and PS1 and at predicted cadences for ZTF and LSST  (see Figure <ref>). To each simulated light curve, we add a magnitude offset corresponding to the mean of the combined SDSS–PS1 light curve. That way, the magnitude distribution of simulated light curves is similar to that of the observed SDSS–PS1 data. For the LSST 10 yr segment (finishing in 2031), we assumed 50 randomly distributed  epochs per year, with the following error model:

\begin{eqnarray}
\label{eq:errorModel}
\sigma_{LSST}(m)^{2} &=& \sigma_{sys}^{2} + \sigma_{rand}^{2} \,\, \mathrm{(mag)}^{2} \\
\sigma_{rand}^{2} &=& (0.04-\gamma)x + \gamma x^{2} \nonumber \\
x &=& 10^{0.4(m-m_{5})} \nonumber
\end{eqnarray}

with  $\sigma_{sys} = 0.005$, $\gamma=0.039$, and $m_{5} = 24.7$ (see Ivezić et al., 2019, Sec. 3.2). For the ZTF 1 yr segment (spring 2019 ZTF DR1, including the data from 2018), we assumed 120 observations (every 3 nights) in $g_{\mathrm{ZTF}}$ and $r_{\mathrm{ZTF}}$, deriving the magnitude-dependent error model by plotting the best mag rms as a function of the best median magnitude for ZTF matches to S82 standard stars in Figure <ref>. We find that the LSST error model (Equation <ref>) with $\gamma = 0.05$, $\sigma_{sys} = 0.005 $, and $m_{5} = 20.8$ adequately describes the ZTF photometric uncertainty.

Best mag rms plotted as a function of magnitude for ZTF nonvariable stars with over 100 observations. We overplot the adopted error model, with $\gamma = 0.05$, $\sigma_{sys} = 0.005 $, and $m_{5} = 20.8$ (see Equation <ref>). The properties of the ZTF photometric uncertainties are largely similar to the PTF uncertainties.

Simulated well-sampled underlying DRW process, one of `master' light curves ($\tau=575$ days, SF$_{\infty} = 0.2$ mag, 4 points day$^{-1}$) shown with  small black dots. To simulate observations, the cadence is degraded (subsampled) to match the ground-based cadence corresponding to real quasar data from SDSS (red), PS1 (green) segments, and simulated LSST (blue) epochs (here we use SDSS–PS1 epochs for quasar dbID=3537034). The orange `error snake' is an envelope marking the standard deviation of the fit to the data using a GP with a DRW kernel (Sec. <ref>).

To mirror observational conditions, we add a Gaussian noise to the true underlying DRW signal with variance defined by photometric uncertainties for corresponding surveys. Figure <ref> illustrates the simulated `master' light curve (black dots, 4 day$^{-1}$), subsampled at SDSS (red), PS1 (green), and LSST (blue) cadence. While PS1 provides a 50% improvement of the SDSS baseline, LSST will nearly triple it. Figure <ref> shows how the simulated distribution of DRW parameters $\sigma$, $\tau$ changes as the SDSS quasar light curves are extended with PS1, ZTF, and LSST data. We quantify the improvement in recovery of DRW parameters with an rms error calculated for $\sigma_{fit}/\sigma_{in}$ and $\tau_{fit} / \tau_{in}$. For ratio $r$,  rms$(r)^{2} = \mathrm{bias}(r)^{2} + \sigma_{G}(r)^{2}$. For $r=\sigma_{fit}/\sigma_{in}$, the rms changes from 0.322 (SDSS only), to 0.273  (SDSS–PS1), to 0.305 (SDSS–PS1–ZTF), to 0.182 (SDSS–PS1–ZTF–LSST). In the future (after more data have been assembled and recalibrated), ZTF will help, but not as dramatically as LSST (which provides a factor of 1.8 improvement in rms errors in $\sigma_{fit}/\sigma_{in}$ and a factor of 1.5 for $\tau_{fit} / \tau_{in}$). Note that, due to larger errors (Figure <ref>), including ZTF causes a widening of the recovered $\tau$ distribution - the rms increases from 0.715 with SDSS to 0.890 with SDSS–PS1–ZTF (see Figure <ref>). Using PS1 data with their excellent  deep photometry (as compared to ZTF or PTF) is the best improvement over existing SDSS results (factor of 1.2 decrease of rms errors for $\sigma$, and a smaller increase for rms errors for $\tau$; factor of 1.25 comparing SDSS–PS1–ZTF, versus factor of 1.07 for SDSS–PS1). For this reason, we use only the SDSS–PS1 portion of quasar light curves as the best trade-off between adding more baseline vs introducing more uncertainty with noisy data.

Ratio of DRW parameters fitted with celerite: $\tau$ and $\sigma$, to the input $\tau_{in} = 575 $d, $\sigma_{in} = 0.2 / \sqrt{2} {\sim} 0.14$  (SF$_{\infty}=0.2$ mag). We simulated 9258 `master' light curves, subsampled at real SDSS or PS1 $r$-band cadence and photometric uncertainties, and simulated ZTF and LSST cadence. To simulate observing conditions, the underlying DRW signal was convolved with a Gaussian noise corresponding to epochal errors. For each light curve, we start with the SDSS segment only, and as we add more segments (PS1, ZTF, LSST), we refit for DRW model parameters with celerite. Thus,	each distribution corresponds to a different segment of the simulated  combined SDSS–PS1–ZTF–LSST light curves. Extending the baseline shifts the distribution of the recovered DRW parameters toward the unbiased regime; the vertical dashed line marks the input matching the output. This corresponds to the top right (well-constrained) portion of Figure <ref>.

In summary, to gauge the improvement in the recovery of DRW parameters due to extending quasar light curves, we simulate $8516$ well-sampled light curves, subsampled at real cadence for SDSS and PS1 segments, and predicted cadences for ZTF and LSST segments. The mean brightness and the photometric uncertainties for each simulated light curve  mimic those in the real SDSS–PS1 data. Importantly, since the true underlying distribution of $\tau$, $\sigma$ is not known, we assume that the DRW parameters for all light curves correspond to the mean of the results of the M10 SDSS-based study: the same input timescale $\tau=575$ days, and amplitude $\sigma_{in}=0.14$ mag.  We find that combining SDSS and PS1 provides a factor of 1.2 decrease in rms errors for $\sigma$, balanced by a modest increase in rms errors for $\tau$.

§ RESULTS: VARIABILITY PARAMETERS FOR S82 QUASARS

We extend the S82 quasar light curves by combining the SDSS $r$-band data with  the PS1 $r$-band data without any photometric offsets. For each quasar, we fit the SDSS and SDSS–PS1 segments  with the DRW model. This yields two sets of DRW parameters per quasar: $(\tau_{SDSS}, \sigma_{SDSS})$, and $(\tau_{SDSS-PS1},\sigma_{SDSS-PS1})$. Because variability is inherent to the quasar, for the remaining analysis, we shift all fitted timescales to the quasar rest frame and implicitly assume that the DRW timescales are considered in the rest frame: $\tau_{\mathrm{RF}} = \tau_{\mathrm{OBS}} / (1+z)$.

Comparison of distributions of the rest-frame variability timescale $\tau_{RF}$ against the  asymptotic variability amplitude SF$_{\infty}$, for M10 SDSS $r$-band (red solid)  and celerite fits using  SDSS (blue dashed) or SDSS–PS1 (green dotted-dashed) segments of combined S82 quasar light curves. Contours show the 1$\sigma$ and 2$\sigma$ levels (enclosing $68.3\%$ and $95.5\%$ of the data). The timescales and SF$_{\infty}$ from M10 and this work overlap, as we recover the same underlying distributions.

In this section, we first correct the fitted $\tau$, $\sigma$ for wavelength dependence. Then we show consistency with the M10 results and consider the trends between DRW parameters and physical quasar properties: black hole mass $M_{\mathrm{BH}}$, absolute $i$-band  magnitude $M_{i}$, or redshift $z$.

§.§ Comparison to M10

The DRW parameters recovered with celerite are  broadly consistent with M10; Figure <ref> shows the rest-frame  $\tau$ and SF$_{\infty}$ distributions for our results for the SDSS segment (blue dashed contours),  SDSS–PS1 combined light curves (green dotted-dashed contours), and  M10 SDSS for $r$-band only (red solid contours). When using exactly the same data as M10 (SDSS), our results agree. The offset of 0.05 dex  between our results and those of M10 results for SDSS, seen in the left panel of Figure <ref>, can be attributed to data cleaning and software differences. The right panel of Figure <ref> shows the same distribution in terms of $K-\hat{\sigma}$ space, orthogonal to $\tau-\sigma$, where $\hat{\sigma} = \sigma\sqrt{2 / \tau}$, and $K = \tau \sqrt{\sigma} 2^{1/4} $.

Comparison of celerite fits using only the  SDSS $r$-band segments of S82 quasars ($\sigma_{\mathrm{SDSS}}, \tau_{\mathrm{SDSS}}$) against M10 results for the SDSS $r$ band ($\sigma_{M10}, \tau_{M10}$) object by object. Note that for each ratio, the median-based bias (bottom left corner) is calculated before taking the logarithm. The small offset ($<0.05 $ dex) can be attributed to data cleaning and software differences. See Figure <ref> for a comparison of rest-frame $\tau$ and SF$_{\infty}$ distributions. This is similar to Figure 3 in M10, except we plot only the $r$-band SDSS results. The right panel shows the comparison in an orthogonal $K-\hat{\sigma}$  space, where $K$ is the direction along the diagonal in the left panel and $\hat{\sigma}$ is perpendicular to the diagonal. For this reason, the right panel has a 10 times smaller scatter along $\hat{\sigma}$ (0.037) than $K$ (0.358).

Ratios of fitted DRW parameters  ($\tau$, $\sigma$), comparing the value of the parameter recovered using the combined light curve length (SDSS–PS1) to the shorter, SDSS-only light curve (labeled S20 in Figure <ref>). The left panel shows the results for S82 quasars using real data, whereas the right panel shows the simulated quasars with realistic cadence, with  $\tau_{in}=575$ days and  SF$_{\infty}=0.2$. The general trend when using the real data (despite having a range of underlying timescales and amplitudes) is similar to that when using simulated data: the diagonal scatter is along the lines of constant $\hat{\sigma}$, and there is much less scatter in the perpendicular direction of $K$ (see Figure <ref>). There is no major change in the shape of the distribution as a function of mean quasar magnitude. The red rectangle marks the outliers belonging to the tail end of the population, with $\log{(\tau_{\mathrm{SDSS-PS1}} /  \tau_{\mathrm{SDSS}})} > 1$   and   $\log{(\sigma_{\mathrm{SDSS-PS1}} / \sigma_{\mathrm{SDSS}})}   > 0.4  )$, discussed in Sec. <ref>.

Same as Figure <ref>, but in $K-\hat{\sigma}$ space, which is orthogonal to the $\tau-\sigma$ space, since $K = \tau \sqrt{\mathrm{SF}_{\infty}} = \tau \sqrt{\sigma} 2^{1/4}$ and $\hat{\sigma} = \mathrm{SF}_{\infty} / \sqrt{\tau} = \sigma \sqrt{2/\tau}$.

§.§ Outliers: possible CLQSO candidates

Figures <ref> and  <ref> show the change in recovered DRW parameters between SDSS and combined SDSS–PS1 light curves. In Figure <ref> the distribution of $f_{\sigma} \equiv \log_{10}{(\sigma_{\mathrm{SDSS-PS1}} / \sigma_{\mathrm{SDSS}})}$ and $f_{\tau} \equiv \log_{10}{(\tau_{\mathrm{SDSS-PS1}} / \tau_{\mathrm{SDSS}} )}$ for real light curves (left) matches the predicted distribution for simulated light curves (right). Analogously, on Figure <ref> these ratios converted to a $K-\hat{\sigma}$ space also show close correspondence between the observed (left) and simulated (right) light curves. Studies show that about 0.1% of quasars will exhibit large variability (in excess of 0.5 mag rms; see Figure 18 in MacLeod et al., 2012). Visual inspection of light curves in the top right region of the left panel of Figure <ref> (marked by the red rectangle) reveals large changes in brightness, similar to those seen in changing-look quasars (CLQSO,see Elitzur et al., 2014, LaMassa et al., 2015, Schawinski et al., 2015, Guo et al., 2016, Ruan et al., 2016, Runnoe et al., 2016, Blanchard et al., 2017, Gezari et al., 2017, Stern et al., 2017, Sheng et al., 2017, Lawrence, 2018, Ross et al., 2018, Stern et al., 2018, Yang et al., 2018, Frederick et al., 2019, MacLeod et al., 2019, Ruan et al., 2019, Trakhtenbrot et al., 2019, Shen et al., 2019, Sheng et al., 2020). The rectangle in Figure <ref> marks the region with  $f_{\tau} > 1$, $f_{\sigma} > 0.4$. There are 28 objects in the simulated sample and 48 objects in the observed sample with that property, which reflects the fact that we do expect quasars to exhibit changes in their brightness as a function of time, with larger variance possible over longer timescales. We further select 40 of these that are brighter than  $20.5$ median PS1 $r$-band magnitude, since PS1 observations correspond to more recent epochs and provide a better indication of the possibility of follow-up. The light curves and properties of these CLQSO candidates are further discussed in Appendix <ref>.

Such large differences in timescale and amplitude of variability can also be inferred directly from the light curves. Consider the difference in magnitude and scatter between the SDSS portion of the light curve (spanning approximately 10 yr between 1998 and 2008) and the PS1 portion (spanning ${\sim}5$ yr between 2009 and 2014; see Figure <ref>). We measure the median magnitudes offset as $\Delta(\mathrm{median}) =  \mathrm{median}(SDSS) - \mathrm{median}(PS1)$, and the scatter difference as $\Delta(\sigma_{G}) = \sigma_{G}(SDSS)-\sigma_{G}(PS1)$. The resulting distributions of  $\Delta(\mathrm{median}) $ and $\Delta(\sigma_{G})$ for S82 quasars are shown in Figure <ref>. Indeed, when plotting $\Delta(\mathrm{median})$ as a function of $f_{\tau}$ and $f_{\sigma}$ there is a gradient indicating that the  CLQSO candidates, outliers in $(f_{\tau}, f_{\sigma})$ space, are also outliers in $\Delta(\mathrm{median})$– $\Delta(\sigma_{G})$ space. Thus, the by-product of extending light curves to recalculate the DRW parameters with increased fidelity is an independent method to discover the CLQSO.

Differences between SDSS and PS1 segments of combined quasar $r$-band light curves. First, we show the difference between the median SDSS and median PS1 portion, plotted as a histogram (top left panel) and cumulative distribution function (bottom left panel). Then, we show difference between the $\sigma_{G}$ calculated for each portion of the light curve ($\sigma_{G}$ is a robust estimate of the standard deviation and is related to the difference between the 75th and 25th percentiles: $\sigma_{G} = 0.7413(Q_{75} - Q_{25})$). The outliers in the median offset space are also outliers in the DRW parameter space (e.g., objects with $\log{(\tau_{\mathrm{SDSS-PS1}} /  \tau_{\mathrm{SDSS}})} > 1$ and $\log{(\sigma_{\mathrm{SDSS-PS1}} / \sigma_{\mathrm{SDSS}})}   > 0.4  )$, and $r > 20.5$ have $\Delta(\mathrm{median}) > 0.1$).

§.§ Rest-frame Wavelength Correction

Objects at cosmological distances are embedded in the Hubble flow due to the expansion of the universe [Riess et al., 2019]. Therefore, light observed from a distant quasar would have been emitted at shorter wavelengths in the quasar's rest frame, $\lambda_{RF} = \lambda_{obs} / (1+z)$, where $z$ is the cosmological redshift. Quasars at different redshifts probe different regions of rest-frame spectra (see Figure 7 in Shen et al., 2019). Thus, before correlating the DRW parameters with quasar properties, we correct $\sigma, \tau$ for the $\lambda_{RF}$ dependence studied by M10 with SDSS $ugriz$ light curves. We plot the DRW parameters in Figure <ref>: SF$_{\infty}$ and $\tau$ as a function of $\lambda_{RF}$. A solid line marks the M10 best-fit power-law trend,

\begin{equation}
\label{eq:lambda}
f \propto \left( \frac{\lambda_{RF}}{4000 \mbox{\AA}} \right)^{B}
\end{equation}

with  $B=-0.479$ and $0.17$ for SF$_{\infty}$ and $\tau$, respectively.

Rest-frame timescale $\tau$ (top panel) and asymptotic SF: SF$_{\infty}$ (bottom panel) as a function of rest-frame wavelength $\lambda_{RF}$. The background contours show the 30% and 70% levels for  M10 SDSS $ urz $  data, and the foreground contours  denote our results using  SDSS (red) and SDSS–PS1 (orange) segments. The red line indicates the best-fit power law to M10 data, with $B=0.17$ and $-0.479$ for $\tau_{RF}$, and SF$_{\infty}$, respectively. As M10 showed, this means that the timescale is almost independent from the bandpass, while the variability amplitude decreases toward redder rest-frame wavelengths.  We take the center of each bandpass to approximate the  observed wavelength; that is, for SDSS $urz$ bandpasses,  $\lambda_{obs} = 3520$, $6250$, and $9110$ $\mbox{\AA}$, respectively, and given the redshift of each quasar, we find $\lambda_{RF}=\lambda_{obs} / (1+z)$.

§.§ Trends with Black Hole Mass, Absolute Luminosity

In the era of large synoptic surveys such as ZTF or LSST, the large increase in the number of discovered quasars means that due to limited observational resources, we can afford a spectroscopic follow-up for only a few percent of AGNs with optical time series [Ivezić et al., 2019]. Therefore, a relationship between the quasar variability parameters ($\tau, \sigma$) and physical properties $M_{BH}$, $M_{i}$ could provide an estimate of the latter for millions of quasars. We inspect correlations between $\tau, \sigma$  and $M_{BH}$, $M_{i}$, using the [Shen et al., 2011] catalog, based on single-epoch SDSS spectra. Here $M_{i}$ is $k$-corrected to  $z=2$, corresponding to the peak of quasar activity [Richards et al., 2006]. For details, see Appendix <ref>.

In Figure <ref>, we examine the distribution of $M_{BH}$, $M_{i}$, as a function of $z$ for S82 quasars. The upward gradient in the top two panels reflects the selection effect that higher-redshift quasars have to be brighter to be included in the magnitude-limited sample (luminosity–redshift degeneracy; see Section 5, Figure 12 in M10, and Dong et al., 2018). Higher-redshift quasars are also more active  and have higher black hole masses due to cosmological downsizing (see McLure & Dunlop, 2004, Babić et al., 2007, Labita et al., 2009). The distribution in the bottom left panel of Figure <ref> is peaked at $z=2,$ which corresponds to the peak of quasar activity.

Distribution of quasars as a function of  redshift, observed $i$-band magnitude, absolute $i$-band magnitude ($k$-corrected to z=2), and virial black hole mass. All data are from [Shen et al., 2011].

Figure <ref> shows the DRW parameters for S82 quasars, $\tau$ and SF$_{\infty}$, plotted as a function of quasar physical properties $M_{BH}$, $M_{i}$, and $z$. The left  panels in Figure <ref> contain a gradient of SF$_{\infty}$  with
$M_{i}$,$z$: brighter quasars have lower variability amplitude, largely independent of black hole mass.

Absolute $i$-band magnitude $M_{i}$  ($k$-corrected to $z=2$) as a function of the virial black hole mass $M_{BH}$ (bottom panels) and redshift $z$ (top panels), colored by the  long-term variability amplitude SF$_{\infty}$ (left panels) or characteristic timescale $\tau$ (right panels). Here $M_{i}$  is a proxy for quasar bolometric luminosity, and the DRW variability parameters are for combined SDSS–PS1 $r$-band data.

We investigate these relations in more detail by fitting $f$ ( $\tau$ or SF$_{\infty}$) as a power-law function of $M_{BH}$, $M_{i}$, and $z$,

\begin{eqnarray}
\label{eq:powlawmodel}
\log_{10}{f} = &A& + B \log_{10}\left( \lambda_{RF} / 4000 \mbox{\AA} \right) + C (M_{i} + 23) \nonumber \\
&+& D \log_{10}{\left( M_{BH} / 10^{9} M_{\odot}  \right)}
\end{eqnarray}

using a Bayesian linear regression method that incorporates measurement uncertainties in all latent variables [Kelly, 2007]. This ansatz is identical to that used by M10, and very similar to the relation used by [Scaringi et al., 2015], since, for black holes, their Equation (1) becomes $\log{t_{b}} = A' \log{M_{BH}} + B' \log{L_{Bol}}+C'$, where $t_{b}$ is the PSD break timescale.

We compare the change in retrieved fit coefficients caused by adding PS1 data to SDSS against the M10 SDSS-only study. Note that M10 fitted a DRW model treating each of the five SDSS bands as a separate light curve, resulting in over 30,000 values of $\tau, \mathrm{SF}_{\infty}$ for 9000 S82 quasars. Grouping the fitted quasar parameters by band, they were correlated to the quasar physical parameters with Equation <ref>. Figure <ref> shows the posterior samples for fitting  Equation <ref> to $f=\mathrm{SF}_{\infty}$ for quasar data separately for each SDSS bandpass.  Each band yields a slightly different fit coefficient. As the fit result, M10 reported the band mean (red vertical dashed line). Since this study uses only $r$-band data, we compare the fit coefficients to M10 SDSS $r$ data (green vertical solid line in Figure <ref>). We show the results of fitting Equation <ref> to new SDSS and SDSS–PS1 parameters in Figures <ref> and  <ref>. First, with $f=\tau$ in Equation <ref> (Figure <ref>), the SDSS–PS1 data confirm M10 for luminosity dependence (the posterior Markov Chain Monte Carlo (MCMC) samples overlap), but the dependence of $\tau$  on $M_{BH}$  is marginally weaker (by 0.007 dex). Second, in Figure <ref> with  $f=\mathrm{SF}_{\infty}$, $\mathrm{SF}_{\infty}$ has a slightly weaker dependence on $M_{BH}$ (by 0.04 dex compared to M10). The difference between the celerite SDSS-only results and M10 can be attributed to data cleaning that was not performed by M10, and software differences. Each distribution from Figures <ref> and  <ref> is summarized in Table <ref>, with the uncertainty in the A,C, and D fit coefficients estimated from  the standard deviation of the posterior samples.

We also searched for a signal of Mg ii variability (see Cackett et al., 2015 for a review). Like H$\alpha$  and  H$\beta$ , Mg ii is a permitted low-ionization line [Yang et al., 2020], but on average, it is being emitted by gas further away from the ionizing source than the  H$\beta$, possibly at the edge of the broad-line region (BLR; Guo et al., 2020). [Ivezić et al., 2004] and [MacLeod et al., 2012] studied the SDSS–POSS sample of quasars and reported the detection of a decrement in the difference between the data and the best-fit model  (residuals) around 2800 $\mbox{\AA}$ when plotting the residuals as a function of $\Delta t$ and $\lambda_{RF}$. We investigate the residuals for $f=\mathrm{SF}_{\infty}$ in Equation <ref>. We find that, using the SDSS data, the decrement in the median($\tau_{\mathrm{RF}}$) around 2800 $\mbox{\AA}$ is visible at ${\sim}3 \sigma$  relative to the smooth model, but by adding the PS1 data, the significance rises to ${\sim}5 \sigma$. We do not see much difference with regard to whether using the subset of 6371 quasars for which M10 had reliable results (listed in Table <ref>) versus  the full set of 8516 quasars fitted with celerite for DRW parameters. The effect is interesting but does not produce a very significant signal; see Appendix <ref> for more details.

Posterior MCMC draws for fitting Equation <ref> with M10 variability amplitude SF$_{\infty}$ against $M_{BH}$, $M_{i}$, and $z$ [Shen et al., 2011]. Since M10 treated the near-simultaneous SDSS $ugriz$ data for 9258 quasars independently for each band, this resulted in  DRW fit parameters for 7014 $u$-,  7408 $g$-,  6871 $r$-, 6814 $i$-, and 5111 $z$-band SDSS quasar light curves that fulfilled the M10 quality of DRW fit selection criteria. The M10 values for SF$_{\infty}$ are corrected to 4000 $\mbox{\AA}$ using Equation <ref>, with the power-law coefficient $B=-0.479$. Each distribution corresponds to a different SDSS band. We compare the results of fitting the SDSS–PS1 $r$ band directly against the M10 results for the SDSS $r$ band (green solid). Note that Table 1 in M10 reported band-averaged values for the A, C, and D coefficients (red vertical dashed  line), while we cite in Table <ref> the mean for the $r$ band (green vertical solid line).

Comparison of best-fit coefficients for Equation <ref> using M10 Results and This Work (S20).

$f$     Source     $A$(offset)     $B(\lambda_{RF})$     $C (M_{i})$     $D (M_{\mathrm{BH}})$
$\tau$     M10, SDSS     $2.5\pm0.027$     $0.17\pm0.02$     $0.03\pm0.009$     $0.178\pm0.027$

S20, SDSS     $2.515\pm0.019$     $0.17\pm0.02$     $0.042\pm0.007$     $0.127\pm0.019$

S20, SDSS–PS1     $2.597\pm0.02$     $0.17\pm0.02$     $0.035\pm0.007$     $0.141\pm0.02$

SF$_{\infty}$     M10, SDSS     $-0.486\pm0.012$     $-0.479\pm0.005$     $0.119\pm0.004$     $0.121\pm0.012$

S20, SDSS     $-0.543\pm0.009$     $-0.479\pm0.005$     $0.125\pm0.003$     $0.104\pm0.008$

S20, SDSS–PS1     $-0.476\pm0.008$     $-0.479\pm0.005$     $0.118\pm0.003$     $0.118\pm0.008$

Here $B$ is fixed to $0.17$ or $-0.479$ from fitting a power law between $\lambda_{RF}$ and $\tau$, SF$_{\infty}$ (see Figure <ref>). Of 8516 quasars with SDSS–PS1 data, for consistency we use here an unbiased subset of 6371 quasars for which M10 had reliable results. For $f=\tau$, $C$ is almost the same between M10 and this work for SDSS–PS1 (rows 1 and 3). However, $D$ based on SDSS–PS1 data is larger than M10 by 0.01 dex (row 3). For $f = \mathrm{SF}_{\infty}$, SDSS–PS1 based C is within 0.01 dex from M10 (rows 4,6), and $D$ based on SDSS–PS1 data is almost identical to M10. When using celerite $\tau$ ($\sigma$) results for all 8516 quasars, the luminosity dependence is unchanged to within 0.01 dex, and the dependence on the black hole mass is stronger by 0.05 dex (0.02 dex), respectively. As shown on Figure <ref>, there is a small  offset between $\log_{10}{(\tau_{\mathrm{S20,SDSS}}/ \tau_{\mathrm{M10,SDSS}})}$  and $\log_{10}{(\sigma_{\mathrm{S20,SDSS}}/ \sigma_{\mathrm{M10,SDSS}})}$, attributed to data cleaning procedures and software differences, which contributes to a shift between $C$ and $D$ parameters for $\tau$ and SF$_{\infty}$ between M10 and S20.

Distribution of MCMC posterior draws fitting Equation <ref> for the characteristic timescale ($f=\tau$) based on SDSS $r$-band results (dashed blue line) and new SDSS–PS1 combined $r$-band results (solid green line). These are considered simultaneously as a function of quasar absolute magnitude $M_{i}$ (left panel) and black hole mass $M_{BH}$ (right hand side panel). Of 9258 spectroscopically confirmed quasars in S82, we employed 8516 that had PS1 matches, of which 6371  fulfill the M10 selection criteria (see M10, Sec 2.2). The results from the SDSS–PS1 light curves are consistent with M10 for the SDSS $r$ band.

Same as Figure <ref>, but fitting the DRW asymptotic amplitude ($f = \mathrm{SF}_{\infty}$ in Equation <ref>)  as a function of absolute magnitude $M_{i}$, and black hole mass $M_{BH}$. New data from PS1 are consistent with the earlier results of M10 on luminosity dependence but support a slightly weaker dependence of $\mathrm{SF}_{\infty}$ on $M_{BH}$ (by 0.06 dex).

§.§ Comparison to other studies: Eddington ratio

The Eddington ratio $(f_{Edd} {=} L_{Bol}/L_{\mathrm{Edd}})$ encodes accretion strength: the proximity of quasar bolometric luminosity to the theoretical Eddington limit, where $L_{\mathrm{Edd}} {=} 1.26 {\times} 10^{38} (M_{\mathrm{BH}} / M_{\odot})$ erg s$^{-1}$ [Shen et al., 2011]. Since $\tau$ and SF$_{\infty}$ depend on $M_{i}$ and $M_{BH}$, we investigate the possibility of the Eddington ratio being the driver of these observed trends. In Figure <ref>, we show $f_{Edd} $ as a function of $M_{i}$, $M_{BH}$, and SF$_{\infty}$. The first two panels depict  $f_{Edd} $ and SF$_{\infty}$ binned as a function of $M_{i}$ and $M_{BH}$. The third panel shows the quasar counts, and the fourth panel shows the bin means (black dots). The means are further binned along  $f_{Edd} $ (as in M10). Combined SDSS–PS1 data support SF$_{\infty}$  being inversely related to  $f_{Edd} $, with a power-law slope of $-0.207 \pm 0.03$, consistent with  $-0.23 \pm 0.03$ reported by M10. Observations are generally consistent with basic predictions from Table <ref>: $A$ increases with increasing  $M_{\mathrm{BH}}$, and  $f_{Edd} $ decreases as $L_{Bol}$ increases, while $\tau$ increases with  $L_{Bol}$. No model from  Table <ref> is rejected.

Absolute $i$-band magnitude $M_{i}$ as a function of black hole mass $M_{BH}$, where color encodes the Eddington ratio $f_{Edd} = L/L_{Edd}$ (first panel), variability amplitude SF$_{\infty}$  (second panel), or quasar count (third panel).
We only plot bins with more than five quasars. The fourth panel shows the median SF$_{\infty}$ as a function of the median $f_{Edd}$ (green crosses),  averaged in bins of $f_{Edd}$ (open circles). The bin width $w$ is found to ensure an equal number of points (crosses) per bin ($N$). The error bars are $\sigma_{y} = 1.25  \sigma_{G}(\mathrm{bin}) / N$, where $\sigma_{G}$ is the robust estimate  of the standard deviation ($\sigma_{G} = 0.7413 (Q_{75}-Q_{25})$). We assume the uncertainty along $f_{Edd}$ as $\sigma_{x} = w/\sqrt{12} $ (see Ivezić et al., 2014). The solid orange line is the best-fit slope,  $-0.207 \pm 0.031$, with the slope uncertainty estimated from the standard deviation of the posterior samples. The best-fit slope agrees with the M10 results ($-0.23 \pm 0.03$), plotted as a dashed magenta line.

Theoretical Predictions of Various Models Concerning Amplitude $A$ and/or Timescale $\tau$ of Variability, Eddington Ratio $f_{Edd}$, Bolometric Luminosity $L_{Bol}$, and Accretion Rate $\dot{m}$.

Theory     Prediction
Standard thin disk [Shakura & Sunyaev, 1973, Netzer, 2013]      $ A \nearrow$ as $f_{Edd} \searrow$ , $\tau \propto L_{Bol}^{1/2}$ (in Caplar et al., 2017)

Strongly inhomogeneous disk [Dexter & Agol, 2011] / local fluctuations [Cai et al., 2018]      $ A \nearrow$ as $f_{Edd} \searrow$

Variations in global accretion rate [Hawkins, 2007, Li & Cao, 2008, Zuo et al., 2012]      $ A \nearrow$ as $M_{BH} \nearrow$,  $ A \searrow$ as  $L_{Bol} \nearrow$

Eddington ratio reflecting AGN age [Martini & Schneider, 2003, Hopkins et al., 2005]     $A \nearrow$ as $f_{Edd} \searrow$

X-ray reprocessing [Kubota & Done, 2018]     $A \nearrow$ as $f_{Edd} \searrow $, and  $A \nearrow$ as $\dot{m} \searrow$

Observations cannot reject any model.

Comparison of Published Results on Correlating Observed Quasar Light-curve Properties: A Measure of Variability Amplitude A, and Timescale of Variability $\tau$ against the Physical Quasar Properties Black Hole mass $M_{\mathrm{BH}}$, Eddington Ratio $f_{\mathrm{Edd}}$, and Bolometric Luminosity $L_{Bol}$.

Measure of Amplitude / Timescale
A vs.                 $\tau$ vs.

$M_{\mathrm{BH}}$     $f_{\mathrm{Edd}}$      $L_{\mathrm{Bol}}$
$M_{\mathrm{BH}}$     $f_{\mathrm{Edd}}$      $L_{\mathrm{Bol}}$

[Wilhite et al., 2008]     SF (ensemble study)     +     -     -     x    x     x

[Kelly et al., 2009]     DRW: $\tau$-decorrelation timescale, $\sigma$       –     0     !–     +     x      $\sim$+

[MacLeod et al., 2010]     DRW: $\tau$,$\sigma$     +     !–     !–     +     x     $\sim$+

[Morganson et al., 2014]a     SF: A, $\gamma$      x     x     -     x     x     +

[Kozłowski et al., 2016]     DRW: $\tau$,$\sigma$      x     x     !–     +     x     $\sim$+

[Simm et al., 2016]b    EV, and PSD (break timescale)     0     –     –      0     0     0

[Caplar et al., 2017]     SF: $\tau$,$\sigma$     $\sim$     x     !–     x     x     +

[Rakshit & Stalin, 2017]c     DRW: $\tau$,$\sigma$     +     !-      $\sim$+     x     x    x

[Sun et al., 2018]d     SF: $\tau$,$\sigma$     x     x    !–     x     +     !+

[Li et al., 2018]e     SF: A, $\gamma$    $\sim$+      –      –     +     x     +

[Sánchez-Sáez et al., 2018]f     SF: A     0     –     x     x     x     x

This work     DRW: $\tau$-decorrelation timescale, $\sigma$     +     !–     +     +     x     $\sim$+

We list the correlations as positive (+), negative (–, i.e. anti-correlation), not found (0), or not studied (x). We further note if the correlation is strong (!), or weak ($\sim$).
aThe 105,783 quasars with SDSS–PS1 sparse data,  SF parametrized as $V(\Delta t | A,\gamma) = A (\Delta t/  1 \mathrm{yr})^{\gamma}$; $\gamma$ is the increase of SF with time lag $\Delta t$
bThe 90 X-ray-selected AGN, PS1 optical data. Variability characterized by normalized excess variance (EV). PSD characterized by break timescale.
cNarrow- and broad-line Seyfert 1 AGNs, $z<0.8$, CRTS data ($>50$ epochs, $5--9$ years baseline), DRW fitted with JAVELIN [Zu et al., 2011], no timescale correlations considered due to short baseline.
dThe 1004 SDSS quasars with  $0.5 \le z \le 0.89$
eThe 119,305 sparse quasar light curves from DECaLS and SDSS; SF as in [Morganson et al., 2014] above.
fThe 1348 QUEST-La Silla quasars, only amplitude of variability and excess variance.

In summary, following the standard procedure of shifting the DRW parameter $\tau$ to the quasar rest-frame, and correcting for wavelength dependence, we compare our results to M10. We show in  Figure <ref> that the new rest-frame distributions of $\tau$,$\sigma$ with SDSS, SDSS–PS1 are similar to those found by M10. There is a small ($<0.05$ dex) offset between the parameters fitted for SDSS light curves between this study and M10, that we attribute to data cleaning and software differences (Figure <ref>). When comparing SDSS and SDSS–PS1 combined data, we expect that there will be a long tail of objects that have a larger timescale and amplitude of variability with an extended baseline, corresponding to a true timescale that is longer than the initially probed time range with SDSS (Figure <ref>). Changing-look quasars belong to the tail end of quasar variability (Sec. <ref>), and we suggest 40 candidates in Appendix  <ref>, as well as a quick diagnostic plot of differences in mean magnitudes and rms scatter between PS1 and SDSS light-curve portions that could also be used to hunt for CLQSOs independently of DRW fitting (Figure <ref>). We explore the correlation between the DRW parameters obtained for SDSS and SDSS–PS1 light curves and quasar physical properties. We find that the general trends are in agreement with M10 at the $2\sigma$ level, with the largest deviations caused by different pipeline procedures, as outlined in Table <ref>.

§ DISCUSSION

§.§ Trends with Eddington ratio

Anticorrelation of  variability amplitude with Eddington ratio  has a variety of possible theoretical explanations. In the thin disk theory [Shakura & Sunyaev, 1973, Frank et al., 2002, Netzer, 2013], the radius of the emission region at a given wavelength increases with Eddington ratio and is inversely proportional to temperature [Rakshit & Stalin, 2017]. Thus, a hotter disk means that the emission observed in a given bandpass is emitted from a larger radius. From causality, a smaller region can be more variable than a larger one. Therefore, a  hotter disk would be less variable at a given wavelength than a colder one, and  the variability amplitude as studied in a particular bandpass (here SDSS $r$ band) would be anticorrelated with Eddington ratio [Fausnaugh et al., 2016, Edelson et al., 2015].

On the other hand, in the strongly inhomogeneous disk model, independent temperature fluctuations in $N$ zones drive the variability [Dexter & Agol, 2011]. In that framework, the inverse trend of variability amplitude against $L/L_{Edd}$  and $L_{Bol}$  can be understood qualitatively if more luminous quasars also have a higher mass accretion rate and thus a greater number of disk inhomogeneities, resulting in smaller flux variability [Simm et al., 2016]. The inhomogeneous disk model was consistent with mean SDSS spectral analysis in Ruan et al., 2014 but was not a preferred explanation for Kokubo, 2015.

Both [Rumbaugh et al., 2018] (with the Dark Energy Survey SF study) and  Sun et al., 2018  (with a low-$z$ subsample of S82 SDSS quasars) confirmed the anticorrelation between quasar variability and luminosity. However, Graham et al., 2020 did not find support for this trend with the sample of extremely variable quasars (EVQs) in the CRTS dataset,  but when selecting for lower-luminosity sources ($M_{V} < -23$), the anticorrelation is recovered. This agrees with an interpretation that a dwindling fuel supply may correspond to higher variability. Furthermore, Sánchez-Sáez et al., 2018 combined the SDSS spectra with the 5 yr light curves of 2345 quasars obtained with the Quasar Equatorial Survey Team (QUEST)-La Silla AGN Variability Survey, and  using the Bayesian parametrization of SF [Schmidt et al., 2010], they also found that the amplitude of variability $A$ is anticorrelated with rest-frame emission wavelength  and Eddington ratio (also see Simm et al., 2016, Rakshit & Stalin, 2017). See Table  <ref> for a comparison of published studies correlating the statistical measures of observed quasar light curve variability (amplitude/timescale),  and the physical quasar properties (black hole mass, bolometric luminosity).

Indeed, $f_{Edd}$ is a proxy for the strength of accretion, which, together with orientation, may be the key to explaining the quasar main sequence (QMS; Shen & Ho, 2014, Marziani et al., 2018). The QMS, defined by so-called Eigenvector-1, is the anticorrelation between the broad-line Feii emission and the strength of the narrow Oiii ($5007$ $\mbox{\AA}$) line [Wang et al., 1996]. An analysis of quasar clustering by Shen & Ho, 2014, later confirmed by Sun & Shen, 2015 with measurements of  black hole mass from the quasar host galaxy stellar dispersion [Ferrarese & Merritt, 2000, Kormendy & Ho, 2013], showed that the entire diversity of quasars in the QMS can be explained by the variation in accretion (affecting $R_{\mathrm{Fe  II}}$, the ratio of the  Feii equivalent width between $4435$ and $4685$ $\mbox{\AA}$ and H$\beta$) or orientation effects (affecting the FWHM of the H$\beta$). However, Panda et al., 2019, Panda et al., 2019 found that these are insufficient, and variations in metallicity, as well as a range of cloud densities and turbulences, are required.  [Jiang et al., 2016] also found that metallicity, and in particular the iron opacity bump, may have a strong influence on the stability of an accretion disk, and thus linking metallicity to AGN variability. This is also consistent with the findings of [Sun et al., 2018]: quasars with high  Feii strength have higher metallicity and have more stable disks.

§.§ Variability Timescales

In the era of changing-look active galaxies (including initially distinct classes of changing-look quasars [LaMassa et al., 2015, MacLeod et al., 2019], changing-look AGNs (CLAGN, see Bianchi et al., 2009, Risaliti et al., 2009, Marchese et al., 2012), and changing-look LINERS [Frederick et al., 2019], to name a few), there is a revived interest in possibly linking the behavior of stellar-sized accreting systems (e.g., black hole binaries) to that of galactic-scale systems (e.g., AGNs, QSOs, LINERS; Noda & Done, 2018, Ruan et al., 2019).

Several relevant timescales are involved, and there are various interlinked mechanisms that could drive the variability. A standard optically thick, geometrically thin, $\alpha$-disk model has a hierarchy of timescales: dynamical, thermal, front, and viscous, with   $t_{\mathrm{dyn}} < t_{\mathrm{th}} < t_{\mathrm{front}}  < t_{\mathrm{visc}} $ [Netzer, 2013, Frank et al., 2002]. We proceed to briefly describe  each timescale, concluding with our interpretation of the mechanism that could drive the variability observed from the data.

The dynamical, or gas orbital, timescale is simply  an inverse of the Keplerian orbital angular frequency $ \Omega$  at radius $R$:

\begin{equation}
t_{\mathrm{dyn}} {\sim}  1 / \Omega = \left( \frac{GM}{R^{3}}\right)^{-1/2}
\end{equation}

The main parameter  describing the accretion disk is $\alpha$, the ratio of the (vertically averaged) total stress to thermal (vertically averaged) pressure:

\begin{equation}
\alpha= \frac{\langle \tau_{r\varphi}  \rangle_{z} }{\langle P \rangle _{z}}
\end{equation}

After [Lasota, 2016],  the hydrodynamical stress tensor (corresponding to  kinematic viscosity $\nu$) is

\begin{equation}
\tau_{r\varphi } = \rho \nu \frac{\partial v_{\varphi}}{\partial R} = \rho \nu \frac{d \Omega}{d \ln{R}} = \frac{3 \rho \nu \Omega}{2}
\end{equation}

so  with  $c_{s}$ (local sound speed) at radius $R$ (isothermal sound speed is $c_{s} = \sqrt{P/\rho}$),

\begin{equation}
\alpha  =   \frac{3 \rho \nu \Omega}{2 P} =  \frac{3 \Omega \nu}{2 c_{s}^{2}}
\end{equation}

This means that a smaller $\alpha$ corresponds to less viscous disks.

The thermal timescale, related to the time needed for readjustment to the thermal equilibrium (derived in detail in Frank et al., 2002), is the ratio of heat content per unit disk area to dissipation rate per unit disk area: $(dE / A) / (dE/dt /  A) = dt $.  The heat content per unit volume is ${\sim} \rho k T / \mu m_{p} {\sim} \rho c_{s}^{2}$, and the heat content per unit area is  ${\sim} \rho c_{s}^{2} / h {\sim} \Sigma c_{s}^{2}$. Meanwhile, the dissipation rate per unit area, $D(R)$, is

\begin{equation}
D(R) = \frac{9}{8} \nu \Sigma R^{-3} G M
\end{equation}

(eq. 4.30 in Frank et al., 2002), so :

\begin{equation}
t_{\mathrm{th}} {\sim} \frac{c_{s}^{2}R^{3}}{G M \nu } = \frac{c_{s}^{2}}{\nu \Omega} = \frac{t_{\mathrm{dyn}}}{\alpha}
\end{equation}

Thus, if the disk is inviscid ($\nu \rightarrow 0$), then $t_{\mathrm{th}}\rightarrow\infty$; i.e., there is no contact with adjacent disk elements.

The cooling and heating fronts propagate through the disk at  $\alpha c_{s} $ [Hameury et al., 2009]; in that description,  with no viscosity, there is no communication between neighboring disk annuli and thus no front propagation [Balbus & Hawley, 1998, Balbus, 2003]. Following [Stern et al., 2018], if we define the disk aspect ratio as $h/R$ with a disk height $h = c_{s} / \Omega$, the characteristic time for front propagation is

\begin{equation}
t_{\mathrm{front}} {\sim} (h/R) ^ {-1} t_{\mathrm{th}}
\end{equation}

The viscous timescale is the characteristic time it would take for a parcel of material to undergo a radial transport due to the viscous torques from the radius $R$ to the black hole [Czerny, 2006]. Note that while viscosity probably has a magnetic origin [Eardley & Lightman, 1975, Grzędzielski et al., 2017], in this simplistic order-of-magnitude estimate, we use a hydrodynamical description of accretion flow.  With $\nu = \eta / \rho$ (kinematic viscosity being the ratio of dynamical viscosity to density), [Frank et al., 2002] showed (chapter 5.2) that

\begin{equation}
t_{\mathrm{visc}} {\sim} R^{2} / \nu {\sim}  R / v_{R} = (h/R)^{-2} t_{\mathrm{th}}
\end{equation}

We can parametrize each timescale for a black hole mass $M_{\mathrm{BH}} = 10^{8} M_{\odot}$, at $R {\sim} 150 r_{g}$, with the gravitational radius $r_{g} = GM_{\mathrm{BH}} / c^{2} {\sim} 4 \, \mathrm{au}$, using Equations (5)-(8) in [Stern et al., 2018]:

\begin{equation}
t_{\mathrm{dyn}} {\sim} 10  \mathrm{days} \left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}} \right)
\left( \frac{R}{150 r_{g}}\right) ^{3/2}
\end{equation}

\begin{equation}
t_{\mathrm{th}}   {\sim} 1 \,\mathrm{year} \left( \frac{\alpha}{0.03}\right)^{-1}
\left( \frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right) \left( \frac{R}{150 r_{g}}\right)^{3/2}
\end{equation}

\begin{eqnarray}
t_{\mathrm{front}} {\sim} 20 \,\mathrm{years} \left( \frac{h/R}{0.05}\right)^{-1}   \left( \frac{\alpha}{0.03}\right)^{-1}  \nonumber  \\
\left( \frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)     \left( \frac{R}{150 r_{g}}\right) ^{3/2}
\end{eqnarray}

\begin{eqnarray}
t_{\mathrm{visc}}  {\sim} 400 \, \mathrm{years} \left( \frac{h/R}{0.05}\right)^{-2}   \left( \frac{\alpha}{0.03}\right)^{-1} \nonumber  \\
\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}} \right)     \left( \frac{R}{150 r_{g}}\right) ^{3/2}
\end{eqnarray}

In summary,  of the  considered timescales, only the thermal and dynamical are short enough to be related to  the observed short-term stochastic variability. It may be that the variability on the scale of days is driven by local changes and that on the longer scale (perhaps hundreds of days) by a different mechanism [Kokubo, 2015]. The other time scales may be more related to the dramatic changes in brightness of the continuum as observed in changing-look AGNs. Indeed, Noda & Done, 2018 favored a change in mass accretion rate, followed by a propagation of the cooling front [Simm et al., 2016, Lawrence, 2018]. Noda & Done, 2018 also suggested that perhaps some short-term variability could be related to the amount of the disk swept by the thermal front propagation due to hydrogen ionization instability, similar to white dwarf systems (see also Ross et al., 2018, Ruan et al., 2019, Śniegowska & Czerny, 2019).

The variability on a several-year timescale could also be explained by the X-ray reprocessing model [Kokubo, 2015, Kubota & Done, 2018], assuming that the AGN UV–optical variability is a result of reprocessing of X-ray or far-UV emission [Krolik et al., 1991].  The idea of X-ray reprocessing  over time has gained more and more support, with evidence from simultaneous X-ray–UV–optical AGN time series [Edelson et al., 2014, McHardy et al., 2018, Zhu et al., 2018]. In particular, the accretion disk blackbody emission is insufficient to explain the broadband AGN spectrum. The total spectral energy distribution with a soft X-ray excess and a hard X-ray tail requires additional sources of emission. A recent model by Kubota & Done, 2018 divides the flow into blackbody emission, a warm Comptonization region (the disk), and a hard X-ray hot Comptonization component (the corona, or a hot material filling the region close to the black hole below the truncation radius).  Since the soft X-rays are correlated with the hard X-rays, at least part of the picture consists of reflection or reprocessing of hard X-rays by the disk [Lawrence, 2018]. This model predicts an increase of variability amplitude (SF$_{\infty}$) with $M_{\mathrm{BH}}$, and adds the insight that the observed slope is due to changes in accretion rate $\dot{m}$, explaining that smaller $\dot{m}$ corresponds to the highest variability. This qualitatively agrees with the picture that a dwindling fuel supply makes the flow more variable. Previous worries about X-ray reprocessing concerned the seemingly insufficient solid angle subtended by the source of the hard X-rays to cause the observed soft X-ray and optical response. This is addressed by realizing that reprocessing could be taking place in the extended region [Gardner & Done, 2017], such as an inflated inner disk (corresponding to a warm Comptonizing region in Kubota & Done, 2018, or even the BLR region, serving as an additional `complex reprocessor'[McHardy et al., 2018]. Also, for Panda et al., 2019, a warm corona helps decrease the dependence of $R_{\mathrm{Fe  II}}$ on $f_{\mathrm{Edd}}$.

Thus, while CLAGN may be related to the state change to Advection-Dominated Accretion Flow [Śniegowska & Czerny, 2019], similar to that of X-ray binaries [Noda & Done, 2018, Ruan et al., 2019], with cooling and heating fronts [Ross et al., 2018], the short-timescale variability requires approximately three distinct emission regions [Kubota & Done, 2018] with an extended reprocessor (such as a diffuse, hot, puffed-up inner disk and BLR; McHardy et al., 2018) that reverberates the rapid hard X-ray variability in soft X-rays to the optical via UV [Fausnaugh et al., 2018]. Some emission (especially soft X-rays) seems to require the warm Comptonizing corona [Kubota & Done, 2018]. The warm corona, coupled with metallicity changes and variation in turbulence level and cloud density, also helps explain the QMS in the optical [Panda et al., 2019, Panda et al., 2019]. Finally, the Kubota & Done, 2018 model, apart from being consistent with other mechanisms [Lawrence, 2018, McHardy et al., 2018, Panda et al., 2019, Ross et al., 2018, Ruan et al., 2019, Śniegowska & Czerny, 2019], explains the observed correlation of variability amplitude with black hole mass as corresponding to variations in mass accretion rate.

§ SUMMARY AND CONCLUSIONS

We model the optical variability of ${\sim} 9000$ S82 quasars as the DRW [Kelly et al., 2009]. The DRW is a GP, described by two parameters: characteristic timescale $\tau$ (representing the decorrelation timescale, or light-curve smoothness) and the asymptotic amplitude SF$_{\infty}$ (which relates to the amplitude of variability). We fit observed and simulated light curves with celerite - a fast GP solver [Foreman-Mackey et al., 2017]. By simulating and fitting DRW light curves, we explore the impact of the ratio of input timescale and the light-curve baseline. We find that the light-curve length needs to be several times larger than the input timescale to allow unbiased timescale retrieval, confirming K17. Motivated by this result, we consider extending SDSS with PS1, PTF, CRTS, and ZTF data.  We calculate appropriate photometric offsets (color terms) to relate PS1 $gri$, PTF $gR$,  CRTS $V$, and ZTF $r$ to the SDSS $r$ band. However, due to larger photometric uncertainties of PTF, ZTF, and CRTS at the faint magnitudes of SDSS quasars, we decided to use only PS1 $r$-band data. Furthermore, the SDSS and PS1 $r$ bands are sufficiently similar that no photometric transformation is required.  Thus, by extending the SDSS $r$-band light curves with PS1 DR2 $r$-band data, we improve upon the fidelity of recovered DRW parameters (e.g., in Figure <ref>, showing a simulated population of $\tau=575$ days, the rms of $\log{(\tau_{fit}/ \mathrm{baseline})}$ decreases from 1.75 dex with SDSS to 1.5 dex with PS1, and in the future, with the inclusion of ZTF and LSST data, it will decrease to $\sim 1$ dex).

We identify 40 objects that exhibit a tenfold increase in variability timescale when using the SDSS–PS1 data set, as compared to the timescale inferred from SDSS alone. Their light curves show characteristics of changing-look quasars (magnitude difference larger than 0.5 mag; MacLeod et al., 2016). Of these, five are confirmed in the literature  [MacLeod et al., 2019, LaMassa et al., 2015]. We recommend spectroscopic follow-up and further monitoring of the brightest targets (see Appendix <ref>).

We investigate the correlation of quasar physical properties, such as black hole mass $M_{\mathrm{BH}}$ and absolute $i$-band magnitude $M_{i}$,  with DRW model parameters. The SDSS–PS1 data, coupled with the Shen et al., 2011 quasar catalog, imply that the damping timescale $\tau$  is correlated with $M_{\mathrm{BH}}$ with a power-law index of $0.141 \pm 0.019$ and almost independent of quasar bolometric luminosity as in M10, Wilhite et al., 2008, and Vanden Berk et al., 2004. The asymptotic variability amplitude SF$_{\infty}$ is correlated with $M_{i}$ (i.e. anticorrelated with luminosity) with a power-law index of $0.118 \pm 0.003$, and correlated with $M_{\mathrm{BH}}$ with a slope of $0.118\pm 0.008$. This can be explained if the driving variable was the Eddington ratio, $f_{Edd}$ [Wilhite et al., 2008]. Indeed, there is an anticorrelation of SF$_{\infty}$ and $f_{\mathrm{Edd}}$, with a power-law slope of $-0.207 \pm 0.031$ (similar to M10).  As suggested by Kubota & Done, 2018, this gradient of SF$_{\infty}$ in the plane of $M_{\mathrm{BH}}$ vs $M_{i}$ could be explained if the lower mass accretion rate corresponds to higher variability, so that when the supply of fuel decreases, the flow becomes less stable and more clumpy and inhomogeneous  [Rakshit & Stalin, 2017, Kokubo, 2015, Dexter & Agol, 2011].  This is also consistent with the X-ray reprocessing model, whereby the hard X-ray variability of the inner disk is reflected/reprocessed by the extended warm Comptonization region (inflated disk) and perhaps a complex reprocessor, including the clouds of the BLR [Kubota & Done, 2018, Panda et al., 2019]. Changes on recovered timescales are too fast to be driven by changes in disk viscosity or thermal front propagation alone; a thermal or dynamical timescale of response to the changes in X-ray emission seems most consistent with our results  [Stern et al., 2018].

More data extending the light curves would help improve the DRW fit coefficients, potentially decreasing the scatter in the observed correlations - for instance, for simulated $\tau$ recovery, an improvement of the rms by a factor of 1.8 (from 0.322 for SDSS only to 0.182 with combined SDSS–PS1–ZTF–LSST). Moreover, given that the uncertainty in black hole mass is one of the biggest sources of error, better measurements of quasar properties would be of high utility [Shen et al., 2011].  This will be possible with the upcoming AGN reverberation mapping campaigns (e.g., SDSS-V black hole mapper), providing better calibration for line width–based methods of estimating black hole masses [Kollmeier et al., 2017]. All quasars in this study were spectroscopically confirmed, but some spectra had low signal-to-noise ratios, resulting in a higher likelihood of incorrect redshift measurement. Better spectroscopy and follow-up of S82 quasars, afforded by SDSS-V panoptic spectroscopy, would not only help improve on the spectrum-based  properties (redshift, absolute magnitude, black hole masses) but also allow the study of spectral changes and further new CLAGN discoveries [MacLeod et al., 2019].

If this study were to be expanded onto a sample of quasars with good photometry over sufficiently long baselines but lacking spectral information, the required physical information on quasars could be obtained by indirect methods of estimating the coarse spectral information from broadband photometry [Kozłowski, 2015]. This would benefit from better catalogs of existing spectroscopically confirmed quasars (SDSS DR14) to improve the calibration, as well as better methods of estimating the redshift based on photometry alone (e.g., photo-$z$; Richards et al., 2015, Yang et al., 2017, Curran & Moss, 2019, Jin et al., 2019).  This will be possible in the short term with the ZTF [Bellm et al., 2019] and in the long term with the LSST [Ivezić et al., 2019]. Occasional coverage adding a few epochs to some quasars may be possible with other surveys (e.g., TESS; Ricker et al., 2014), but to improve the statistics of an entire sample of S82 quasars would require longer baselines. Combining SDSS and PS1 with LSST would provide an unprecedented 35 yr baseline, which, assuming timescales below 1000 days, is over 10 times longer, allowing unbiased DRW parameter retrieval. This, coupled with correlations with quasar properties, would provide an estimate of black hole masses and bolometric luminosities for millions of quasars [Ivezić et al., 2019].

§ ACKNOWLEDGMENTS

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. The SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at
the University of Utah. The SDSS website is www.sdss.org. The SDSS-IV is managed by the Astrophysical Research Consortium for the  Participating Institutions of the SDSS Collaboration, including the  Brazilian Participation Group, the Carnegie Institution for Science, the Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, the Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, the Kavli Institute for the Physics and Mathematics of the Universe (IPMU) /
University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory,
Leibniz Institut für Astrophysik Potsdam (AIP),  Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, the University of Notre Dame, Observatário Nacional / MCTI, the Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México, the University of Arizona, the University of Colorado Boulder, the University of Oxford, the University of Portsmouth,
the University of Utah, the University of Virginia, the University of Washington, the University of Wisconsin, Vanderbilt University, and Yale University.

The CSS survey is funded by the National Aeronautics and Space Administration under grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program.  The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182.

Based on observations obtained with the Samuel Oschin 48 inch telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. The ZTF is supported by the National Science Foundation under grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

§ MEASURING QUASAR PROPERTIES

In this work, we employ black hole masses, bolometric luminosities, and $k$-corrections from the [Shen et al., 2011] catalog based on single-epoch SDSS spectra. Here we explain the choices made in the difficult art of estimating each of these quasar physical properties.

It is nontrivial to measure the mass of black holes living in the centers of active galaxies, even provided with a detailed spectrum. The most common  approach to estimate black hole masses in AGNs is to assume that the BLR is virialized,

\begin{equation}
M_{\mathrm{BH}} = f \frac{ R\Delta V^{2} }{G} = f M_{\mathrm{vir}}
\end{equation}

where $f$ is a constant of order unity, $R$ is the size of the BLR (estimated from  emission line lag $\Delta t$ as $R = c \Delta t$), $\Delta V$ is the virial velocity, and $G$ is the gravitational constant [Shen et al., 2008].  From reverberation mapping (RM) studies (e.g., Shen et al., 2019), we know that continuum luminosity $L$ is related to the size of the BLR region as $R \propto L^{\gamma}$ [Vestergaard & Peterson, 2006], with $\gamma$ very close to $1/2$ (e.g., Bentz et al., 2009 found from RM studies $\gamma = 0.519 \pm 0.06$). Thus, we find  that $R \Delta V^{2} \propto L^{\gamma} \Delta V^{2} \equiv \mu$. The virial velocity $\Delta V$ is usually estimated from the width of the broad emission lines (or line dispersion).  In the absence of a quasar spectrum, there are alternative methods using a conversion of the broadband photometry into monochromatic fluxes in the vicinity of reverberating lines (e.g., Kozłowski, 2015, used in Kozłowski, 2017 to estimate black hole masses for 280,000 AGNs). Depending on the redshift, different rest-frame calibrated emission lines shift into the observed passband: broad H$\alpha$ at  $6562\mbox{\AA}$, H$\beta$ at  $5100 \mbox{\AA}$, Mg ii at $2800\mbox{\AA}$, and C iv at $1350\mbox{\AA}$ (see Figure 7 in Shen et al., 2019, and Vestergaard, 2002). Some authors even separately consider  C iv- and  Mg ii-based  black hole mass estimates. We refer the reader to Shen et al., 2008, who described in detail the various biases and inherent assumptions of virial black hole mass measurements.

Another important quasar property, bolometric luminosity, is most often estimated from the absolute $i$-band magnitude, $M_{i}$ (see Shen et al., 2008, Figure 2). Here $M_{i}$ is derived from the observed $i$-band magnitude by correcting for Galactic extinction and the fact that at different redshifts, different portions of the spectral energy distribution are observed by the telescope filter bandpass. The latter, known as $k$-correction $K(z)$ [Oke & Sandage, 1968],  is defined as $m_{\mathrm{intrinsic}} = m_{\mathrm{observed}} - K(z)$. In the early 2000s, the common approach was to $k$-correct to a redshift of zero, but as [Richards et al., 2006] pointed out, since the distribution of quasars peaks at a redshift of 2, for most quasars, correcting to a redshift of zero required shifting the observed spectrum into the far infrared. Moreover, the procedure was to correct separately for the continuum and emission line contributions, assuming a particular spectral shape (e.g., power law  $f_{\nu} \propto \nu^{\alpha}$, with $\alpha=-0.5$; see Vanden Berk et al., 2001, Richards et al., 2006, Schneider et al., 2010).  This introduces a larger error for $K(z=0)$ than for $K(z=2)$ if the assumed spectral shape $\alpha=-0.5$ is far from the real spectral index. In early 2010s, after  Richards et al., 2006, Wisotzki, 2000, Blanton et al., 2003,  the practice started shifting toward $k$-correcting to a redshift of 2,  and including custom quasar spectral shapes, as reflected by the content of the [Shen et al., 2011] quasar catalog. Thus, in this study, we use  the absolute $i$-band magnitude $k$-corrected to $z=2$: $M_{i}(z=2)$.

These methods were used to create catalogs of quasar properties derived from spectra. Since quasars are variable at the ${\sim}0.2$ mag level, the ideal is to use a single-epoch calibrated spectrum to estimate the continuum luminosity and find virial black hole masses using relationships based on the monochromatic fluxes and broad line widths described above. A glance at the available quasar catalogs reveals that, given any SDSS data release, there is indeed first a catalog of basic quasar properties (redshift and photometry; e.g., Schneider et al., 2007, Schneider et al., 2010), and more detailed catalogs containing black hole masses and bolometric luminosities  follow (e.g., Shen et al., 2008, Shen et al., 2011). More recently, once the SDSS DR12 Quasar Catalog [Pâris et al., 2017] was released,  K17 followed using SDSS photometry as a proxy for monochromatic luminosities. Chen et al., 2018 added a detailed analysis of continuum luminosities in the  H$\alpha$, H$\beta$ regions for low-redshift quasars. Using the spectra from the Chinese LAMOST survey, Dong et al., 2018 also sought to estimate virial black hole masses, and the results, while consistent with Shen et al., 2011, suffered from the necessity of pegging the noncalibrated spectra to the SDSS photometry, which was taken at a different epoch. Thus, even though the SDSS DR12 Quasar Catalog of Pâris et al., 2018  is the most recent, like Pâris et al., 2017, it lacks black hole masses and bolometric luminosities, and there is no recent work that reanalyzed the spectral data. Therefore, we use black hole mass estimates and monochromatic luminosities from Shen et al., 2011 based directly on single-epoch spectra.

§ CLQSO CANDIDATES

Based on the DRW model  parameters $\tau$, $\sigma$ fitted with celerite using the SDSS and PS1 data, we find that there are quasars for which there is a pronounced difference between $\tau$, $\sigma$  inferred from combined SDSS–PS1 data vs just SDSS. Specifically, Figure <ref> shows that there are objects where $f_{\sigma} = \log_{10}{\left( \sigma_{\mathrm{SDSS-PS1}} / \sigma_{\mathrm{SDSS}} \right)} > 0.4 $ and $f_{\tau} = \log_{10}{\left( \tau_{\mathrm{SDSS-PS1}} / \tau_{\mathrm{SDSS}} \right)} > 1$ (a tenfold increase in $\tau$ and over twofold increase in $\sigma$). Visual inspection of objects simultaneously satisfying $f_{\sigma}> 0.4$ and $f_{\tau}> 1 $ shows that these underwent a significant ($>0.5$ mag) change in brightness between the SDSS (baseline 1998 - 2008) and PS1 DR2 observations (2009-2014; see Figure <ref>). Thus, DRW fitting could also be a way of finding changing-look quasar (and AGN) candidates. Figures <ref>– <ref> show the SDSS–PS1 $r$-band light curves of 40 CLQSO candidates with median PS1 brightnesses larger  than 20.5 mag. The open circles indicate day-averaged epochs (see Sec. <ref>). Table <ref> contains the basic physical parameters for these quasars. Some quasars show a downward trend, like turn-off CLQSOs (e.g., 123909, 1412379, and 1644710), while some are seen in a brightening stage, like turn-on CLQSOs (e.g., 1976348, 221006, 4069419, and 4205621). Quasar 612585 has the largest amount of auxiliary multiwavelength coverage– X-ray from XMM-Newton, UV from the Galaxy Evolution Explorer (GALEX), and IR from the UKIRT Infrared Deep Sky Survey (UKIDSS), VHS, the Wide-field Infrared Survey Explorer (WISE)–and has been analyzed as part of the X-ray-targeted sample of S82 quasars, S82X, by [LaMassa et al., 2016]. Quasar 751557 was previously identified by [MacLeod et al., 2019] as a CLQSO candidate, with detailed Magellan spectroscopy described therein.  Two quasars, 1003694 and 1299803, have WISE data in the S82X catalog [LaMassa et al., 2016]. Quasars 612585 and 3633437 have X-ray detections in the 3XMM DR5 catalog [Rosen et al., 2016], but they have no matches in the Chandra point-source catalog (second release; Evans et al., 2010, Evans & Civano, 2018). There are no matches against the unified radio catalog of [Kimball & Ivezić, 2008], which includes FIRST and NVSS data. We especially recommend spectroscopic follow-up of the brightest targets: 1976348 (mean 17.8 mag, turn-on; top left panel in Figure <ref>) and 2104791 (mean 18.4 mag, turn-off; left column, third row panel in Figure <ref>).

Outliers in the space of the recovered DRW parameters between SDSS and SDSS–PS1, as well as median offsets. Page 1/4 (continued on Figure <ref>).

As Figure <ref>, page 2/4 (continued on Figure <ref>).

As Figure <ref>, page 3/4 (continued on Figure <ref>).

As Figure <ref>, page 4/4.

CLQSO Candidates.
dbID     SDSSJID     $\alpha$     $\delta$     Redshift     $\log_{10}{(L_{\mathrm{Bol}})}$     $\log_{10}{(M_{\mathrm{BH}}/ M_{\odot})}$     $f_{\mathrm{Edd}}$     $\Delta(\mathrm{mag})$     $\Delta(\sigma_{G})$     Median PS1
123909     [SDSS J001626.54+003632.4]001626.54+003632.4     4.111     0.609     3.24     46.57     9.47     -1.0     -0.29     -0.054     20.48

8442     [SDSS J001731.70+004910.1]001731.70+004910.1     4.382     0.819     2.43     46.61     9.09     -0.58     -0.36     -0.107     20.36

4069419     [SDSS J003359.39+000230.0]003359.39+000230.0     8.497     0.042     1.64     45.95     9.05     -1.21     0.22     0.017     20.2

221006     [SDSS J005142.20+002129.0]005142.20+002129.0     12.926     0.358     1.55     45.95     8.24     -0.39     0.21     -0.01     19.97

257776     [SDSS J005513.15-005621.2]005513.15-005621.2     13.805     -0.939     3.61     47.13     9.58     -0.54     -0.32     -0.242     19.53

612585a     [SDSS J010812.00-000516.5]010812.00-000516.5     17.05     -0.088     1.0     45.52     9.06     -1.64     -0.37     0.006     20.5

1003694b     [SDSS J012114.19-010310.8]012114.19-010310.8     20.309     -1.053     1.89     46.59     8.83     -0.34     0.28     0.032     19.11

1299803b     [SDSS J014303.23-004354.0]014303.23-004354.0     25.763     -0.732     0.53     45.78     8.68     -1.0     -0.18     -0.008     18.56

1644710     [SDSS J021259.00-000550.1]021259.00-000550.1     33.246     -0.097     0.81     45.67     8.38     -0.81     -0.28     0.022     19.7

1730482     [SDSS J021529.02-005314.9]021529.02-005314.9     33.871     -0.887     1.37     45.98     8.8     -0.92     0.66     -0.009     19.57

2104791*     [SDSS J022239.83+000022.5]022239.83+000022.5     35.666     0.006     0.99     46.28     9.33     -1.16     -0.31     0.013     18.66

2061101     [SDSS J022505.06+001733.2]022505.06+001733.2     36.271     0.293     2.42     46.38     8.09     0.2     -0.34     -0.016     20.15

2006852     [SDSS J023917.86-001916.8]023917.86-001916.8     39.824     -0.321     1.41     46.07     8.73     -0.76     -0.45     -0.013     20.24

2503955     [SDSS J025316.46+010759.7]025316.46+010759.7     43.319     1.133     1.03     46.31     8.94     -0.73     -0.39     -0.071     19.07

2484608     [SDSS J025654.42-011455.4]025654.42-011455.4     44.227     -1.249     0.54     45.57     8.48     -1.01     -0.43     -0.023     19.65

3052176     [SDSS J030504.07+011324.5]030504.07+011324.5     46.267     1.223     0.61     45.29     9.2     -2.01     -0.35     0.016     20.09

3096136     [SDSS J031401.11+011131.6]031401.11+011131.6     48.505     1.192     1.31     46.06     9.21     -1.25     -0.27     -0.181     20.09

3147102     [SDSS J031846.13-005622.8]031846.13-005622.8     49.692     -0.94     2.12     46.57     9.07     -0.6     0.47     -0.078     19.07

3781306     [SDSS J032745.74+005217.2]032745.74+005217.2     51.941     0.871     1.16     45.77     N/A     N/A     0.19     0.032     19.94

3858587     [SDSS J032825.19-003252.3]032825.19-003252.3     52.105     -0.548     0.77     45.61     8.68     -1.17     0.33     0.031     19.39

3810874     [SDSS J033047.73+004859.4]033047.73+004859.4     52.699     0.816     0.86     45.8     8.35     -0.65     -0.1     -0.177     19.8

3739253     [SDSS J033059.05+010952.0]033059.05+010952.0     52.746     1.164     0.56     45.39     8.16     -0.86     -0.27     0.01     19.72

3791907     [SDSS J033431.17-000904.0]033431.17-000904.0     53.63     -0.151     1.64     46.01     9.12     -1.21     0.28     0.02     20.29

4913626     [SDSS J034512.62+002245.7]034512.62+002245.7     56.303     0.379     0.42     45.5     8.81     -1.41     -0.48     -0.055     19.43

4205621     [SDSS J203932.41-001818.3]203932.41-001818.3     309.885     -0.305     1.58     46.21     8.66     -0.55     0.17     -0.126     19.77

3319354     [SDSS J204952.62+011306.6]204952.62+011306.6     312.469     1.219     1.09     46.14     9.52     -1.48     -0.91     -0.048     20.08

3633437c     [SDSS J205105.02-005847.5]205105.02-005847.5     312.771     -0.98     0.54     45.34     8.47     -1.23     -0.51     0.071     20.34

1835106     [SDSS J215015.05-005331.4]215015.05-005331.4     327.563     -0.892     1.9     46.33     9.21     -0.98     0.69     0.035     19.28

1901056     [SDSS J215055.51-001739.4]215055.51-001739.4     327.731     -0.294     1.54     46.26     8.6     -0.44     -0.19     0.009     20.46

1976348*     [SDSS J215841.40-001507.7]215841.40-001507.7     329.673     -0.252     1.46     46.92     9.39     -0.57     0.39     -0.051     17.56

1446022     [SDSS J220535.23+000756.3]220535.23+000756.3     331.397     0.132     1.69     46.45     9.25     -0.9     0.17     -0.039     19.34

1378415     [SDSS J221347.32+001928.4]221347.32+001928.4     333.447     0.325     2.31     46.41     8.59     -0.29     0.23     -0.113     20.02

1412379     [SDSS J221831.58-004548.9]221831.58-004548.9     334.632     -0.764     1.23     46.15     9.48     -1.43     -0.61     -0.057     20.04

1124333     [SDSS J222918.25-004003.6]222918.25-004003.6     337.326     -0.668     1.16     45.81     8.35     -0.64     -0.28     -0.003     20.12

751557d     [SDSS J225240.37+010958.7]225240.37+010958.7     343.168     1.166     0.53     45.32     8.88     -1.66     0.39     -0.536     19.67

467617     [SDSS J231032.17-011449.5]231032.17-011449.5     347.634     -1.247     1.82     46.03     7.86     0.06     0.23     0.024     20.18

568312     [SDSS J231953.07-010139.0]231953.07-010139.0     349.971     -1.028     1.15     45.6     8.29     -0.79     0.33     0.005     19.79

292959     [SDSS J232030.97-004039.2]232030.97-004039.2     350.129     -0.678     1.72     46.69     9.39     -0.8     -0.22     -0.009     19.07

3976336     [SDSS J235213.27-004326.3]235213.27-004326.3     358.055     -0.724     0.9     45.64     8.85     -1.3     -0.18     -0.08     19.92

3946479     [SDSS J235248.71-001518.4]235248.71-001518.4     358.203     -0.255     1.34     45.8     9.01     -1.3     -0.07     -0.025     20.26

Catalog information from [Shen et al., 2011] concerning DR7 name (dbID), SDSSJID, location $\alpha$ and $\delta$ (in degrees, J2000), distance (spectrum-based redshift),  physical parameters (bolometric luminosity  $L_{\mathrm{Bol}}$ erg s$^{-1}$, black hole mass $M_{\mathrm{BH}}/ M_{\odot} $, Eddington ratio $f_{\mathrm{Edd}} = L_{\mathrm{Bol}} / L_{\mathrm{Edd}}$),  and light-curve properties (the difference between median SDSS and PS1 magnitudes $\Delta(\mathrm{mag})$, the difference in scatter between SDSS and PS1 segments $\Delta(\sigma_{G})$, and the median PS1 magnitude; see Figure <ref>).
*Recommended for follow-up
aS82X [LaMassa et al., 2016], XMM-Newton, GALEX UV, UKIDSS, VHS, WISE
bS82X [LaMassa et al., 2019], WISE
cIn 3XMM DR5 X-ray catalog [Rosen et al., 2016]
dM19, CLQSO candidate, Magellan follow-up

§ MG II VARIABILITY

We searched for the dependency of variability parameters on other physical properties beyond the black hole mass and quasar luminosity employed in Equation <ref>. The quasar optical spectrum has certain strong emission lines, depending on the redshift. In particular, the Mg ii emission line lags the continuum and is less variable [Reichert et al., 1994]. The Mg ii line  is also an important virial black hole mass estimator for quasars; [McLure & Jarvis, 2002] noticed that the FWHM of the Mg ii doublet trails that of the  H$\beta$ line (also see  Shen, 2013). Following  [Ivezić et al., 2004] and [MacLeod et al., 2012],  we investigated the residuals after fitting the model ($f$=SF$_{\infty}$ in Equation  <ref>) to the SDSS–PS1 data as a function of wavelength. The left panel of Figure <ref> shows the SF$_{\infty}$model residuals in the $\lambda_{\mathrm{RF}}$ vs $\tau_{\mathrm{RF}}$ space. The right hand side panel shows the residuals marginalized along $\tau_{\mathrm{RF}}$. The decrement around 2800 $\mbox{\AA}$ in the right panel is more pronounced when using combined SDSS–PS1 light curves. We see that both aggregates based on raw data (blue dots) or on medians (black dots) agree; the two statistical methods show a small (${\sim}5\%$), but statistically significant detection of a depression in the 2800 $\mbox{\AA}$  region, as expected from theory.

The residuals (z) after fitting Equation <ref> to the quasar SF$_{\infty}$, using SDSS–PS1 quasar light curves. The left panel shows the residuals as a function of  $\lambda_{\mathrm{RF}}$ and $\tau_{\mathrm{RF}}$. The right panel shows the median of the residuals marginalized along $\tau_{\mathrm{RF}}$, with an error on the median based on the scatter in the $\lambda_{\mathrm{RF}}$ bins, $\sigma_{median} = 1.25 \mathrm{RMS}_{\mathrm{bin}} / \sqrt{N_{\mathrm{bin}}}$, where $\mathrm{RMS}_{\mathrm{bin}}$ is estimated from the robust Gaussian interquartile-based $\sigma_{G}$ [Ivezić et al., 2014]. Black dots are medians based on the averages plotted on the left ($N_{\mathrm{bin}}$ is between 7 and 20, depending on the number of nonempty bins in the left panel), and blue dots are medians based on the raw data ($N_{\mathrm{bin}}$ varies as a function of  $\lambda_{\mathrm{RF}}$, from a few to 800, peaking at about 2500 $\mbox{\AA}$, reflecting the fact that the quasar distribution peaks at redshift $z=2$).  We overplot the composite quasar spectrum from Vanden Berk et al., 2001, and mark the location of the Mg ii 2798 $\mbox{\AA}$ line with a vertical dashed red line.

[AlSayyad, 2016]
AlSayyad, Y. 2016, PhD thesis, University of Washington.

[Annis et al., 2014]
Annis, J., Soares-Santos, M., Strauss, M. A., et al. 2014, , 794,
120, 10.1088/0004-637X/794/2/120

[Aranzana et al., 2018]
Aranzana, E., Körding, E., Uttley, P., Scaringi, S., & Bloemen,
S. 2018, , 476, 2501, 10.1093/mnras/sty413

[Babić et al., 2007]
Babić, A., Miller, L., Jarvis, M. J., et al. 2007, , 474, 755,

[Balbus, 2003]
Balbus, S. A. 2003, , 41, 555,

[Balbus & Hawley, 1998]
Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1,

[Bauer et al., 2009]
Bauer, A., Baltay, C., Coppi, P., et al. 2009, , 696, 1241,

[Bellm et al., 2019]
Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, , 131,
018002, 10.1088/1538-3873/aaecbe

[Bentz et al., 2009]
Bentz, M. C., Peterson, B. M., Netzer, H., Pogge, R. W., &
Vestergaard, M. 2009, , 697, 160, 10.1088/0004-637X/697/1/160

[Bianchi et al., 2009]
Bianchi, S., Piconcelli, E., Chiaberge, M., et al. 2009, , 695,
781, 10.1088/0004-637X/695/1/781

[Blanchard et al., 2017]
Blanchard, P. K., Nicholl, M., Berger, E., et al. 2017, The Astrophysical
Journal, 843, 106, 10.3847/1538-4357/aa77f7

[Blanton et al., 2003]
Blanton, M. R., Lin, H., Lupton, R. H., et al. 2003, , 125, 2276,

[Borucki et al., 2010]
Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977,

[Cackett et al., 2015]
Cackett, E. M., Gültekin, K., Bentz, M. C., et al. 2015, , 810,
86, 10.1088/0004-637X/810/2/86

[Cai et al., 2016]
Cai, Z.-Y., Wang, J.-X., Gu, W.-M., et al. 2016, , 826, 7,

[Cai et al., 2018]
Cai, Z.-Y., Wang, J.-X., Zhu, F.-F., et al. 2018, The Astrophysical Journal,
855, 117, 10.3847/1538-4357/aab091

[Caplar et al., 2017]
Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2017, , 834, 111,

[Chambers et al., 2016]
Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv
e-prints, arXiv:1612.05560.

[Charisi et al., 2016]
Charisi, M., Bartos, I., Haiman, Z., et al. 2016, , 463, 2145,

[Chen et al., 2018]
Chen, Z.-F., Pan, D.-S., Pang, T.-T., & Huang, Y. 2018, , 234,
16, 10.3847/1538-4365/aa9d90

[Curran & Moss, 2019]
Curran, S. J., & Moss, J. P. 2019, , 629, A56,

[Czerny, 2006]
Czerny, B. 2006, in Astronomical Society of the Pacific Conference Series,
Vol. 360, AGN Variability from X-Rays to Radio Waves, ed. C. M. Gaskell,
I. M. McHardy, B. M. Peterson, & S. G. Sergeev, 265

[Dexter & Agol, 2011]
Dexter, J., & Agol, E. 2011, , 727, L24,

[Dexter & Begelman, 2019]
Dexter, J., & Begelman, M. C. 2019, , 483, L17,

[Dong et al., 2018]
Dong, X. Y., Wu, X.-B., Ai, Y. L., et al. 2018, , 155, 189,

[Drake et al., 2009]
Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, , 696,
870, 10.1088/0004-637X/696/1/870

[Eardley & Lightman, 1975]
Eardley, D. M., & Lightman, A. P. 1975, , 200, 187,

[Edelson et al., 2014]
Edelson, R., Vaughan, S., Malkan, M., et al. 2014, , 795, 2,

[Edelson et al., 2015]
Edelson, R., Gelbord, J. M., Horne, K., et al. 2015, , 806, 129,

[Elitzur et al., 2014]
Elitzur, M., Ho, L. C., & Trump, J. R. 2014, , 438, 3340,

[Evans & Civano, 2018]
Evans, I. N., & Civano, F. 2018, Astronomy and Geophysics, 59, 2.17,

[Evans et al., 2010]
Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, The
Astrophysical Journal Supplement Series, 189, 37,

[Fausnaugh et al., 2016]
Fausnaugh, M. M., Denney, K. D., Barth, A. J., et al. 2016, , 821,
56, 10.3847/0004-637X/821/1/56

[Fausnaugh et al., 2018]
Fausnaugh, M. M., Starkey, D. A., Horne, K., et al. 2018, , 854,
107, 10.3847/1538-4357/aaaa2b

[Ferrarese & Merritt, 2000]
Ferrarese, L., & Merritt, D. 2000, The Astrophysical Journal, 539, L9,

[Flewelling, 2018]
Flewelling, H. 2018, in American Astronomical Society Meeting Abstracts, Vol.
231, American Astronomical Society Meeting Abstracts 231, 436.01

[Flewelling et al., 2020]
Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2020, ,
251, 7, 10.3847/1538-4365/abb82d

[Foreman-Mackey et al., 2017]
Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, ,
154, 220, 10.3847/1538-3881/aa9332

[Frank et al., 2002]
Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in
Astrophysics: Third Edition (Cambridge University Press), 398

[Frederick et al., 2019]
Frederick, S., Gezari, S., Graham, M. J., et al. 2019, , 883, 31,

[Gardner & Done, 2017]
Gardner, E., & Done, C. 2017, Monthly Notices of the Royal Astronomical
Society, 470, 3591, 10.1093/mnras/stx946

[Gezari et al., 2017]
Gezari, S., Hung, T., Cenko, S. B., et al. 2017, The Astrophysical Journal,
835, 144, 10.3847/1538-4357/835/2/144

[Graham et al., 2015]
Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, , 518,
74, 10.1038/nature14143

[Graham et al., 2020]
Graham, M. J., Ross, N. P., Stern, D., et al. 2020, , 491, 4925,

[Grzędzielski et al., 2017]
Grzędzielski, M., Janiuk, A., Czerny, B., & Wu, Q. 2017, ,
603, A110, 10.1051/0004-6361/201629672

[Guo et al., 2017]
Guo, H., Wang, J., Cai, Z., & Sun, M. 2017, , 847, 132,

[Guo et al., 2016]
Guo, H., Malkan, M. A., Gu, M., et al. 2016, The Astrophysical Journal, 826,
186, 10.3847/0004-637x/826/2/186

[Guo et al., 2020]
Guo, H., Shen, Y., He, Z., et al. 2020, , 888, 58,

[Hameury et al., 2009]
Hameury, J.-M., Viallet, M., & Lasota, J.-P. 2009, , 496, 413,

[Hawkins, 2007]
Hawkins, M. R. S. 2007, Astronomy and Astrophysics, 462, 581,

[Hernitschek et al., 2016]
Hernitschek, N., Schlafly, E. F., Sesar, B., et al. 2016, The Astrophysical
Journal, 817, 73

[Hopkins et al., 2005]
Hopkins, P. F., Hernquist, L., Martini, P., et al. 2005, , 625,
L71, 10.1086/431146

[Ivezić et al., 2014]
Ivezić, Ž., Connolly, A. J., VanderPlas, J. T., & Gray, A.
2014, Statistics, Data Mining, and Machine Learning in Astronomy

[Ivezić et al., 2004]
Ivezić, Ž., Lupton, R. H., Juric, M., et al. 2004, in IAU
Symposium, Vol. 222, The Interplay Among Black Holes, Stars and ISM in
Galactic Nuclei, ed. T. Storchi-Bergmann, L. C. Ho, & H. R. Schmitt,

[Ivezić et al., 2007]
Ivezić, Ž., Smith, J. A., Miknaitis, G., et al. 2007, ,
134, 973, 10.1086/519976

[Ivezić et al., 2019]
Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, , 873,
111, 10.3847/1538-4357/ab042c

[Jiang et al., 2016]
Jiang, Y.-F., Davis, S. W., & Stone, J. M. 2016, , 827, 10,

[Jin et al., 2019]
Jin, X., Zhang, Y., Zhang, J., et al. 2019, , 485, 4539,

[Kasliwal et al., 2015]
Kasliwal, V. P., Vogeley, M. S., & Richards, G. T. 2015, , 451,
4328, 10.1093/mnras/stv1230

[Kelly, 2007]
Kelly, B. C. 2007, , 665, 1489, 10.1086/519947

[Kelly et al., 2009]
Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, The Astrophysical
Journal, 698, 895

[Kelly et al., 2014]
Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., &
Uttley, P. 2014, , 788, 33, 10.1088/0004-637X/788/1/33

[Kelly et al., 2013]
Kelly, B. C., Treu, T., Malkan, M., Pancoast, A., & Woo, J.-H. 2013,
, 779, 187, 10.1088/0004-637X/779/2/187

[Kimball & Ivezić, 2008]
Kimball, A. E., & Ivezić, Ž. 2008, The Astronomical Journal,
136, 684, 10.1088/0004-6256/136/2/684

[Kokubo, 2015]
Kokubo, M. 2015, , 449, 94, 10.1093/mnras/stv241

[Kollmeier et al., 2017]
Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv
e-prints, arXiv:1711.03234.

[Kormendy & Ho, 2013]
Kormendy, J., & Ho, L. C. 2013, , 51, 511,

[Kozłowski, 2015]
Kozłowski, S. 2015, , 65, 251.

[Kozłowski, 2016]
Kozłowski, S. 2016, , 459, 2787, 10.1093/mnras/stw819

[Kozłowski, 2017]
Kozłowski, S. 2017, , 228, 9, 10.3847/1538-4365/228/1/9

[Kozłowski et al., 2016]
Kozłowski, S., Kochanek, C. S., Ashby, M. L. N., et al. 2016, ,
817, 119, 10.3847/0004-637X/817/2/119

[Kozłowski et al., 2010]
Kozłowski, S., Kochanek, C. S., Udalski, A., et al. 2010, , 708, 927

[Kozłowski, Szymon, 2017]
Kozłowski, Szymon. 2017, , 597, A128,

[Krolik et al., 1991]
Krolik, J. H., Horne, K., Kallman, T. R., et al. 1991, , 371, 541,

[Kubota & Done, 2018]
Kubota, A., & Done, C. 2018, , 480, 1247,

[Labita et al., 2009]
Labita, M., Decarli, R., Treves, A., & Falomo, R. 2009, , 399,
2099, 10.1111/j.1365-2966.2009.15400.x

[LaMassa et al., 2019]
LaMassa, S. M., Georgakakis, A., Vivek, M., et al. 2019, , 876, 50,

[LaMassa et al., 2015]
LaMassa, S. M., Cales, S., Moran, E. C., et al. 2015, , 800, 144,

[LaMassa et al., 2016]
LaMassa, S. M., Urry, C. M., Cappelluti, N., et al. 2016, , 817,
172, 10.3847/0004-637X/817/2/172

[Lasota, 2016]
Lasota, J.-P. 2016, in Astrophysics and Space Science Library, Vol. 440,
Astrophysics of Black Holes: From Fundamental Aspects to Latest Developments,
ed. C. Bambi, 1

[Lawrence, 2018]
Lawrence, A. 2018, Nature Astronomy, 2, 102,

[Li & Cao, 2008]
Li, S.-L., & Cao, X. 2008, , 387, L41,

[Li et al., 2018]
Li, Z., McGreer, I. D., Wu, X.-B., Fan, X., & Yang, Q. 2018, ,
861, 6, 10.3847/1538-4357/aac6ce

[Lira et al., 2015]
Lira, P., Arévalo, P., Uttley, P., McHardy, I. M. M., & Videla,
L. 2015, , 454, 368, 10.1093/mnras/stv1945

[MacLeod et al., 2010]
MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, The
Astrophysical Journal, 721, 1014

[MacLeod et al., 2011]
MacLeod, C. L., Brooks, K., Ivezić, Ž., et al. 2011, The
Astrophysical Journal, 728, 26

[MacLeod et al., 2012]
MacLeod, C. L., Ivezić, Ž., Sesar, B., et al. 2012, The
Astrophysical Journal, 753, 106

[MacLeod et al., 2016]
MacLeod, C. L., Ross, N. P., Lawrence, A., et al. 2016, , 457,
389, 10.1093/mnras/stv2997

[MacLeod et al., 2019]
MacLeod, C. L., Green, P. J., Anderson, S. F., et al. 2019, , 874,
8, 10.3847/1538-4357/ab05e2

[Marchese et al., 2012]
Marchese, E., Braito, V., Della Ceca, R., Caccianiga, A., & Severgnini, P.
2012, Monthly Notices of the Royal Astronomical Society, 421, 1803,

[Martini & Schneider, 2003]
Martini, P., & Schneider, D. P. 2003, , 597, L109,

[Marziani et al., 2018]
Marziani, P., Dultzin, D., Sulentic, J. W., et al. 2018, Frontiers in
Astronomy and Space Sciences, 5, 6, 10.3389/fspas.2018.00006

[Masci et al., 2019]
Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, , 131,
018003, 10.1088/1538-3873/aae8ac

[McGreer et al., 2018]
McGreer, I. D., Fan, X., Jiang, L., & Cai, Z. 2018, , 155, 131,

[McGreer et al., 2013]
McGreer, I. D., Jiang, L., Fan, X., et al. 2013, , 768, 105,

[McHardy et al., 2018]
McHardy, I. M., Connolly, S. D., Horne, K., et al. 2018, , 480,
2881, 10.1093/mnras/sty1983

[McLure & Dunlop, 2004]
McLure, R. J., & Dunlop, J. S. 2004, , 352, 1390,

[McLure & Jarvis, 2002]
McLure, R. J., & Jarvis, M. J. 2002, , 337, 109,

[Morganson et al., 2014]
Morganson, E., Burgett, W. S., Chambers, K. C., et al. 2014, , 784,
92, 10.1088/0004-637X/784/2/92

[Mushotzky et al., 2011]
Mushotzky, R. F., Edelson, R., Baumgartner, W., & Gandhi, P. 2011, The
Astrophysical Journal, 743, L12, 10.1088/2041-8205/743/1/l12

[Netzer, 2013]
Netzer, H. 2013, The Physics and Evolution of Active Galactic Nuclei
(Cambridge University Press)

[Noda & Done, 2018]
Noda, H., & Done, C. 2018, , 480, 3898,

[Oke & Sandage, 1968]
Oke, J. B., & Sandage, A. 1968, , 154, 21, 10.1086/149737

[Palanque-Delabrouille et al., 2013]
Palanque-Delabrouille, N., Magneville, C., Yèche, C., et al. 2013,
, 551, A29, 10.1051/0004-6361/201220379

[Panda et al., 2019]
Panda, S., Czerny, B., Done, C., & Kubota, A. 2019a,
, 875, 133, 10.3847/1538-4357/ab11cb

[Panda et al., 2019]
Panda, S., Marziani, P., & Czerny, B. 2019b, , 882, 79,

[Pâris et al., 2017]
Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, , 597, A79,

[Pâris et al., 2018]
Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, , 613,
A51, 10.1051/0004-6361/201732445

[Rakshit & Stalin, 2017]
Rakshit, S., & Stalin, C. S. 2017, , 842, 96,

[Rasmussen & Williams, 2006]
Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for
Machine Learning (Cambridge, MA: MIT Press)

[Rau et al., 2009]
Rau, A., Kulkarni, S. R., Law, N. M., et al. 2009, , 121, 1334,

[Reichert et al., 1994]
Reichert, G. A., Rodriguez-Pascual, P. M., Alloin, D., et al. 1994,
, 425, 582, 10.1086/174007

[Richards et al., 2006]
Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, , 131, 2766,

[Richards et al., 2015]
Richards, G. T., Myers, A. D., Peters, C. M., et al. 2015, , 219,
39, 10.1088/0067-0049/219/2/39

[Ricker et al., 2014]
Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in ,
Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and
Millimeter Wave, 914320

[Riess et al., 2019]
Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D.
2019, , 876, 85, 10.3847/1538-4357/ab1422

[Risaliti et al., 2009]
Risaliti, G., Miniutti, G., Elvis, M., et al. 2009, The Astrophysical
Journal, 696, 160, 10.1088/0004-637x/696/1/160

[Rosen et al., 2016]
Rosen, S. R., Webb, N. A., Watson, M. G., et al. 2016, , 590, A1,

[Ross et al., 2013]
Ross, N. P., McGreer, I. D., White, M., et al. 2013, , 773,

[Ross et al., 2018]
Ross, N. P., Ford, K. E. S., Graham, M., et al. 2018, , 480,
4468, 10.1093/mnras/sty2002

[Ruan et al., 2014]
Ruan, J. J., Anderson, S. F., Dexter, J., & Agol, E. 2014, The Astrophysical
Journal, 783, 105, 10.1088/0004-637x/783/2/105

[Ruan et al., 2019]
Ruan, J. J., Anderson, S. F., Eracleous, M., et al. 2019, , 883,
76, 10.3847/1538-4357/ab3c1a

[Ruan et al., 2016]
Ruan, J. J., Anderson, S. F., Cales, S. L., et al. 2016, , 826,
188, 10.3847/0004-637X/826/2/188

[Rumbaugh et al., 2018]
Rumbaugh, N., Shen, Y., Morganson, E., et al. 2018, , 854, 160,

[Runnoe et al., 2016]
Runnoe, J. C., Cales, S., Ruan, J. J., et al. 2016, , 455, 1691,

[Rybicki & Press, 1992]
Rybicki, G. B., & Press, W. H. 1992, , 398, 169,

[Sánchez-Sáez et al., 2018]
Sánchez-Sáez, P., Lira, P., Mejía-Restrepo, J., et al.
2018, , 864, 87, 10.3847/1538-4357/aad7f9

[Scaringi et al., 2015]
Scaringi, S., Maccarone, T. J., Kording, E., et al. 2015, Science
Advances, 1, e1500686, 10.1126/sciadv.1500686

[Schawinski et al., 2015]
Schawinski, K., Koss, M., Berney, S., & Sartori, L. F. 2015, Monthly Notices
of the Royal Astronomical Society, 451, 2517, 10.1093/mnras/stv1136

[Schmidt et al., 2010]
Schmidt, K. B., Marshall, P. J., Rix, H.-W., et al. 2010, , 714,
1194, 10.1088/0004-637X/714/2/1194

[Schneider et al., 2007]
Schneider, D. P., Hall, P. B., Richards, G. T., et al. 2007, , 134,
102, 10.1086/518474

[Schneider et al., 2008]
—. 2008, VizieR Online Data Catalog, 7252

[Schneider et al., 2010]
Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, , 139,
2360, 10.1088/0004-6256/139/6/2360

[Sesar et al., 2007]
Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, , 134, 2236

[Shakura & Sunyaev, 1973]
Shakura, N. I., & Sunyaev, R. A. 1973, , 24, 337

[Shen, 2013]
Shen, Y. 2013, Bulletin of the Astronomical Society of India, 41, 61.

[Shen et al., 2008]
Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T., &
Schneider, D. P. 2008, , 680, 169, 10.1086/587475

[Shen & Ho, 2014]
Shen, Y., & Ho, L. C. 2014, , 513, 210, 10.1038/nature13712

[Shen et al., 2011]
Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, , 194, 45,

[Shen et al., 2019]
Shen, Y., Hall, P. B., Horne, K., et al. 2019, , 241, 34,

[Sheng et al., 2017]
Sheng, Z., Wang, T., Jiang, N., et al. 2017, The Astrophysical Journal, 846,
L7, 10.3847/2041-8213/aa85de

[Sheng et al., 2020]
Sheng, Z., Wang, T., Jiang, N., et al. 2020, , 889, 46,

[Simm et al., 2016]
Simm, T., Salvato, M., Saglia, R., et al. 2016, , 585, A129,

[Smith et al., 2018]
Smith, K. L., Mushotzky, R. F., Boyd, P. T., et al. 2018, , 857,
141, 10.3847/1538-4357/aab88d

[Śniegowska & Czerny, 2019]
Śniegowska, M., & Czerny, B. 2019, arXiv e-prints, arXiv:1904.06767.

[Stern et al., 2017]
Stern, D., Graham, M. J., Arav, N., et al. 2017, , 839, 106,

[Stern et al., 2018]
Stern, D., McKernan, B., Graham, M. J., et al. 2018, , 864, 27,

[Suberlak et al., 2017]
Suberlak, K., Ivezić, Ž., MacLeod, C. L., Graham, M., &
Sesar, B. 2017, , 472, 4870, 10.1093/mnras/stx2310

[Sun & Shen, 2015]
Sun, J., & Shen, Y. 2015, , 804, L15,

[Sun et al., 2018]
Sun, M., Xue, Y., Wang, J., Cai, Z., & Guo, H. 2018, , 866, 74,

[Trakhtenbrot et al., 2019]
Trakhtenbrot, B., Arcavi, I., Ricci, C., et al. 2019, Nature Astronomy,
3, 242, 10.1038/s41550-018-0661-3

[Vanden Berk et al., 2001]
Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, , 122,
549, 10.1086/321167

[Vanden Berk et al., 2004]
Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ,
601, 692, 10.1086/380563

[Vestergaard, 2002]
Vestergaard, M. 2002, , 571, 733, 10.1086/340045

[Vestergaard & Peterson, 2006]
Vestergaard, M., & Peterson, B. M. 2006, , 641, 689,

[Wang et al., 1996]
Wang, T., Brinkmann, W., & Bergeron, J. 1996, , 309, 81

[Wilhite et al., 2008]
Wilhite, B. C., Brunner, R. J., Grier, C. J., Schneider, D. P., &
vanden Berk, D. E. 2008, , 383, 1232,

[Wisotzki, 2000]
Wisotzki, L. 2000, , 353, 861

[Yang et al., 2017]
Yang, J., Fan, X., Wu, X.-B., et al. 2017, , 153, 184,

[Yang et al., 2018]
Yang, Q., Wu, X.-B., Fan, X., et al. 2018, The Astrophysical Journal, 862,
109, 10.3847/1538-4357/aaca3a

[Yang et al., 2020]
Yang, Q., Shen, Y., Chen, Y.-C., et al. 2020, , 493, 5773,

[Zhu et al., 2018]
Zhu, F.-F., Wang, J.-X., Cai, Z.-Y., et al. 2018, The Astrophysical Journal,
860, 29, 10.3847/1538-4357/aac2d7

[Zu et al., 2013]
Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ,
765, 106, 10.1088/0004-637X/765/2/106

[Zu et al., 2011]
Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, , 735, 80

[Zuo et al., 2012]
Zuo, W., Wu, X.-B., Liu, Y.-Q., & Jiao, C.-L. 2012, The Astrophysical Journal,
758, 104, 10.1088/0004-637x/758/2/104