The universal power spectrum of Quasars in optical wavelengths: Break timescale scales directly with both black hole mass and accretion rate

Explore Supermassive Black


Introduction
The flux variability of quasars was recognised early on as one of the few tools to study the structure and dynamics of their central engine and immediate environs, as reviewed by Cackett et al. (2021).Different studies of the continuum optical variability have concluded that at least part of the variations must originate in the accretion disc, as opposed to representing only a reprocessed response to the highly variable X-ray emission, both on long and short timescales (e.g.Krolik et al. 1991;Arévalo et al. 2008;Ai et al. 2010;Edelson et al. 2015;Lira et al. 2015;Smith et al. 2018).Therefore, optical variations are a viable tool to constrain the behaviour of the accretion disc.Several authors have made significant progress in correlating variability properties with basic quasars parameters, with increasingly large samples and longer or more frequently sampled light curves and more complex methods (Bauer et al. 2009;Caplar et al. 2017;Sánchez-Sáez et al. 2018;Li et al. 2018;Luo et al. 2020;Tachibana et al. 2020;Arévalo et al. 2023;Tang et al. 2023).In particular, the search for a characteristic timescale of variations has been hampered by the long lightcurves that would be required to measure it (e.g.Kozłowski 2017).However, Burke et al. (2021) recently measured a strong, positive correlation between the damping timescale and black hole mass, as obtained from the fit of Damped Random Walk models to high-quality, light curves of quasars.In this work, we aim to establish whether the characteristic timescale also depends on another fundamental quasar parameter, namely the accretion rate, quantified through the Eddington ratio R Edd .
The characteristic variability timescale, measured as a bend in the power spectrum of the lightcurves or as the de-correlation timescale of a Damped-Random-Walk (DRW) model (Kelly et al. 2009) can depend on the rest-frame wavelength of the lightcurves.In particular, longer characteristic timescales are predicted for longer wavelengths if the variations are produced on a local characteristic timescale (e.g. the radius-dependent thermal timescale of the accretion disc) and different wavelengths are preferentially emitted at different radii in the disc.Some observational results favour this scenario, for example, Sun et al. (2014) found that the amplitude of variability depends on the wavelength on short timescales but not on long timescales, which means that the power spectrum shape must depend on the emitted wavelength.Xin et al. (2020) find that optical and UV light curves are not perfectly correlated, using Galaxy Evolution Explorer and the Catalina Real-time Transient Survey light curves of over 1300 quasars, implying that the power spectra in different optical/UV bands are not simply scaled by amplitude.More recently, Stone et al. (2022) found that the characteristic timescale of variations depends on the emitted wavelength, from a sample of 190 quasars with decades-long light curves.Therefore, to refine the search for correlations between characteristic timescales and physical parameters, the light curves should necessarily probe the same rest-frame wavelength for all objects.
In this work we explore the optical power spectra of quasars in a single band and a narrow redshift range (see Sec. 2), to avoid potential dependencies of the power spectrum shape on the restframe wavelength of emission.We calculate the median power spectra for bins in mass, M, and R Edd , to establish the dependence of the power spectrum break frequency and normalisation on both parameters.The variance is estimated using the Mexican Hat power spectrum (Arévalo et al. 2012), which is well suited to light curves with gaps and uneven sampling (see Sec. 3) and the resulting power spectra are presented in Sec. 4. Model fitting to a bending power law model, with different slopes, is done by simulating light curves with specific model parameters and folding them through the same sampling and variance estimation process as the real data, as described in Sec. 5. Best-fitting models and scaling parameters were searched by minimising the difference between data and model power spectra.With this minimisation we show that the scaling relations of amplitude and characteristic timescale of the variability with mass and R Edd are largely independent of the model and select the better-fitting power spectrum model parameters in Sec. 6.A basic interpretation of the scaling relations in terms of accretion disc parameters is presented in Sec.7 and the conclusions are summarised in Sec. 8.

Data
The selection of quasars included in this work and the light curves used are the same as those presented in Arévalo et al. (2023).In summary, we use g-band light curves obtained by the Zwicky Transient Facility (ZTF, Masci et al. 2019) through their public data release DR14 of the 4770 quasars that: i) are in the redshift range 0.6<z<0.7;ii) have mass and R Edd estimated through Sloan Digital Sky Survey (SDSS) spectra by Rakshit et al. (2020) with reported statistical errors on mass less than 0.2 dex; iii) have a mean g-band ZTF magnitude less than 20, and that iv) their ZTF light curves, restricted to data obtained with a single CCD, excluding bad nights and multiple observations in the same night, have at least 90 data points and are at least 900 days long.Some objects have 2 light curves because they were observed with two different CCDs on different days and both light curves conform to the minimum standards.The work described below used the resulting 5433 valid light curves, which have a median length of 1547 days and a median number of points of 233.We note that since all galaxies are at z>0.6 they appear almost point-like in the ZTF images.Therefore the PSF photometry on the science images, which is used to produce the light curves of the data release, is sufficiently accurate.

Estimation of the Variance
The procedure to obtain the variance of the light curves follows the procedure described in Arévalo et al. (2023).Below we briefly summarise the method.
We estimated the variability power density on 4 wellseparated timescales with the Mexican hat filter (Arévalo et al. 2012) and multiplied them by the corresponding frequency to obtain the variance at the four timescales.Since the light curves are mean-subtracted and divided by the mean flux, these variances are dimensionless and represent the amplitude of variability relative to the mean flux, squared.The timescales chosen were 30, 75, 150 and 300 days in the rest-frame of the quasars.The additional variance produced by photometric errors of the light curve was subtracted from the measured variances of each light curve using simulated light curves of pure noise, constructed from the error bars of the real light curves.We further correct these noise estimates by carrying out the same procedure of variance estimation and noise subtraction for ZTF lightcurves of non-variable stars, which should result in net variances distributed around 0. We found small deviations from 0 for the median net variances, mostly towards negative values.We fitted the relation between median net variances and apparent magnitude with fourth-order polynomials and added this "missing variance" as a function of magnitude to the net variance of each quasar, for each timescale.These were always small corrections as can be seen in appendix B of Arévalo et al. (2023).
We split the sample of quasars on a grid of M and R Edd in order to measure the power spectrum dependence on these two basic parameters.The ranges in M and R Edd used for the grid are displayed in the line and column headings, respectively, in Table 1 and the table contains the number of light curves included in each bin.The analysis below will only use bins with 10 or more lightcurves, in order to obtain reliable estimates of the uncertainty of the variances, by bootstrapping the measured variances within each bin as described in Arévalo et al. (2023).

Power Spectrum of the data
Light curves can be characterised through the power spectrum, which quantifies the amount of variance found on different timescales or equivalently, different frequencies f = 1/t.Here we collect the variance estimates in four timescales to construct a low-resolution power spectrum for each object, covering one decade in frequency.We note that whenever we plot power spectra we will multiply the power density by the corresponding frequency so that the quantity plotted is Power × frequency.
We split the quasar sample in bins of both R Edd and mass as described in Sec.3 and plot the median power spectrum of each M−R Edd bin in Figs. 2 and 3.The error bars in the plots represent the standard deviation of the medians obtained for 1000 random re-samples of the data in each bin.We note that since all the quasars in the sample are at approximately the same redshift (0.6 ≤ z ≤ 0.7) and all the light curves were observed with the same filter, they sample the same rest-frame wavelength λ ∼ 2900Å.This allows us to determine the dependence of variability amplitude on M and R Edd without having to account for the dependence of amplitude on the emitted wavelength (e.g.Sun et al. 2014;Sánchez-Sáez et al. 2018;Stone et al. 2022).
Figure 1 compares the dependence of the power spectra on R Edd , for objects within a narrow range in mass (left), to the dependence of the power spectra on mass, for quasars within a narrow range in R Edd (right).Evidently, the power spectrum depends on both parameters but in different ways: the Eddington ratio af- fects mostly the normalisation and the black hole mass affects mostly the slope, in the frequency range sampled by the light curves.
Figure 2 highlights the dependence of the power spectra on R Edd .On each panel, the power spectra at three or more accretion rate levels are shown in different colours, and the different panels correspond to different mass bins.The most noticeable difference between the power spectra is the normalisation, with higher accretion rates having increasingly lower normalisation.In each panel, or equivalently, black hole mass bin, the power spectra of different R Edd levels are nearly uniformly spaced.Noting that the bins are split uniformly in log(R Edd ), this spacing implies that log(var) scales approximately linearly with log(R Edd ).For each mass bin, the power spectra show lower amplitudes for higher accretion rates, and the decrease is almost independent of frequency (in each plot, the power spectra are almost parallel to each other).However, a small steepening of the power spectra at higher accretion rates is apparent.This implies that the shape of the power spectrum depends to some extent on R Edd , either via slope or the bend frequency.
Figure 3 shows the dependence of the power spectra on black hole mass.Each panel contains one bin in R Edd to compare the power spectra for different masses, coded by colour, and similar R Edd .In all but the first panel the mass dependence becomes obvious, where higher masses correspond to increasingly steeper slopes.In the first panel, all slopes are similarly steep, but we note that the range of masses included is smaller and they are all in the high mass end of the sample.A similar steepening of the power spectra with increasing mass was observed by Caplar et al. (2017), for a large sample of about 28.000 quasars with a range of redshifts up to z ∼ 2.2, observed in the r-band by the Palomar Transient Factory survey.This dependence of the power spectrum slope on the black hole mass can be reconciled with a universal power spectrum shape for all quasars if this shape corresponds to a bending power law and the bend frequency scales with mass, with a higher bend frequency for lower black hole masses, as has been shown by Burke et al. (e.g. 2021).Therefore, if we plotted the frequency axes in Fig. 3 relative to the bend frequency of each mass bin (i.e.f / f b ), then the spectra of lower mass objects would shift to the left and that of higher mass objects would shift to the right.In principle, these shifts could make the differently coloured spectra align onto a single power spectrum, with a bending power-law shape.Burke et al. (2021) found that the bend timescale of the power spectrum scales with mass approximately as t b ∝ M 0.5 , and that the characteristic timescale for a 10 8 M⊙ black hole is about 200 days.Following this prescription we re-scaled the frequencies of our power spectra by multiplying them by their corresponding value of (M/M 8 ) 0.5 .The result of this scaling is shown in Fig. 4, colour coded by R Edd as in Fig. 2 but this time including all the mass bins in the same panel.There is a good alignment of the power spectra for different masses and equal R Edd , consistent with a universal power spectrum shape where the bend frequency scales with mass and the normalisation scales with R Edd , although with some scatter.The black solid line represents a bending power law model with α L = 0 and α H = −2 with the bend at [200 days] −1 , comparable to the DRW model with the bend frequency at the location found by Burke et al. (2021), plotted only as a guideline.The comparison between spectral shapes cannot be done directly in this plot, however, since the light curve sampling and variance estimation methods distort the spectra to some extent.Therefore, we will model the spectral shape using a "folded" model as described below.

Modelling of the Power Spectrum
Power spectra of AGN, in the optical as well as the X-ray regime, resemble bending or broken power laws, with at least one characteristic break frequency (e.g.McHardy et al. 2006;Burke et al. 2021;Stone et al. 2022, and references therein).This shape includes the power spectrum of a Damped-Random-Walk (DRW) model, often used for fitting optical power spectra of quasars (e.g. Kelly et al. 2009;Kozłowski et al. 2010;MacLeod et al. 2010MacLeod et al. , 2012;;Sánchez-Sáez et al. 2018;Suberlak et al. 2021;Burke et al. 2021;Stone et al. 2022), when the low-frequency power-law slope is 0 and the high-frequency slope is -2.This particular model is shown by the solid black line in Fig. 4, where the bend frequency was set at 1/(200 days) following the findings of Burke et al. (2021), and with an arbitrary normalisation.We note that although this model is close to the shape and apparent bend frequency of our power spectra, it does not reproduce perfectly the measured slopes.Other authors have also found deviations from the DRW model for other samples of quasars and other light curve sets (e.g.Mushotzky et al. 2011;Zu et al. 2013;Kasliwal et al. 2015;Caplar et al. 2017;Sánchez-Sáez et al. 2018;Stone et al. 2022).We will therefore use a more general bending power law model that allows for different power law slopes and is often used to model quasar power spectra in the X-ray regime (e.g.McHardy et al. 2004): In this model, frequencies below the break approximate a power law of slope α L and above the break a power law of slope α H .
As discussed in Sec. 4, for a given mass, the logarithm of the variance scales approximately linearly with the logarithm of R Edd .The black hole mass in turn appears to modify the break frequency, so that the frequency range probed, i.e. 1/300 to 1/30 days −1 covers part of the power spectra below the break for the lowest mass bins and above the break for the highest mass bins,  The power spectra depend on both parameters but in different ways, the mass changes the slope (or location of the bend in the power spectrum as described in Sec.4) and the accretion rate modifies mostly the normalisation, although some change in slope is apparent.Frequencies in are in units of day −1 .
again showing a uniform separation of log(variance) for a uniform separation of log(M).Additionally, there is a weaker dependence of the power spectrum shape on R Edd .Therefore, we will attempt to fit all 104 power spectral points shown in Fig. 4 together, with a bending power law model where the bend frequency and normalisation have the following dependencies: where f b and B will be in units of days −1 .We note that A1 and B will be the normalisation of the power spectrum and the bend frequency, respectively, for quasars of 10 8 M ⊙ and R Edd = 0.1, which are the most typical in our sample.
It is unclear whether a universal power spectrum would have equal variances (i.e.P × f ) at the bend frequency or at zero frequency, for a given R Edd and different masses.If the lowfrequency slope is α L = −1 then both cases are equivalent, but if the low-frequency slope is flatter then the choices are to align the low-frequency spectra or to match the variance at the break.We will explore the case where the universal power spectrum implies equal variance at the bend frequency.In this case, the normalisation A1 has an additional dependence on the break frequency: With this scaling, the variance, P × f will be equal to A2 × B at the bend frequency for all mass bins when R Edd = 0.1 and for the other R Edd ranges it will scale with R Edd only, i.e.
We will fit different power spectrum models to obtain the values of A2, B, C, D and F that best describe the data.

PDS model fitting with simulated lightcurves
The sampling pattern of the lightcurves and the method used to obtain the variance inevitably affect the measurements to some degree.Therefore, we resort to simulated light curves using the model defined in Eq.1 and time samplings taken from the real observations, to compare the measured variances with a 'folded' model.The light curves were simulated following the prescription of Timmer & Koenig (1995), with a total length of 1200 days and a generation sampling of 0.1 days, both in the quasar rest-frame.We allowed the model parameters to depend on the physical properties M and R Edd as described in Eqs. 2, 3 and 4. We searched a grid in the parameters B, C, D and F.
We note that the dependence of the normalisation on R Edd is set by the parameter F, but the overall normalisation is arbitrary.We set it by minimising the difference between the sum of the data variances and the sum of the model variances, to obtain the best-fitting value of A2 for each set of the other parameters.The steps in C, D and F were linear and the steps on B were multiplicative.Different step sizes were used to refine the search for the best fit, with a final step size of 1.2 in B, 0.05 in C, D and F.
For each set of model parameters, we calculated the break frequency f b and normalisation A for each M−R Edd bin, and created 100-300 simulated light curves with a fixed set of 30 real sampling patterns re-scaled to the quasars rest-frame times.These light curves were then used to measure the variance using the Mexican Hat method at the same four timescales as was done with the real data.We then averaged the logarithm of the resulting variances to obtain the 'folded model' variance for each M−R Edd bin and timescale.
The figure of merit, f om, recorded for each set of parameters was where ṽ ar is the folded model for one M−R Edd bin and timescale and var is the variance for the corresponding bin and timescale of the real data.The errors considered, err, are the errors on var of the real data, obtained through bootstrapping as described above.

PDS Model comparison
Optical light curves of quasars are often modelled as a DRW (e.g Kelly et al. 2009;Kozłowski et al. 2010;MacLeod et al. 2010;Burke et al. 2021), which produces a power spectrum slope of 0 below a break frequency and a slope of -2 above.To compare with previous works we include a model with a similar shape, i.e. α L = 0 and α H = −2.Guo et al. (2017)  They conclude that this effect indeed can explain the observed scatter and that the steepest the low-frequency slope can be is α L = −1.3,so we will include also steeper values of α L = −1 and α L = −0.5, i.e. within the acceptable range according to that study.Additionally, as in Fig. 4 the measured high-frequency power spectra appear steeper than the DRW model plotted in black, we will also use steeper high-frequency slopes α H = −2.5 and α H = −3.
We initially explored a grid in the parameters with B = 0.002 − 0.02 in multiplicative steps of 1.5, C = −0.8− −0.3 in steps of 0.1, D = −0.6 − 0.0 in steps of 0.1 and F = −1.0− 0.0 in steps of 0.05, creating 100 simulated light curves per M−R Edd bin to produce the folded model.If the minimum f om fell on a limiting value of one or more parameters, the fit was repeated with a wider parameter range.Table 2 shows the best f om obtained for each pair of slopes, together with the best-fitting values of all the parameters.
Considering the low resolution of the grid, the best-fitting parameters are fairly consistent among the different models, except for the bend frequency B, which is higher for the models with steeper high-frequency slopes, and consequently, the normalisation at the bend A2 differs.Although there is not much difference between the f om values between models with equal α L , there is a significant difference between models with different α H , where α H = −2.5 and α H = −3 fit consistently better than α H = −2, so we will not consider the latter models further.
A second round of simulations was done for the better-fitting models: all three models with α H = −2.5 and α L = −0.5, 0 for α H = −3.These fits consisted of significantly more iterations, with 300 simulations per M−R Edd bin and per set of parameters, with steps of 1.2 in B and 0.05 in C, D and F. The best-fitting values for all parameters and the minimum f om are presented in Table 3.Although models with f om values around 2 do not appear to be very good representations of the data, they are suf-  .In most R Edd bins, the power spectra steepen for higher masses, with the exception of the first panel, where all masses are high and all spectra are steep.
ficiently good if the errors on M and R Edd are considered, as discussed below in Sec.6.2.Therefore, we did not introduce further parameters to the model. .
Figure B.5 shows the minimum f om/N for fixed values of D on differently coloured lines, as a function of F. We recall that D determines the dependence of the bend frequency f b on R Edd as f b ∝R Edd D and F determines the dependence of the normalisation of the power spectrum on R Edd as A ∝ R Edd F .Since both parameters relate to the dependence of the power spectrum on R Edd they are partially degenerate, which produces the flat profiles of f om/N vs D and also f om vs F to a lesser extent.It is clear, however, that extreme values of D produce overall worse fits, for any value of F, and D = 0 (i.e.no dependence of f b on R Edd ) produces significantly worse fits.
The best-fitting values of parameters C, D and F can be used to scale the frequency and variance of the power spectra of f om/N 0 -2.5 1.5 0.0052 -0.65 -0.30 -0.50 2.0 -0.5 -2.5 1.4 0.0052 -0.60 -0.30 -0.45 1.9 -1 -2.5 1.1 0.0058 -0.all M−R Edd bins to that of the fiducial bin, i.e.M = 10 8 M ⊙ , R Edd =0.1.This scaling results in the power spectrum shown in large dots in Fig. 5.One of the best-fitting 'folded' models, with α L = −1 and α L = −3, is over-plotted in small black dots.The model follows the data quite closely and shows some scatter, produced by the sampling of the light curves and the stochastic nature of the 'folding' process.Figure 6 shows the difference between the model and data points, divided by the error on the data points, colour-coded by mass.Although some bins have a large deviation from the model, they are not preferentially of one R Edd or mass, as can be seen in the colour codes of both figures.

Error bounds on the fitted parameters
To estimate the confidence intervals of the parameters we start by considering the parameter values where f om = f om min + ∆, where ∆ is set to produce a given confidence interval assuming a χ 2 distribution with 4 degrees of freedom since 4 parameters are being searched simultaneously over a grid.For 90% confidence intervals, ∆ = 7.8 which equates to an increase in f om/N of 7.8/104=0.075.This value of ∆ is too small however, given that the folded model is calculated through a stochastic process that produces a larger uncertainty on the f om min of each realisation of the fitting routine.In other words, different realisations of the same model calculation, with the same model parameters result in slightly different values of the f om/N, with a scatter larger than this ∆ = 0.075.In Fig. 7 we show the distribution of values of f om/N obtained for the best-fitting parameters of each model in Table 3 for 100 realisations of each.For each realisation, 300 lightcurves were generated for each M−R Edd bin, to make the f om values comparable to the second round of simulations.Evidently, each model has a distribution of values of f om/N for the best-fitting parameters, with a standard deviation of 0.1-0.2.Therefore we will consider as acceptable values all parameter combinations that return a f om/N equal or lower than the 90% upper bound of this distribution of f om/N, for each model.The upper 90% bounds of f om/N are listed in the fourth column in Table 5 and represent an addition of 0.2 to 0.6 to f om min /N.These limit values of f om/N are marked by the grey horizontal lines in

Note on the figure of merit
The f om has been calculated using only the uncertainty in the variance of the binned data, not the uncertainty on the M or R Edd of the data bins.The latter errors should not be included for fitting purposes because their propagation to the error on the variance is model-dependent, so including them could bias the -2.5 0.0035 0.0055 -0.7 -0.53 -0.46 -0.12 -0.68 -0.3 -0.5 -2.5 0.0042 0.0058 -0.68 -0.54 -0.36 -0.16 -0.62 -0.42 -1 -2.5 0.0044 0.0081 -0.71 -0.51 -0.47 -0.02 -0.76 -0.31 -0.5 -3 0.007 0.0097 -0.63 -0.47 -0.4 -0.19 -0.56 -0.32 -1 -3 0.0095 0.0131 -0.67 -0.5 -0.41 -0.12 -0.64 -0.32  fit towards parameter values where the dependence of the variance on mass and R Edd is larger and therefore the errors are also larger and the f om becomes smaller.We note that all the fitting described above, as well as all the errors used for the figures, represent only the scatter on the measured variances in each data bin.The error on the mass and R Edd of the bins does affect the figure of merit since they will add to the difference between the model and data points but are not considered in the errors in the denominator in Eq. 5. To estimate how much the f om could improve if these errors were considered, we took the best-fitting parameters for each model, propagated the error on mass and R Edd measured from the scatter of these values in the bin, to an error on the model variance and combined these errors quadratically with the error on the variance from the data.Replacing these combined errors in the denominator of equation 5 reduced the figures of merit from the values in Table 3 to the values listed in the fifth column in Table 5 ( f om scatter /N).
Alternatively, we can use the uncertainties in the values of M and R Edd for individual objects to estimate the errors on M and R Edd of each bin, not just the scatter in reported values as above.We assumed a systematic error of 0.3 dex in the single epoch mass estimate and an error of 50% in R Edd as noted by Rakshit et al. (2020) to be relevant errors for their catalogue.We calculated dM = M/ √ N lc and dR Redd = 0.5R Redd / √ N lc , where N lc is the number of light curves in the relevant M−R Edd bin and dM and dR Redd will be the errors in the mass and R Edd of the bin, respectively.Propagating these errors and then combining them with the uncertainty in the variance of each bin produced the values listed in the sixth column in Table 5 ( f om ind /N).These errors on mass and R Edd of the individual measurements in the catalogue can therefore explain the remaining difference between the model and data points.

Discussion
The data presented above show that the optical power spectra, at around 2900 Å in the rest-frame, conform with a universal bending power-law spectrum, where the bend frequency scales inversely with both M and R Edd and the normalisation scales inversely with R Edd .Comparing the dependence of the power spectra on R Edd (Fig. 2) to the dependence of the power spectra on mass (Fig. 3) we see their very distinct effects: the R Edd mostly changes the normalisation while the mass has the largest effect on the location of the break in the power spectrum.As a result, the mass determines which part of the power spectrum is being sampled by the light curves, resulting in different slopes visible in the frequency range probed.An immediate consequence of this is that the power spectrum cannot be simply a function of luminosity, which is the product of the two.For example, a lowmass, high-accretion rate object will have a flatter power spectrum with lower normalisation than a high-mass, low-accretion rate object of the same luminosity.Below we discuss the implications of these results in terms of the accretion disc structure and other considerations for comparing these results to other data sets.

Interpretation of the scaling relations
The orbital timescale at a given location in the disc, for example at the ISCO, scales linearly with black hole mass.Re-scaling the break timescales by mass, so that they can be interpreted in terms of the orbital timescales at a given radius in units of R g , results in t b /t IS CO ∝ M 0.55 R Edd 0.35 /M = (R Edd /M) 0.35 M −0.1 , where we have used the best fitting parameters for the model with α L = −1 and α H = −3.Alternatively, a value of C = −0.65 is allowed within the uncertainty range for the same model.With this value, the relation becomes t b /t IS CO ∝ M 0.65 R Edd 0.35 /M = (R Edd /M) 0.35 , i.e. a bend timescale that scales linearly with the orbital timescale of the innermost orbit times a function of the combination R Edd /M.In the standard thin disc model (Shakura & Sunyaev 1973) the maximal surface temperature T is proportional to (R Edd /M) 0.25 , suggesting a connection between the characteristic timescale and the disc temperature.
The size of the emitting region of a given wavelength λ * can be estimated by the largest radius in the disc where the temperature is above T * = 0.29/λ * with λ * in cm and T in K.The orbital timescale at this radius is t * orb = 2πr * 3/2 R g /c or equivalently t * orb = r * 3/2 t IS CO , where the dimensionless radius r = R/R g and R g = GM/c 2 .Assuming a radial temperature dependence of This relation between the orbital timescale at the outer edge of the optical emitting region, t * orb , and M and R Edd is slightly steeper than the relation we find between these quantities and the power spectrum bend timescale We note that re-replacing t IS CO by M in Eq. 7 results in t * orb ∝ λ * 2 (R Edd M) 0.5 ∝ λ * 2 L 0.5 bol , as calculated for example by MacLeod et al. (2010).Our fits to the data can be expressed in terms of bol R Edd 0.2 , i.e. the difference between C and D implies that the bend timescale cannot depend on L bol alone.We note that although both parameters C and D have associated errors, the allowed ranges for C and D values do not overlap, for any of the fitted models, as shown in Table 4.
Returning to the standard thin disc model, if we characterise the optical emitting region by a light-weighted radius, rather than the outermost radius, the relation between orbital timescale and mass and R Edd changes somewhat.We used thin disc model prescriptions for the emitted flux as a function of mass, R Edd and radius to obtain the black body temperature as a function of radius, T (r).We start from the flux emitted by a thin accretion disc (e.g.eq.23 in Spruit 1995): and assume that the mass accretion rate Ṁ in cgs units corresponds to: Replacing this value in equation 8 to obtain a relation in terms of the mass in solar units m = M/M ⊙ and the Eddington ratio R Edd =L/L Edd gives: 1.26 × 10 38 mR Edd ηc2 (10) where r = R/R g , R g = GM/c 2 = 1.5 × 10 5 m cm and σ is the Stefan Boltzmann constant.Then the energy flux per unit frequency interval emitted by (both sides of) a ring of area dA = 2πrdrR 2 g at a given frequency ν is These annular fluxes were used to calculate a light-weighted average radius as For these calculations we assumed f = 1 − √ r min /r, r min = 6 and η = 0.1.The outer radius of the emission region, in units of the gravitational radius R g , was set at r out = 3000.
For a given mass, the relation between the light-weighted radius r lw and R Edd is slightly flatter than the relation of the outermost radius and R Edd .In Fig. 8 we show the orbital timescale, i.e. t orb = 2πr 3/2 lw R g /c, with r lw in units of R g , as a function of both mass and R Edd , calculated for parameter ranges similar to those of our data.For a given R Edd , or colour in the top panel in Fig. 8, the relation between orbital timescale and mass follows closely a t orb ∝ M 0.65 relation, plotted as a blue solid line, and different masses show approximately parallel relations.The bottom panel displays the relation between t orb and R Edd which, for a given mass, closely follows a t orb ∝ R Edd 0.35 relation, plotted in the solid line, with approximately parallel relations for other masses.Therefore, the orbital timescale at a light-weighted radius for the optical emitting region scales with M and R Edd in approximately the same way as the bend timescale.
As can be seen in Fig. 8 the orbital timescale at the lightweighted radius of 2900 Å, for the standard thin disc, and for M = 10 8 M ⊙ and R Edd = 0.1 is about 10 days.The bend timescale t b = 1/ f b that we find for these black hole parameters is between 90 and 200 days, depending on the power spectrum model.The ratio between the bend and orbital timescales is therefore 9 to 20.As pointed out before by several authors (e.g. Kelly et al. 2009;MacLeod et al. 2010;Burke et al. 2021), the bend timescale could correspond to the thermal timescale of a thin disc, which is expected to be around 10 times the orbital timescale, independently of the accretion rate and most disc parameters (e.g Spruit 1995, and references therein).
We note however that this modelling with a standard thin disc produces a bolometric correction factor k 5100 = L bol /L 5100 that depends on the disc temperature, and therefore also on the ratio R Edd /M.The values of R Edd used for the fitting have been calculated with a constant bolometric correction factor k 5100 = 9.26 (Rakshit et al. 2020).In Appendix A we explore the effect of correcting the values of R Edd to the values predicted by the thin disc model, given the L 5100 and mass of each quasar, to produce a more self-consistent model comparison.The result is that with this correction, the relation between f b and mass is stronger ( f b ∝ M −0.66 instead of f b ∝ M −0.55 ) and the relation of f b to the "true" R Edd is weaker, but still of the same sign, i.e. f b ∝ R Edd −0.25 model ), for the model with α L = −1 and α H = −3, still fairly consistent with the model scalings described above.
Apart from the scaling of the bend frequency, we also find that the normalisation of the power spectrum at the break scales with R Edd approximately as A ∝ R Edd −0.45 .In the interpretation described above this would mean that smaller optical emitting regions characteristic of lower R Edd sources and therefore also lower temperature discs, are not only characterised by shorter timescales but also by larger amplitude of fluctuations.This scaling might help in constraining accretion disc models that incorporate prescriptions for its variations in luminosity such as the inhomogeneous disc model with thermal perturbations in clumps proposed by Dexter & Agol (2011) and expanded by Cai et al. (2016) or the Corona Heated Accretion Disc model of Sun et al. (2020), perhaps by having fewer independently-varying regions that emit in the observed wavelength for cooler discs.
Incidentally, relating the bend frequency to the size of the emitting region implies that the bend frequency should be wavelength-dependent, where shorter wavelengths, being produced by a smaller region, would have higher bend frequencies.Stone et al. (2023) recently found a relation between rest-frame wavelength and the de-correlation timescale of the DRW model, τ DRW , but with a weaker dependence of τ DRW ∝ λ 0.4+/−0.1 , from very long (20 years) light curves and DRW modelling of 190 quasars.As they argue, this weaker dependence of timescale with λ can be expected, for example, for a radially-averaged thermal timescale, as described for the model of Sun et al. (2020).
Finally, we note that, as has been shown in every power spectrum analysis before (e.g de Vries et al. 2005;Caplar et al. 2017;Burke et al. 2021;Stone et al. 2022), there is significant variability power on timescales much shorter than the bend timescale.In other words, the variability is not simply concentrated on one characteristic timescale but, instead, there is variability on a continuum of timescales, with smaller amplitudes for increasingly shorter timescales.In the interpretation described above, this could be achieved if all the optical emitting regions produced variations on their local thermal timescales.In this case, the smaller radii in the disc would produce shorter timescale fluctuations and would modulate a smaller fraction of the 2900 Å emission, both by having a smaller area and by having a higher temperature, so that a larger fraction of their luminosity is emitted at shorter wavelengths than observed. 1

Other implications
We note that, in our modelling, the bend frequency of the power spectrum scales both with mass and with R Edd but the dependence on R Edd is about half as strong as the dependence on mass.Therefore, for smaller or more heterogeneous samples, the dependence on R Edd might not be noticeable.One effect of not accounting for this dependence is that it will 'leak' into the dependence of f b on mass, for example, if mass and R Edd are anticorrelated in the sample.In the most extreme case, where ob-  jects in a sample have M ∝R Edd −1 then, for example, 3 , flattening the dependence on mass.In less well anti-correlated samples, such as the one used in this paper, the effect would be weaker but still present.
The best fitting parameters and f om obtained by fixing D = 0 (i.e.neglecting any dependence of the bend frequency on R Edd ) are summarised in Table 6.Comparing Table 6 to Table 3, where D was a free parameter exemplifies the effect described in the previous paragraph, where neglecting the dependence of f b on R Edd results in a shallower dependence of f b on mass, i.e. a smaller value of |C|.We note also that all these models fit the data significantly worse than models with free D. However, these fits are useful to show the overall dependence of variance as a function of accretion rate.In all models, this relation results in either P × f ∝R Edd −0.75 or P × f ∝R Edd −0.7 .Although part of this dependence results from the scaling of power spectrum normalisation with R Edd and the rest from the dependence of the bend frequency on R Edd , this relation is useful to predict the variance expected for equal mass and different accretion rates, at least at timescales shorter than the bend in the power spectrum.

The shape of the Power Spectrum
The universal power spectrum model fits our data well, particularly for a high-frequency slope of α H = −2.5, which allows low-frequency slopes of α L = −0.5 or 0, or for α H = −3 and α L = −1, all with similar goodness of fit.et al. (2011) using Kepler light curves, among others, all of which showed high-frequency power spectra that are steeper than α H = −2.Regarding the low-frequency slope, a value of α L = −1 was preferred in our analysis only for the model with α H = −3, while for α H = −2.5 any low-frequency slope (i.e.α L = 0, −0.5, −1) was acceptable.This result probably relates to the location of the model bend, which is at a higher frequency in the case of α H = −3 and therefore the low frequency slope is better constrained in this model.We note that a similarly steep lowfrequency slope of α L = −1 was found by Simm et al. (2016) from Pan-STARRS-1 data of 90 quasars using CARMA (Kelly et al. 2014) modelling.
All our models with α H = −2.5 gave bend frequencies with overlapping confidence ranges, and the models with α H = −3 resulted in significantly higher bend frequencies.This result might indicate that a continuously bending power-law is a better description of the true power spectrum than a single-bend powerlaw model since steeper slopes fit better with bends at higher frequencies.

Conclusion
We present a power spectral analysis of optical light curves taken in the g-band by the ZTF project using 5433 light curves for 4770 individual quasars at 0.6<z<0.7.This narrow redshift range allows us to probe the dependence of the power spectral parameters on mass and R Edd without having to correct for uncertain dependencies of the power spectrum on the rest-frame wavelength of emission.The light curves have a median length of 1547 days and a median of 233 good data points per light curve.The masses and R Edd of all the quasars used in this study were obtained by Rakshit et al. (2020) by modelling SDSS spectra homogeneously.
We used the Mexican hat method (Arévalo et al. 2012) to estimate the variance of the light curves at 4 different timescales, in the range of 30 to 300 days in the rest-frame of the quasars (already presented in Arévalo et al. 2023).With these variance estimates, we construct low-resolution power spectra for 26 bins in a grid of mass and R Edd , with 7.5 < log M/M ⊙ < 9.5 and −2 < R Edd < 0. Each bin had a minimum of 14 lightcurves.The power spectra can be approximately described by a power-law shape with a negative slope.When separated in groups of similar R Edd the power spectra show a clear dependence on mass, where higher mass objects show steeper power spectra and the spacing between high-frequency log(power) grows approximately linearly with the separation in log(M).When split in groups of similar mass, the power spectra of different R Edd appear almost parallel, with higher normalisation for lower R Edd .In this case, too the spacing between the logarithm of the power spectra scales approximately linearly with the separation in log(R Edd ).There is still a small dependence of the slope on R Edd but this is less marked than the dependence on mass.
The power spectrum models tested were smoothly bending power-law models, characterised by a bend frequency f b , lowand high-frequency power-law slopes α L and α H and a normalisation.The f b was allowed to depend on the mass and R Edd and the normalisation on R Edd , based on the findings described in the previous paragraph.We note that the power spectrum model used had the same variance at the bend frequency for objects of different mass and equal R Edd .For low-frequency slopes α L −1 this normalisation implies that the power at 0 frequency is different for objects of different masses.The best-fitting bend timescales, for M = 10 8 M ⊙ and R Edd = 0.1, range from 96 +9 −20 days for the model with α H = −3 and α L = −1 to 192 +46 −10 days for the model with α H = −2.5 and α L = −0.5.The latter model is the closest to the DRW used in Burke et al. (2021) and for this model, our bend timescales are consistent with their results.The best fitting slope pairs were α H = −3, α L = −1 and α H = −2.5 with either α L = −0 or α L = −0.5.DRW-type high-frequency slopes α H = −2 gave significantly worse fits.
We find, for the first time, a dependence of the break frequency on R Edd in the optical power spectra of quasars, with the bend frequency scaling as f b ∝ M −0.67−−0.5 ×R Edd −0.41−−0.12, where the ranges in the exponent values correspond to a powerspectral model with α L = −1 and α H = −3, but they are similar to the ranges of the other well-fitting models.As the exponents for the dependence on mass and on R Edd are different, the bend timescale cannot be simply a function of bolometric luminosity, which is the product of both parameters.Also as the exponent of the dependence on R Edd is significantly different from 0, the bend timescale is not a function of mass alone.Interestingly, the higher R Edd objects show less high-frequency power, this sense is opposite to the behaviour of the X-ray power spectra of quasars, where higher accretion rate objects have more high-frequency variability (e.g.McHardy et al. 2006, and references therein).This can be explained if the bend timescale corresponds to the thermal timescale of the light-weighted radius of the emitting region of the observed in the context of the standard thin disc model.In this case, for a given mass and different accretion rates, the characteristic radius increases with temperature and therefore with accretion rate and so the relevant timescale becomes longer.A simple calculation with thin disc prescriptions shows that this light-weighted radius for emission at 2900 Å has an orbital timescale that scales with M 0.65 and R Edd 0.35 , similar to the scaling we find for t b = 1/ f b , at least for masses and accretion rates in the same range as the data.
The dependence of f b on mass goes from f b ∝ M −0.65 to f b ∝ M −0.55 depending on the power spectrum model and is slightly steeper than that previously published by Burke et al. (2021).This might be partly due to the dependence of f b on R Edd and the anti-correlation between mass and R Edd that usually exists in flux-limited quasar samples which, combined, results in a flatter mass dependence if the dependence of f b on R Edd is not included.Alternatively, a flatter mass dependence could result from a positive relation between bend timescale and wavelength.High-mass objects are preferentially found at higher redshifts, so their rest-frame wavelength is shorter than for lower-mass objects in the same sample observed with the same filter.If the bend timescale indeed scales with wavelength as predicted by the thin disc model, then the high mass, high redshift objects would appear to have shorter bend timescales than if they were observed in the same rest-frame wavelength, reducing the difference with the low mass, low redshift objects.
Measuring the wavelength dependence of the bend timescale is challenging since the effect is subtle, the uncertainties on the bend timescale are large, and the dependence on the other black hole parameters must be controlled in a sufficiently accurate manner.Such a correlation, however, has been found using the characteristic timescale of variability measured with the Structure Function, by MacLeod et al. (2010, and references therein) where the reported dependence is very weak.A similar analysis but estimating the characteristic timescale with the DRW model (τ DRW ) and limiting the study to objects whose τ DRW is much shorter than the length of the light curves was presented recently by Stone et al. (2022).This later work found a stronger correlation with rest-frame wavelength but still weaker than predicted by the thin disc model.Long-term, well-sampled light curves of many more objects, in multiple optical bands and to high depths, as will be provided in the future by the Vera Rubin LSST (Ivezić et al. 2019) will likely shed more light on this issue.
It is also possible that not only the characteristic timescale of variability depends on wavelength but that the slope of the power spectrum might depend on this parameter as well.Results in this direction will be presented in the future by Patel et al. (in prep.) who calculate the ratio of variances on different timescales (i.e. the power spectrum slope) as a function of rest-frame wavelength for objects of equal mass and R Edd .In addition, a replication of the work presented in this paper but using the other available ZTF band (r) will show whether the fitted power spectrum parameters depend strongly on wavelength.

Fig. 1 .
Fig. 1.Low-resolution power spectra for different values of R Edd within a narrow range in black hole mass (left) and different black hole masses but similar values of R Edd (right).The markers represent the median variance (i.e.Power× frequency) of all quasars falling in each M−R Edd bin and temporal frequency, and the error bars represent the root-mean-squared (rms) of the median for each bin obtained by bootstrapping the samples.The power spectra depend on both parameters but in different ways, the mass changes the slope (or location of the bend in the power spectrum as described in Sec.4) and the accretion rate modifies mostly the normalisation, although some change in slope is apparent.Frequencies in are in units of day −1 .

Fig. 3 .
Fig. 3. Low-resolution power spectra (plotted as Power × frequency) for different R Edd bins in different panels and different black hole masses in different colours.The markers represent the median variance of all quasars falling in each M−R Edd bin and temporal frequency, and the error bars represent the rms of the median for each bin obtained by bootstrapping the samples.Frequencies in are in units of day −1

Fig. 4 .
Fig.4.Power spectra (plotted as Power × frequency) of all the sources, binned by mass and Eddington ratio.The variances have been calculated for timescales of 30, 75, 150 and 300 days, in this plot the corresponding frequencies f = 1/t have been multiplied by (M/M 8 ) 0.5 to put them on a common x-axis scaled to the bend timescale approximately as measured byBurke et al. (2021).The colour represents the different bins in R Edd .The black solid line represents a bending power law model with α L = 0 and α H = −2, similar to the DRW model, with a bend timescale of 200 days for a 10 8 M ⊙ black hole.(We note that as power is multiplied by frequency in the plot, the part of the black line with α L = 0 appears rising with frequency).Frequencies are in units of day −1

Fig. 5 .
Fig. 5.The large dots, colour-coded by R Edd represent the power spectra of the data (plotted as Power × frequency), after the frequencies and amplitudes have been re-scaled according to the mass and R Edd of each bin and the best-fitting scaling parameters C=-0.55,D=-0.35,F=-0.4.The black dots represent a 'folded' model constructed from a bending power law model with α L = −1 and α H = −3.Frequencies are in units of day −1 .

Fig. 6 .
Fig. 6.Residuals of the model and data plotted in Fig.5, color coded by mass.Frequencies are in units of day −1 .
Distributions of values of f om/N obtained for 100 different realisations of each model, using only the best-fitting parameters in each case.The different models are characterised by the low and highfrequency slopes of their power spectra and these are written in the legend.

Fig. 8 .
Fig. 8. Standard thin disc model predictions for the orbital timescale at the light-weighted radius for emission at 2900Å, as a function of black hole mass and R Edd .The top panel shows the timescales as a function of M, color coded by R Edd , together with a solid line displaying a t orb ∝ M 0.65 relation.The bottom panel shows the same model points as a function of R Edd , color coded by mass and the solid line displays a t orb ∝R Edd 0.35 relation.

P
figure of merit

Table 1 .
Number of light curves included in each M−R Edd bin.The line and column headings display the ranges in M and R Edd , respectively, that limit each bin.
modelled the scatter in observed DRW parameters of SDSS Stripe 82 quasars by allowing a low-frequency slope that is steeper than the DRW model.P. Arévalo et al.: Scaling of the quasar optical power spectra Fig. 2. Similar to Fig. 2 but arranged in separate panels for different mass bins, and different colours used to differentiate by accretion rate level.

Table 2 .
Best-fitting parameters and corresponding figure of merit for the 9 power spectrum models tested, with low and high-frequency slopes given in the first two columns.The parameters are explained in Eqs.Eqs.1,2,3,4 and N is the number of power spectrum points, i.e. the number of timescales times the number of M−R Edd bins.First, coarser grid run.

Table 3 .
Best-fitting parameters and corresponding figure of merit for the 5 better-fitting power spectrum models, with low and highfrequency slopes given in the first two columns.The parameters are explained in Eqs.1,2,3,4 and N is the number of power spectral points, i.e. the number of timescales times the number of M−R Edd bins.Second, finer grid run.

Table 4 .
Lower and upper limits for the power spectrum model parameters, for each combination of low and high-frequency slopes written in the first two columns.The break frequencies B are in units of days −1 .

Table 5 .
om 90 /N f om scatter /N f om ind /N Figure of merit for each model in Table3, f om/N corresponds to the best fitting model, where errors in the denominator only consider the scatter in variances in the data bins; f om 90 /N are the upper 90% bounds of the distribution of f om/N values for different realisations of each model with its best-fitting parameters, described in Sec.6.1; f om scatter /N and f om ind /N are re-scaled versions of f om/N, incorporating also the errors on mass and R Edd , first using only the scatter of these parameters within each bin ( f om scatter /N) and then propagating the errors on individual measurements of both values ( f om ind /N) as described in Sec.6.2.
.Arévalo et al.:Scaling of the quasar optical power spectra P Models with a flatter high-frequency slope of α H = −2, similar to the DRW model, resulted in significantly worse fits.This last result is consistent with previous findings by Caplar et al. (2017) who used lightcurves from the Palomar Transient Factory (PTF/iPTF) surveys as well as Stone et al. (2022) using light curves from SDSS, PanSTARRS-1 and the Dark Energy Survey and Mushotzky