Universality in the random walk structure function of luminous quasi-stellar objects

Rapidly growing black holes are surrounded by accretion disks that make them the brightest objects in the Universe. Their brightness is known to be variable, but the causes of this are not implied by simple disk models and still debated. Due to the small size of accretion disks and their great distance, there are no resolved images addressing the puzzle. In this work, we study the dependence of their variability on luminosity, wavelength and orbital/thermal timescale. We use over 5,000 of the most luminous such objects with light curves of almost nightly cadence from >5 years of observations by the NASA/ATLAS project, which provides 2 billion magnitude pairs for a structure function analysis. When time is expressed in units of orbital or thermal timescale in thin-disk models, we find a universal structure function, independent of luminosity and wavelength, supporting the model of magneto-rotational instabilities as a main cause. Over a >1 dex range in time, the fractional variability amplitude follows log(A/A0)≃1/2×log(Δt/tth)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log (A/{A}_{0})\simeq 1/2\times \log (\Delta t/{t}_{{{{\rm{th}}}}})$$\end{document}. Deviations from the universality may hold clues as to the structure and orientation of disks. A study of the emission variability of roughly 5,000 of the brightest quasars supports the presence of an optically thick yet geometrically thin accretion disk, and may provide a means of measuring the size and inclination of the disk.

dissipates orbital energy.It is now clear that Keplerian disks with weak magnetic fields are prone to magneto-rotational instability (MRI) 5 .MRI is a fundamental instability converting differential rotation into turbulent fluid motions, feeding an energy cascade, and explaining the anomalously high viscosity needed for luminous QSO disks 6 .
Accretion disks are an elegant explanation for the high luminosity and compactness of QSOs.The standard model of optically thick and geometrically thin disks predicts their temperature profiles and emission 7 .As they are too small to be spatially resolved by current facilities, we seek clues about their properties from the variability of their emission, which is ubiquitous on all timescales.While MRI must be one source of short-term fluctuations, another must be the rapidly fluctuating heating from the X-ray corona near the black hole [8][9][10] .The "lamppost" model exploits the latter and relates the sizes of accretion disks to time lags between light curves of different passbands probing regions of different temperature and radius in the disk [9][10][11][12][13][14][15] .
Monitoring of QSO samples has revealed that fluctuations can be described by a 'damped random walk' model [16][17][18] : brightness changes tend to be larger when the time delay between any two measurements is longer.However, this model provides only a behavioural analogy and no physical insight.Mathematically, the light curve of an object can be statistically described by a variability "structure function" (SF): the typical brightness difference is expressed as a function of time interval 19,20 .The SFs of QSOs are characteristically different from other types of variable objects, so that QSOs can even be identified in a time-domain survey by the SF alone 21 .
Measurements of the SF agree on two trends: fluctuations on a fixed timescale get stronger in less luminous QSOs and in bluer passbands 17,[22][23][24][25][26] .Thus, the SF can be parametrised as where amplitudes A of fractional flux changes are in terms of magnitude, and L is the measured luminosity of the QSO.Both the time separation between the two measurements, ∆t, and the wavelength of observed light, λ, are expressed in the QSO restframe with cosmological redshift and time dilation removed.Outside of the randomwalk regime, the amplitudes may be suppressed at shortest timescales [27][28][29][30][31] for reasons that are not yet established.They tend to level off above a damping or decorrelation timescale 17,31,32 , which has been shown to be driven by the thermal timescale of the accretion disk 32 .Past works agree that B L < 0 and B λ < 0 but disagree on the strength of the trends.Some complicating factor may not be controlled for in existing studies.The SF slope γ = d log A/d log ∆t at fixed (L, λ) is in the range 20,22,23,25 of 0.2 to 0.6.B L = d log A/d log L at fixed (∆t, λ) describes the dependence on disk luminosity (hence size) and is reported in the range 17,22,23,25,26 of −0.33 to −0.16.Finally, B λ = d log A/d log λ at fixed (∆t, L) describes the trend with rest-frame wavelength and thus location of the fluctuating source within the disk; measurements range 17,23,25 from −0.75 to −0.48, although observations of non-linear relations also include slope values outside this range.
Here, we present a new analysis of the variability structure function of the ∼ 5, 000 brightest QSOs in the sky using light curves with unprecedented sampling.The sample has a median black-hole mass of ∼ 2 billion solar masses according to an analysis of spectra from the Sloan Digital Sky Survey (SDSS) 33 .The light curves are from the 0.5-metre NASA/ATLAS telescope (Asteroid Terrestrial-impact Last Alert System 34 ), which was built to find dangerous asteroids that might threaten the Earth.ATLAS is designed to image the entire sky seen from its location in Hawaii four times every night, weather permitting.With up to 1,200 nights per QSO per band, we get ∼ 2.1 billion magnitude pairs to constrain the QSO SF.This is one of the best data sets for studying QSO variability before the new Vera C. Rubin Telescope operates.

Results
The new light curves include observations from 2015 to 2021 in two passbands, cyan and orange (centred on 533 and 679 nm, respectively).We limit the sample to redshift 0.5 < z < 3.5 (z < 2.4 in the cyan band), where QSOs are so luminous that their host galaxies will not affect the measurements.We reject radio-loud QSOs, where variability may be enhanced by shocks in jets and mask the pure accretion disk behaviour 35 .We also reject gravitationally lensed QSOs and QSOs with close neighbours as seen by the ESA Gaia mission 36 to avoid blended light curves.
Previous work typically fitted a global model for the SF, such as Equation 1, to a whole data set.The richness of this new data set offers the alternative to split the data into subsets without incurring the penalty of noisy statistics.We first split our data into twenty bins in ∆t and fit the trends with luminosity and wavelength separately for each bin, using the equation Thus, we get twenty different estimates for B L and B λ , independently for each ∆t bin.
As we can now study the trends on different timescales, we may see different regimes of behaviour.We use the measured variability on timescales of 1 day to estimate the noise in the observations.At the long end, window effects from the finite length of the light curves may cause noise in the parameters 37 .We restrict the wavelength range to longer than 130 nm (avoiding flux from the strong Ly-α emission line) and shorter than 300 nm (avoiding a previously seen upturn in variability 17,22,24 , which is currently unexplained and discussed further below).We use two types of statistical analysis, simple bin averages, and a computationally more intense bootstrap analysis.
Figure 1 shows the resulting parameters for a range in ∆t from 10 to 337 days.Both methods of analysis produce consistent estimates.The estimates for B L and B λ show a strong trend with ∆t.It is immediately clear that a global fit will produce estimates for B L and B λ that depend on the fitting range in ∆t and on the distribution of pairs in the range, and thus on subtleties of survey cadence.We suggest this fact explains the variety of estimates obtained in the past.We demonstrate how results depend on the range of behaviour forced into a global fit in Figure 2, where fiducial global estimates for B L and B λ are obtained by choosing different global fitting ranges for ∆t.
Finally, we improve the global fit by restricting it to the random-walk regime.Given that our sample comprises the most luminous known QSOs, we expect damping timescales well in excess of a year, which do not affect our analysis.Having controlled for ∆t by fitting the trends in each bin separately (Figure 1), we reveal two regimes in behaviour: at short ∆t both B L and B λ drift with ∆t, whereas at longer ∆t they both remain constant.Thus, at longer ∆t the choice of fitting range will have little impact.Here, we find fit parameters of γ ≈ 1/2, B L ≈ −1/4 and B λ ≈ −1.Past estimates of B L cover values from a lower range 17,26 of −0.33 to −0.29, depending on the chosen parametrisation, to an upper limit 23 of −0.16 and include values very close 22,25 to −1/4.

Changing the clock for a universal structure function
So far, B L and B λ have described how variability amplitudes behave at fixed ∆t.We now look at the timescales, on which disks of different L, observed at different λ, show a fixed amplitude A and thus potentially a universal structure function if their slopes γ are the same as well.We find these timescales to be parametrised as: The first study in this direction empirically found C L = 0.47 ± 0.06 and argued that the fractional optical variability per disk orbital or thermal timescale is likely constant 38 .This argument is based on the well-known fact that the predicted temperature profile in thin disks 7 , T (R) ∝ R −3/4 , sets a mean radius for the origin of radiation and thus a mean orbital timescale t orb of the emitting surface, given λ and a size, and thus L, of the disk 17 .Note, that for a fixed viscosity of the disk material, the thermal timescale t th is proportional to the orbital timescale and we get 38 : One later study 24 reported C L ≃ 0.41 ± 0.04.Our result shows that the concept of a constant fractional variability per orbital or thermal timescale 38 is not only consistent with measurements of C L but also our C λ .Thus, we rephrase the whole SF from Equation 1 with time in new units, choosing those of thermal timescale given its role for the damping timescale 32 : This form is expected for fluctuations arising from MRI, where the equations of motion scale with the orbital timescale (and thus the thermal timescale) 6,39 .It also applies in the corona-heated accretion-disk reprocessing (CHAR) model 40 , which is motivated by the concern that a lamppost will not create sufficient amplitude by heating alone.
The CHAR model proposes that the corona couples directly to the magnetic field in the disk and drives stronger fluctuations than expected from a lamppost corona alone.We determine t th for each QSO and passband given its (L, λ). Figure 3 shows the resulting SF expressed both in natural time and with a t th -scaled clock.Here, we increased the time resolution but reduced noise by grouping QSOs into subsamples split by (L, λ).The wide (L, λ)-range of our sample causes a wide range of amplitudes at fixed ∆t (left panel), but at fixed ∆t/t th (right panel) we find indeed a tight range, supporting the concept of a universal structure function with a constant slope and intercept once expressed in suitable units of time.The relation holds over more than 1 dex in ∆t/t th , limited by the length of our light curves at long timescales.At short timescales, the amplitude suppression relative to the random walk remains.The CHAR model actually predicts a suppression with a break timescale proportional to the thermal timescale, which would imply a fixed break timescale in Figure 3, right panel.However, we see a spread of timescales such that the break depends more strongly on luminosity than the thermal timescale, which will be investigated in future work.For our final global fit of the SF, we avoid the suppressed regime by using only data at log A > −1.4.We find a bootstrap slope with natural clock time of γ = 0.503 ± 0.001 and with time in units of thermal timescale of γ th = 0.510 ± 0.002, extremely consistent with the γ = 1/2 of a random walk (Table 1 summarises all fit parameters).

Outlook
The transition from a regular clock to an orbital or a thermal time-scaled clock has standardised the random walk of QSOs as the wide dispersion observed among objects with a variety of luminosities and at a variety of restframe wavelengths has disappeared.Instead, the SF of QSOs appears now as an intrinsically tight relation.
A universal SF opens a window to characterise QSOs in more detail by deviations from the mean relation using future higher-precision measurements.Such deviations are expected to arise from any effect that varies the orbital and thermal timescale estimates.These effects include the dependence of the observed luminosity on disk orientation and black-hole spin 41 , which are not independently known for individual objects; the true temperature profile of the accretion disk, which is still debated as it may change from T ∝ R −3/4 in an inner part driven by viscosity to T ∝ R −1/2 in an outer part dominated by irradiation from an X-ray corona acting as a lamppost 39 ; and the actual interior structure of the disk, which may well change from thin to slim disks as the accretion rate approaches the Eddington rate 39,42,43 .
Deviations from a simple picture are generally evident at wavelengths of λ > 300 nm.Here, we measure an upturn in A, relative to a trend of B λ ≈ −1, which has been seen before 17,22,24,44 (see also Figure 4, right panel).It suggests that fluctuations in the outer, cooler parts of QSO accretion disks are to be enhanced by an additional mechanism, which has not yet been explained.Also, we excluded λ < 130 nm, as we cannot disentangle light from the disk continuum and light from the broad Ly α line.The hotter inner disks are harder to study as a wide swath of spectrum in the UV domain is absorbed by intervening gas along the line-of-sight to QSOs.
The SF suppression at short ∆t has been observed across a range in luminosity: observations of low-luminosity AGN with the Kepler space mission have revealed breaks in the SF on time scales of days to weeks [27][28][29] , although it has also been argued that there "remain significant levels of spacecraft-induced effects in the standard pipeline reduction of the Kepler data" 45 .Steeper slope changes have been seen in a recent study of luminous QSOs on time scales of 10 to 15 days 31 , while our break time scales range from ∼ 15 to ∼ 50 days.The reasons behind these breaks are still under debate: either the structure function is intrinsically different such as in the CHAR model 40 , or a "smearing" effect may be caused by synchronised emission from different parts of the disk reaching the observer after different time lags 30 .
The forthcoming 10-year Legacy Survey in Space and Time (LSST [46][47][48] ) planned at the Vera C. Rubin Telescope for the years 2024-34 will eclipse this data set with better data: more passbands will probe a wider temperature range in the accretion disk.Deeper photometry will probe QSOs at much fainter luminosity.This, and the reduced photometric errors in LSST will provide low-noise structure functions for individual objects and reveal how tight and universal the structure function of QSO variability really is.LSST should also be a game changer for analysing subtleties in the shortterm breaks in the SF.

Conclusion
We have determined with unprecedented precision the parameters (γ, B L , B λ ) of the variability structure function for the UV continuum emission of luminous QSOs.It is consistent with a universal random walk once time is expressed in units of the orbital or thermal timescales of the emitting material.The mean amplitudes of fluctuations can be described by the simple law log(A/A 0 ) = 1/2 × log(∆t/t th ), where a single fit parameter A 0 represents the overall QSO population.Note, that the thermal timescale has recently been shown to drive the damping timescale of QSO variability as well 32 , although this fact may not have affected the random walk at shorter timescales.
The result appears consistent with the temperature profiles of thin-disk models 7 , at least at wavelengths from 130 to 300 nm.Deviations observed at longer wavelengths that show the cooler, outer parts of disks 17,22,24,44 still need to be explained.It also appears consistent with an origin of the fluctuations in magneto-rotational instabilities (MRI) in differentially rotating accretion disks that are made of magnetised plasma; there, the equations of motion scale with the orbital timescale 6,39 .MRI is suggested to produce variability with a random walk phenomenology 49 and γ ≈ 0.4 to 0.5.Disks could thus be ordered in a 1-parameter family with the mean thermal timescale at a chosen wavelength as an ordering parameter.The structure function of QSOs may then be a universal relation apart from the short-term suppression.It appears intrinsically tight and would be broadened in observation by disk orientation and black hole spin, as well as host galaxy dust extinction affecting the inferred luminosity.It would also be broadened by intrinsic variety in vertical disk structure driven by accretion rates.
The role of the lamppost for driving QSO variability is likely confined to a contribution of modest amplitude on shorter timescales.If the CHAR model also described the physics of disk variability correctly, then it would be hard to disentangle genuine lamppost reverberation from outwards-propagating magnetic fields; and if the velocity of the latter dropped below the speed of light 40 , the sizes of disks inferred by a reverberation analysis could be overestimated.
Future data from the 10-year LSST project will shed light on these questions and more subtle properties with its extremely precise measurements for a nearly all-sky QSO sample of unprecedented size, depth, variety and completeness.

QSO sample and data cleaning
We combined the Million Quasars Catalog 50 (MILLIQUAS v7.1 update) with the Gaia eDR3 catalogue 36 and select 6,163 spectroscopically confirmed, non-lensed QSOs with Gaia magnitude Rp < 17.5, redshift 0.5 < z < 3.5, declination δ > −45 • , Galactic foreground reddening 51 of E(B − V ) SFD < 0.15 and Gaia BpRp Excess Factor < 1.3 (indicative of single, unblended, point sources).We also required that the MILLIQUAS and Gaia positions are coincident within 0.3 ′′ , which excludes blended objects of similar SED such as the multiple images of lensed QSOs, and that sources have no Gaia neighbours within 15 ′′ to keep the ATLAS photometry unaffected.We confirmed with recent lists of lensed QSOs 52,53 that none of these are included in our sample.Our sample has 6,115 sources in common with NASA/ATLAS.The ATLAS photometric system is consistent with other reliable sources of photometry such as Pan-STARRS and Gaia 54 .
The selection is based on Gaia eDR3 mean photometry, which is arguably the bestcalibrated among large astronomical datasets 55 ; the Gaia observations extend from mid-2014 to mid-2017 and overlap with the period of the ATLAS lightcurves (LCs).Thus, there will be no noticeable selection bias in favour of QSOs that might be brighter during an early selection epoch but fainter when the LCs are observed 56 .
We further remove six QSOs from the sample, which have large flux errors throughout their LCs (90-percentile flux errors in the orange band of σ fν ,o > 150µJy or in the cyan band LC of σ fν ,c > 85µJy).On individual images, the faintest objects in the sample have signal-to-noise ratios (SNR) of > 11 in the orange band and > 13 in the cyan band.Up to four images are available per band and night.
We reject radio-loud QSOs by cross-matching the sample with catalogues from the Faint Images of the Radio Sky at Twenty-cm 57 (FIRST) Survey, the NRAO VLA Sky Survey 58 (NVSS), and the Sydney University Molonglo Sky Survey 59 (SUMSS).We use matching radii between the MILLIQUAS and FIRST, NVSS, and SUMSS coordinates of 3", 12", and 11", respectively.1,242 QSOs are matched with at least one radio catalogue.If a QSO is matched to multiple radio sources in NVSS or SUMSS, the radio detection with the closest separation is used.We apply a radio loudness criterion of 0.4 × (m o − t NVSS/SUMSS ) > 1.3 for matches with NVSS and SUMSS.Radio sources below this threshold are labelled radio-intermediate.QSOs with matches in FIRST but not in NVSS turned out to be all radio-intermediate.The sample contains 775 radioloud and 467 radio-intermediate objects, and 4,848 radio non-detections.Thus, our final sample contains 5,315 QSOs, of which 3,607 have single-epoch virial black hole mass measurements derived from SDSS spectra 33 .
We clean the ATLAS LCs first by excluding all observations with large errors of log(σ fν ,o ) > −3.94−0.12×(mo −16.5) and log(σ fν ,c ) > −4.17−0.10×(mc −16.5), which constitute about 5 percent of all observations; and further by comparing each measurement with other observations within ±7 days and removing outliers with a 2σclipping to reject spurious variability, e.g., due to temporary chance blending with faint asteroids.If there is only one observation within ±7 days, it is retained.

Luminosity
We estimate QSO luminosities at restframe 3,000 Å (L 3000 ) using a linear interpolation of the Gaia Bp and Rp magnitudes.At redshifts of 0.5 < z < 0.7 and z > 1.56, the restframe 3,000 Å point is extrapolated beyond the pivotal wavelengths of the Gaia filters.We correct for foreground dust extinction in the Milky Way 51 using bandpass coefficients 60 of R Bp = 3.378 and R Rp = 2.035.Luminosity distances are derived in a flat ΛCDM cosmology with Ω m = 0.3 and a Hubble-Lemaître constant of H 0 = 70 km sec −1 Mpc −1 .Luminosity errors from this procedure are expected to be < 20% and thus much smaller than those introduced by host galaxy dust and the orientation of the disk.Luminosity estimates from broadband photometry does not distinguish between contributions from the accretion disk continuum and those from line emission, especially from the Fe complex that is strong in the vicinity of 3,000 Å .Comparing our L 3000 estimates with those obtained by spectral fitting on the available subsample 33 , we find an average offset of 0.07 dex and reduce our L 3000 values to correct for line contributions.Bolometric luminosities, as used by some authors 24 , are derived with log(L bol /L 3000 ) = log(f BC × λ) = log(5.15× 3000 Å) = 4.189, using f BC = 5.15 as the bolometric correction factor 61 .

Thermal timescale
For each QSO and passband, we estimate the characteristic (defined in the following) thermal timescale of the emitting material in the accretion disk.We combine the equation for t th 39 and an equation for the scale length of the disk given L, λ values 62 , and get: where h p is the Planck constant, α = 0.1 63 is the assumed viscosity parameter, and i = 45 • is the mean inclination.We assume a radiative efficiency of η = L bol /( Ṁ c 2 ) = 0.1, where Ṁ is the mass accretion rate; changes in η due to variations in black-hole spin affect the inner cutoff of the accretion disk.For a given Ṁ , such a variation has little effect on the outer disk, the overall L 3000 , and the thermal timescales at the wavelengths observed here.Instead, it mostly affects the emission of higher-energy photons and L bol .

Variability structure function
We adopt the following definition of the variability structure function (SF) 64 : where m i and m j are any two observations and σ is the magnitude error due to noise.We find that the magnitude errors quoted by ATLAS underestimate the noise and instead derive σ 2 from the apparent intra-day variability, which we expect to be negligible for radio-quiet QSOs.We determine < σ 2 > as a function of observed magnitude, m obs , such that the SF is anchored at A ≈ 0 for ∆t < 1 day.The required noise model is shown in Figure 5 and given by with (n 0 , n 1 ) = (−12.411,0.585) and (−12.428,0.573) for the orange and cyan passband, respectively (pure Poisson noise suggests n 1 = 0.4).
Informed by the pair sampling distribution of the LCs, we reduce window effects from the finite extent of the LCs by ignoring pairs with ∆t obs ≥ 0.47max(∆t obs ).This cut is applied in the observed frame, while the rest of the paper uses the restframe ∆t = ∆t obs /(1 + z), corrected for cosmological time dilation.

Binning the time axis
We bin the time axis in two different ways depending on purpose.A first case is used for the relation log A(∆t) in Equation 2and Figures 1, 2 and 4. A second case is used for fits over the final fitting range and for log A(∆t/t th ) in Equation 6, Figure 3 and Table 1.
In the first case, we split the ∆t axis into 20 bins.The first ∆t bin contains all pairs with ∆t < 1 day, which is used only for noise analysis.Further bins are chosen to balance the number of pairs among the bins.The average number of pairs per ∆t bin are ∼103 million and ∼8 million for the orange and cyan passband, respectively.We split the z axis into 9 bins, thus grouping QSOs with observations of similar restframe wavelengths.In the cyan passband, we exclude QSOs at 2.4 < z < 3.5 to avoid contamination from the Ly-α line.In the orange passband, we use this z range as a single, highest-z, bin.The remaining eight bins cover the range of 0.5 < z < 2.4 with roughly balanced QSO numbers.In each z bin, QSOs are further split into ten L 3000 bins, each with almost equal QSOs numbers.In total, there are 1,800 bins in (z, L 3000 , ∆t).In each bin, the mean value of z and L 3000 and the mid-point of ∆t is used for the analysis and figures.
In the second case, we split the ∆t axis into 100 bins.The shortest ∆t < 1 day bin is exactly the same as in the first binning method, and further bins are balancing the number of pairs among the bins for both passbands.The QSOs are grouped into 16 bins along the t th axis with an equal number of QSOs per bin, which is 2,950 and 4,134 in orange and cyan, respectively.In the ∆t/t th axis, we split the observed pairs into 100 bins of constant number, independently for each t th group in each passband after excluding pairs with ∆t < 1 and ∆t > 337 days.Combining both passbands, there are in total 3200 bins for (t th , ∆t) and separately 3200 bins for (t th , ∆t/t th ).In each bin, the mid-point of ∆t and ∆t/t th is used for the analysis and figures.

Simple bin average structure function
For each individual QSO, we calculate its SF with Equation 8 from all available pairs in each ∆t or ∆t/t th bin.We use a 3σ-clipped mean of magnitude differences as < ∆m >.Then, for each (z, L 3000 , ∆t) or (t th , ∆t) or (t th , ∆t/t th ) bin we calculate the median squared amplitude, A 2 , and the standard error of the median (SEMED), which ensures that QSOs with negative noise-corrected amplitudes A 2 , due to an oversubtraction by the σ 2 noise model, are still included in the statistics.

Bootstrapping the structure function
Using bootstrapping, we can account better for the contribution of each QSO.Here, we select random pairs in each of the (z, L 3000 , ∆t) or (t th , ∆t) or (t th , ∆t/t th ) bins.We randomly pick a QSO N times from all the QSOs in that bin, where N is the total number of QSOs in that bin.For each picked QSO, we randomly select 300 distinct pairs in the orange and 30 distinct pairs in the cyan passband, and then calculate the SF using Equation 8.This process is repeated 100 times to get the bootstrap distribution, from which we calculate the mean and standard deviation for each bin.

SF fitting
We fit a structure function of the form log(A(∆t)) = B 0,∆t + B L,∆t log( L 3000 10 42.7 erg/s/ Å) + B λ,∆t log( using a Levenberg-Marquardt least-squares fit 65 to the (z, L 3000 ) bins for each ∆t bin separately.Here, log A and its error may be results of bin averages or from the bootstrap method.λ rf is the pivotal wavelength of a passband divided by 1 + z.The normalisation factors for λ rf and L 3000 are chosen to be close to the median of the sample.Data with λ rf > 3000 Å are excluded from the fit as it deviates from a linear relation.In Figure 4, we show an example for one ∆t bin: the left and middle panels show the two passbands.The right panel shows the intercepts of the fits at log(L 3000 )− 42.7 for different λ rf .We also compare the upturn at longer wavelengths with results from the literature 24 .The B 0,∆t , B L,∆t , and B λ,∆t from each ∆t bin are shown in Figure 1.

Global parameter estimates using broad ∆t ranges
We fit a global model of the form to the SF data in the (z, L 3000 , ∆t) bins with λ rf < 3000 Å Ẇe test the robustness of global parameter averages by varying the time span of data that is included in the global fit; we vary the lower edge ∆t 0 , while leaving the upper edge at ∆t=337 days unchanged.Figure 2 shows the parameter estimates as f (∆t 0 ).
The structure function as f (∆t) and f (∆t/t th ) For Figure 3 we split the sample for each passband separately into 16 groups ranked by t th .We then fit Equation 11 to the points with log A > −1.4 using the errors for weighting (see Table 1 and left panel of Figure 3 for results).Similarly, we fit the bins with log A > −1.4 to Equation 6(see Table 1 and right panel of Figure 3 for results).
For reference, we find log A 0 = −1.185when fixing γ = 1/2 in the bootstrap method.Note an upwards deviation at λ > 300 nm, which is consistent with the dotted line from one past work 24 and also seen in other works 22,44 .
Figure 5: The variability amplitude, π 2 < ∆m > 2 , at ∆t < 1 day as a function of observed magnitude.Gray symbols are individual QSOs.Dotted lines show the 3σclipped mean within fine magnitude bins.Solid lines show the fit of Equation 9.

Figure 1 :
Figure 1: Coefficients of the structure function from Equation 10 in narrow time intervals, comparing bootstrap results and simple bin averages.This figure contains the result of 5315 QSOs.The bootstrap error bars are the standard deviation of the distribution while the error bars of the simple bin averages are the standard error of the median (SEMED).Data points and error bars are based on ∼ 2 billion magnitude pairs from 5315 QSOs.

Figure 2 :
Figure 2: Global coefficient estimates from Equation11(bootstrap case).∆t 0 is the short end of the fitting window, while the long end of the window (∆t end ) is fixed at 337 days.The crosses and their error bars are the means and the standard deviations of the bootstrap distributions.When short timescales are included, the estimates drift, as non-constant behaviour is included into the fit.Restricting the window to the longest timescales, shrinks the data set and makes the fits noisier (data points and error bars are based on ∼ 2 billion magnitude pairs from 5315 QSOs).

Figure 3 :
Figure 3: Variability amplitudes log A vs. log ∆t (left) and vs. log(∆t/t th ) (right).The colour encodes the thermal timescale of the measurement from blue for the shortest to red for the longest t th .The best-fit lines to the data at log A > −1.4 have slopes of γ = 0.503 (left) and γ th = 0.510 (right), consistent with a random walk.Solid lines show the fit range while the dotted lines are extrapolated.At short ∆t, the variability appears suppressed and this effect reaches to longer ∆t in larger (here: redder) disks.

Figure 4 :
Figure 4: One example time separation bin, ∆t rest = [125; 141] days.The 5,315 individual QSOs are points, while large symbols and their error bars are bootstrap mean and standard deviation of the variability amplitudes in bins of luminosity L 3000 and redshift z, lines are fits to bootstrap mean values (left: orange passband; centre: cyan passband).Right: intercepts and errors of the fits at fixed L 3000 for different rest-frame wavelengths, where either passband is seen to sample the same underlying behaviour.Solid lines are fit ranges, dashed lines are extrapolations shown for clarity.Note an upwards deviation at λ > 300 nm, which is consistent with the dotted line from one past work24 and also seen in other works22,44 .

Table 1 :
Coefficients for final fit results.