© The Authors 2024

1. Introduction

The flux variability of quasars was recognised early on as one of the few tools with which 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 light curves 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 the black hole mass, as obtained from the fit of damped random walk (DRW) models to some 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 light curves or as the de-correlation timescale of a DRW model (Kelly et al. 2009), can depend on the rest-frame wavelength of the light curves. 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 the 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 Sect. 2) to avoid potential dependencies of the power spectrum shape on the rest-frame wavelength of emission. We calculated the median power spectra for bins in the mass, M, and R_Edd, to establish the dependence of the power spectrum break frequency and normalisation on both parameters. The variance was 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 Sect. 3) and the resulting power spectra are presented in Sect. 4. Model fitting to a bending power law model, with different slopes, was 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 Sect. 5. Best-fitting models and scaling parameters were searched by minimising the difference between the data and the model power spectra. With this minimisation we show that the scaling relations of the amplitude and the 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 Sect. 6. A basic interpretation of the scaling relations in terms of the accretion disc parameters is presented in Sect. 7 and the conclusions are summarised in Sect. 8.

2. 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 had their mass and R_Edd estimated through Sloan Digital Sky Survey (SDSS) spectra by Rakshit et al. (2020), with reported statistical errors on the masses lower 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 two 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 233 points. We note that since all galaxies are at z > 0.6 they appear almost point-like in the ZTF images. The PSF photometry on the science images, which is used to produce the light curves of the data release, is therefore sufficiently accurate.

3. 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 four well-separated 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 corrected these noise estimates by carrying out the same procedure of variance estimation and noise subtraction for ZTF light curves 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 ten or more light curves, 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).

4. 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 Sect. 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 the 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 mass range (left), to the dependence of the power spectra on mass, for quasars within a narrow R_Edd range (right). Evidently, the power spectrum depends on both parameters but in different ways: the Eddington ratio affects mostly the normalisation and the black hole mass affects mostly the slope, in the frequency range sampled by the light curves.

Fig. 1. Low-resolution power spectra for different values of REdd within a narrow range in black hole mass (left) and different black hole masses but similar values of REdd (right). The markers represent the median variance (i.e., power × frequency) of all quasars falling in each M − REdd 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 the location of the bend in the power spectrum, as described in Sect. 4) and the accretion rate modifies mostly the normalisation, although some change in the slope is apparent. Frequencies are in units of day−1.

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 an 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 the slope or the bend frequency.

Fig. 2. Low-resolution power spectra (plotted as power × frequency) for different mass bins in different panels and different accretion rate levels in different colours. The markers represent the median variance of all quasars falling in each M − REdd bin and temporal frequency, and the error bars represent the rms of the median for each bin obtained by bootstrapping the samples. Frequencies are in units of day−1.

Figure 3 shows the dependence of the power spectra on the 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 the 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.

Fig. 3. Similar to Fig. 2 but arranged in separate panels for different REdd bins, and with different colours used to differentiate by mass. In most REdd 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.

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. (2021). Therefore, if we plotted the frequency axes in Fig. 2 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) find 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⁸ 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₈)^0.5. The result of this scaling is shown in Fig. 4, colour coded by R_Edd as in Fig. 3 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 solid black line represents a bending power law model with α_L = 0 and α_H = −2 with the bend at [200 days]⁻¹, 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 modelled the spectral shape using a “folded” model, as is described below.

Fig. 4. Power spectra (plotted as power × frequency) of all the sources, binned by mass and Eddington ratio. The variances were calculated for timescales of 30, 75, 150, and 300 days, and in this plot the corresponding frequencies, f = 1/t, were multiplied by (M/M8)0.5 to put them on a common x-axis scaled to the bend timescale approximately, as was measured by Burke et al. (2021). The colour represents the different bins in REdd. The solid black 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 108 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 to rise with frequency.) Frequencies are in units of day−1.

5. 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 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. 2010, 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 therefore used 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):

$$\begin{array}{c}\hfill P(f)=\frac{A{f}^{{\alpha }_{\mathrm{L}}}}{1+{(f/f_{\mathrm{b}})}^{{\alpha }_{\mathrm{L}}-{\alpha }_{\mathrm{H}}}}.\end{array}$$(1)

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 is discussed in Sect. 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, that is, 1/300–1/30 days⁻¹ covers part of the power spectra below the break for the lowest-mass bins and above the break for the highest-mass bins, 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 attempted 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:

$$\begin{array}{cc}& log{f}_{\mathrm{b}}=logB+C\times logM/M_8+D\times log{R}_{\mathrm{Edd}}/0.1\hfill \end{array}$$(2)

$$\begin{array}{cc}& logA=logA1+F\times log{R}_{\mathrm{Edd}}/0.1,\hfill \end{array}$$(3)

where f_b and B are in units of days⁻¹. We note that A1 and B are the normalisation of the power spectrum and the bend frequency, respectively, for quasars of 10⁸ 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 low-frequency 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 explored 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:

$$\begin{array}{c}\hfill A1=A2\times B/f_{\mathrm{b}}\times f_{\mathrm{b}}^{-a_{\mathrm{L}}}.\end{array}$$(4)

With this scaling, the variance, P × f, is 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 scales with R_Edd only, that is, P(f_b)×f_b = A2 × B× (R_Edd/0.1)^F.

We fitted 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 light curves. The sampling pattern of the light curves and the method used to obtain the variance inevitably affect the measurements to some degree. Therefore, we resorted 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 is described in Eqs. (2)–(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, fom, recorded for each set of parameters was

$$\begin{array}{c}\hfill \mathrm{fom}={\displaystyle \sum _{\mathrm{bins},\mathrm{timescales}}\frac{{(\stackrel{\sim }{\mathrm{var}}-\mathrm{var})}^2}{{\mathrm{err}}^2},}\end{array}$$(5)

where $\stackrel{\sim }{\mathrm{var}}$ 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 was described above.

6. 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) 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. They conclude that this effect can indeed explain the observed scatter and that the steepest the low-frequency slope can be is α_L = −1.3, so we also included steeper values of α_L = −1 and α_L = −0.5, in other words 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 also used 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 to −0.3 in steps of 0.1, D = −0.6 to 0.0 in steps of 0.1, and F = −1.0 to 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 fom fell on a limiting value of one or more parameters, the fit was repeated with a wider parameter range. Table 2 shows the best fom 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 fom 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 did 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 the parameters and the minimum fom are presented in Table 3. Although models with fom values of around two do not appear to be very good representations of the data, they are sufficiently good if the errors for M and R_Edd are considered, as is discussed below in Sect. 6.2. Therefore, we did not introduce further parameters to the model.

The minimum fom/N as a function of the bend frequency, B, is shown in Fig. B.1; in other words, for each value of B, the fom/N was marginalised over all the other parameters. The different panels show the results for the five different power spectrum models shown in Table 3. Evidently, the bend frequency shows a well-defined minimum in fom/N, within the accessible frequency range. The best-fitting bend frequencies are in the range, f_b = 0.005 − 0.01 days⁻¹, for M = 10⁸ M_⊙ and R_Edd = 0.1, being higher for steeper power law models. Similar plots for the other parameters are shown in Figs. B.2–B.4. Although these show shallower minima, they do show a constrained value different from 0 for all parameters.

The best-fitting values for the mass dependence of the bend frequency, f_b ∝ M^C, range from C = −0.65 for α_L = 0, α_H = −2.5 to C = −0.55 for α_L = −1, α_H = −3, with C = −0.6 for the power spectrum models with intermediate slopes. The best-fitting value for the dependence of the bend frequency on R_Edd, $f_{\mathrm{b}}\propto R_{\mathrm{Edd}}^D$ are either D = −0.3 or D = −0.35 for all models. Upper and lower limits for all the parameters, calculated as explained below in Sect. 6.1, are listed in Table 4. Notably, all models produce consistent values of C, in the sense that all their uncertainty ranges overlap, and also of D, while, for each model, C and D are significantly different, that is, the scaling of the bend frequency with mass is different to its scaling on R_Edd. Finally, the scaling of the normalisation on R_Edd, $A\propto R_{\text{Edd}}^F$, has best-fitting values of F = −0.4 to F = −0.5 and the uncertainty ranges for this parameter overlap for all models.

Figure B.5 shows the minimum fom/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_{\text{b}}\propto R_{\text{Edd}}^D$ and F determines the dependence of the normalisation of the power spectrum on R_Edd as $A\propto R_{\text{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 fom/N versus D and also fom versus F to a lesser extent. It is clear, however, that extreme values of D produce worse overall 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 all M − R_Edd bins to that of the fiducial bin, that is, M = 10⁸ 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 and 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.

Fig. 5. The large dots, colour-coded by REdd, 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 REdd of each bin and the best-fitting scaling parameters, C = −0.55, D = −0.35, and 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. Residuals of the model and data plotted in Fig. 5, colour-coded by mass. Frequencies are in units of day−1.

6.1. Error bounds on the fitted parameters

To estimate the confidence intervals of the parameters, we started by considering the parameter values where fom = fom_min + Δ, where Δ was set to produce a given confidence interval assuming a χ² distribution with four degrees of freedom, since four parameters were being searched simultaneously over a grid. For 90% confidence intervals, Δ = 7.8, which equates to an increase in fom/N of 7.8/104 = 0.075.

This value of Δ is too small however, given that the folded model was calculated through a stochastic process that produces a larger uncertainty on the fom_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 fom/N, with a scatter larger than this Δ = 0.075. In Fig. 7 we show the distribution of values of fom/N obtained for the best-fitting parameters of each model in Table 3 for 100 realisations of each. For each realisation, 300 light curves were generated for each M − R_Edd bin, to make the fom values comparable to the second round of simulations. Evidently, each model has a distribution of values of fom/N for the best-fitting parameters, with a standard deviation of 0.1–0.2. Therefore we consider as acceptable values all parameter combinations that return a fom/N equal or lower than the 90% upper bound of this distribution of fom/N, for each model. The upper 90% bounds of fom/N are listed in the fourth column in Table 5 and represent an addition of 0.2–0.6 to fom_min/N. These limit values of fom/N are marked by the horizontal grey lines in Figs. B.1–B.4.

Fig. 7. Distributions of values of fom/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 high-frequency slopes of their power spectra and these are written in the legend.

6.2. Note on the figure of merit

The fom was 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 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 fom becomes smaller. We note that all the fitting described above, as well as all the errors used for the figures, represents 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 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 fom 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 Eq. (5) reduced the figures of merit from the values in Table 3 to the values listed in the fifth column in Table 5 (fom_scatter/N).

Alternatively, we can use the uncertainties in the values of M and R_Edd for individual objects to estimate the errors in 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 was noted by Rakshit et al. (2020) to be relevant errors for their catalogue. We calculated $\mathrm{d}M=M/\sqrt{N_{\mathrm{lc}}}$ and $\mathrm{d}{R}_{\mathrm{Redd}}=0.5R_{\mathrm{Redd}}/\sqrt{N_{\mathrm{lc}}}$, where N_lc is the number of light curves in the relevant M − R_Edd bin and dM and dR_Redd are 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 (fom_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.

7. 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. 3) to the dependence of the power spectra on mass (Fig. 2) 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 low-mass, 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 datasets.

7.1. Interpretation of the scaling relations

The orbital timescale at a given location in the disc, for example at the ISCO, scales linearly with the 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_{\mathrm{b}}/t_{\mathrm{ISCO}}\propto M^{0.55}{R_{\mathrm{Edd}}}^{0.35}/M={(R_{\mathrm{Edd}}/M)}^{0.35}{M}^{-0.1}$, where we 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_{\mathrm{b}}/t_{\mathrm{ISCO}}\propto M^{0.65}{R_{\mathrm{Edd}}}^{0.35}/M={(R_{\mathrm{Edd}}/M)}^{0.35}$, in other words 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_{\text{orb}}^{*}=2\pi r^{*3/2}{R}_{\text{g}}/c$ or equivalently $t_{\text{orb}}^{*}=r^{*3/2}{t}_{\text{ISCO}}$, where the dimensionless radius, r = R/R_g and R_g = GM/c². Assuming a radial temperature dependence of T(r) = T_maxr^−3/4 gives

$$\begin{array}{c}\hfill t_{\mathrm{orb}}^{\ast }={(T^{\ast }/T_{\max })}^{3/2\times -4/3}{t}_{\mathrm{ISCO}}={(T^{\ast }/T_{\max })}^{-2}{t}_{\mathrm{ISCO}}.\end{array}$$(6)

If T_max ∝ (R_Edd/M)^0.25 then

$$\begin{array}{c}\hfill t_{\mathrm{orb}}^{\ast }\propto T^{\ast -2}{(R_{\mathrm{Edd}}/M)}^{0.5}{t}_{\mathrm{ISCO}}\propto {\lambda }^{\ast 2}{(R_{\mathrm{Edd}}/M)}^{0.5}{t}_{\mathrm{ISCO}}.\end{array}$$(7)

This relation between the orbital timescale at the outer edge of the optical emitting region, $t_{\text{orb}}^{*}$, and M and R_Edd is slightly steeper than the relation we find between these quantities and the power spectrum bend timescale, $t_{\mathrm{b}}\propto {R_{\mathrm{Edd}}}^{-D}{M}^{-C-1}{t}_{\mathrm{ISCO}}={(R_{\mathrm{Edd}}/M)}^{0.35}{M}^{-0.1}{t}_{\mathrm{ISCO}}$.

We note that re-replacing t_ISCO by M in Eq. (7) results in $t_{\mathrm{orb}}^{\ast }\propto {\lambda }^{\ast 2}{(R_{\mathrm{Edd}}M)}^{0.5}\propto {\lambda }^{\ast 2}{L}_{\mathrm{bol}}^{0.5}$, as was calculated for example by MacLeod et al. (2010). Our fits to the data can be expressed in terms of L_bol ∝ M × R_Edd as $t_{\mathrm{b}}\propto R_{\mathrm{Edd}}^{-D}{M}^{-C}\propto L_{\mathrm{bol}}^{-D}{M}^{D-C}=L_{\mathrm{bol}}^{0.35}{M}^{-0.2}$ or $t_{\mathrm{b}}\propto L_{\mathrm{bol}}^{0.55}$ $R_{\mathrm{Edd}}^{0.2}$, that is, 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 is 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 the orbital timescale and the 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 started from the flux emitted by a thin accretion disc (e.g., Eq. (23) in Spruit et al. 1995):

$$\begin{array}{c}\hfill \sigma T^4=F_{\mathrm{acc}}=\frac{3GM\dot{M}f}{8\pi R^3}\end{array}$$(8)

and assumed that the mass accretion rate, Ṁ, in cgs units corresponds to

$$\begin{array}{c}\hfill \dot{M}=\frac{L}{\eta c^2}=\frac{L_{\mathrm{Edd}}{R}_{\mathrm{Edd}}}{\eta c^2}=\frac{1.26\times {10}^{38}m{R}_{\mathrm{Edd}}}{\eta c^2}.\end{array}$$(9)

Replacing this value in Eq. (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

$$\begin{array}{cc}& \sigma T^4=\frac{3GMf}{8\pi R_{\mathrm{g}}^3{r}^3}\frac{1.26\times {10}^{38}m{R}_{\mathrm{Edd}}}{\eta c^2}\hfill \end{array}$$(10)

$$\begin{array}{cc}& \sigma T^4=\frac{3f1.26\times {10}^{38}m{R}_{\mathrm{Edd}}}{8\pi R_{\mathrm{g}}^2{r}^3\eta }=\frac{3f5.6\times {10}^{27}{R}_{\mathrm{Edd}}}{8\pi r^3\eta m}\hfill \end{array}$$(11)

$$\begin{array}{cc}& T(r)={\left(\frac{1.2\times {10}^{31}f}{\eta r^3}\frac{R_{\mathrm{Edd}}}{m}\right)}^{0.25}\mathrm{K},\hfill \end{array}$$(12)

where r = R/R_g, R_g = GM/c² = 1.5 × 10⁵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 $\text{d}A=2\pi r\text{d}r{R}_{\text{g}}^2$ at a given frequency, ν, is

$$\begin{array}{c}\hfill \mathrm{d}L=\mathrm{d}A\times \frac{4\pi h{\nu }^3}{c^2({\mathrm{e}}^{h\nu /kT(r)}-1)}.\end{array}$$(13)

These annular fluxes were used to calculate a light-weighted average radius as

$$\begin{array}{c}\hfill r_{\mathrm{lw}}=\frac{{\int }_{r_{\min }}^{r_{\mathrm{out}}}r\mathrm{d}L}{{\int }_{r_{\min }}^{r_{\mathrm{out}}}\mathrm{d}L}=\frac{{\int }_{r_{\min }}^{r_{\mathrm{out}}}{({\mathrm{e}}^{h\nu /kT(r)}-1)}^{-1}{r}^2\mathrm{d}r}{{\int }_{r_{\min }}^{r_{\mathrm{out}}}{({\mathrm{e}}^{h\nu /kT(r)}-1)}^{-1}r\mathrm{d}r}.\end{array}$$(14)

For these calculations we assumed $f=1-\sqrt{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, that is, $t_{\mathrm{orb}}=2\pi r_{\mathrm{lw}}^{3/2}{R}_{\mathrm{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 the 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 solid blue 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_{\mathrm{orb}}\propto {R_{\mathrm{Edd}}}^{0.35}$ relation, plotted as a 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.

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 REdd. The top panel shows the timescales as a function of M, colour-coded by REdd, together with a solid line displaying a torb ∝ M0.65 relation. The bottom panel shows the same model points as a function of REdd, colour-coded by mass, and the solid line displays a $t_{\mathrm{orb}}\propto {R_{\mathrm{Edd}}}^{0.35}$ relation.

As can be seen in Fig. 8, the orbital timescale at the light-weighted radius of 2900 Å, for the standard thin disc, and for M = 10⁸ M_⊙ and R_Edd = 0.1, is about ten 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 nine to 20. As has been 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 ten times the orbital timescale, independently of the accretion rate and most disc parameters (e.g., Spruit et al. 1995, and references therein).

We note however that this modelling with a standard thin disc produces a bolometric correction factor, k₅₁₀₀ = L_bol/L₅₁₀₀, 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₅₁₀₀ = 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₅₁₀₀ 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, that is, $f_{\mathrm{b}}\propto {R_{\mathrm{Edd}}}_{\mathrm{model}}^{-0.25}$, 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\propto R_{\mathrm{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 a larger amplitude of fluctuations. This scaling might help in constraining accretion disc models that incorporate prescriptions for variations in luminosity such as the inhomogeneous disc model with thermal perturbations in clumps proposed by Dexter & Agol (2011) and expanded upon 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 the 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 year) 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 is 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¹.

7.2. 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 anti-correlated in the sample. In the most extreme case, where objects in a sample have $M\propto R_{\mathrm{Edd}}^{-1}$, then, for example, f_b ∝ M^−0.6 $R_{\mathrm{Edd}}^{-0.3}$ would result in f_b ∝ M^−0.6M^+0.3 = M^−0.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 fom 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, that is, 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 the accretion rate. In all models, this relation results in either $P\times f\propto R_{\mathrm{Edd}}^{-0.75}$ or $P\times f\propto R_{\mathrm{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.

7.3. 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. Models with a flatter high-frequency slope of α_H = −2, similar to the DRW model, result in significantly worse fits. This last result is consistent with previous findings by Caplar et al. (2017) who used light curves 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 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 low-frequency 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 give bend frequencies with overlapping confidence ranges, and the models with α_H = −3 result 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 power law model, since steeper slopes fit better with bends at higher frequencies.

8. 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 four different timescales, in the range of 30–300 days in the rest-frame of the quasars (already presented in Arévalo et al. 2023). With these variance estimates, we constructed 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 light curves. 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 into 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, low- and 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⁸ M_⊙ and R_Edd = 0.1, range from ${96}_{-20}^{+9}$ days for the model with α_H = −3 and α_L = −1 to ${192}_{-10}^{+46}$ 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, give 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_{\mathrm{b}}\propto M^{-0.67\mathrm{to}-0.5}\times {R_{\mathrm{Edd}}}^{-0.41\mathrm{to}-0.12}$, where the ranges of the exponent values correspond to a power-spectral 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 is the opposite of 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 wavelength, 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 the 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_{\mathrm{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 is 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 the rest-frame wavelength but still weaker than that predicted by the thin disc model. Long-term, well-sampled light curves of many more objects, in multiple optical bands and to great 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.

Acknowledgments

The authors would like to thank Prof. H. Netzer for providing his model results. We also thank the anonymous referee for their useful comments and suggestions. The authors acknowledge financial support from the National Agency for Research and Development (ANID) grants: Millennium Science Initiative Program ICN12_12009 (PSS,LHG), and NCN19_058 (PA, PL, EL, SB); FONDECYT Regular 1201748 (PL); FONDECYT Postdoctorado 3200250 (PSS); Programa de Becas/Doctorado Nacional 21212344 (SB), 21200718 (EL) and 21222298 (PP), and from the Max-Planck Society through a Partner Group grant (PA, SB). Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant No. AST-2034437 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and UW.

References

Appendix A: Re-calculating R_Edd according to the standard thin disc model

According to the standard thin disc model, the surface temperature of the disc is a function of mass and R_Edd so that the luminosity at a fixed wavelength represents a different fraction of the total disc luminosity depending on these parameters. The bolometric correction used by Rakshit et al. (2020) to calculate the R_Edd in their catalogue has a fixed value, and therefore will over- or underestimate the R_Edd for colder and hotter discs, respectively, according to this model. In particular, for the redshift range used in our work, the bolometric luminosity was estimated from the monochromatic luminosity at 5100 Å and a bolometric correction factor, k₅₁₀₀ = 9.26, that is, R_Edd, cat = L₅₁₀₀ × k₅₁₀₀/L_Edd. In this section we use the subscript cat to refer to the R_Edd obtained from the catalogue, which is the one used throughout this paper, to distinguish it from R_Edd, model, that is, the Eddington ratio used as an input in the disc model.

The dependence of the bolometric correction factor, k₅₁₀₀, on mass, accretion rate, and black hole spin is described in detail by Netzer (2019). For a simple analysis, we fitted their results for a spin parameter, a = 0.7, kindly provided by Dr. Netzer. This value of a produces values of k_{5100, Netzer}, approximately halfway between a = 0.998 and a = 0, −1. With these we can convert the R_Edd values as they would appear in the catalogue of Rakshit et al. (2020) to the R_Edd values of the model as

$$\begin{array}{c}\hfill R_{\mathrm{Edd},\mathrm{model}}=\frac{L_{5100}\times k_{5100,\mathrm{Netzer}}}{L_{\mathrm{Edd}}}=R_{\mathrm{Edd},\mathrm{cat}}\frac{k_{5100,\mathrm{Netzer}}(M,R_{\mathrm{Edd},\mathrm{cat}})}{k_{5100}}\end{array}$$(A.1)

We note that for the fits we evaluated the model values of k_{5100, Netzer} as a function of R_Edd, cat, not of R_Edd, model, in order to replace the values of R_Edd, cat in the power spectral fitting functions. To do this, we multiplied the luminosity at 5100 Å (i.e. the x axis in the third panel in Figure 2 of Netzer 2019) by k₅₁₀₀/(1.26e38 × M) and fitted k_{5100, Netzer} as a function of this new variable. The best-fitting linear models to log(R_Edd, cat) vs log(k_{5100, Netzer}), for a spin parameter a = 0.7, were

$$\begin{array}{cc}\hfill log(k_{5100,\mathrm{Netzer}})=0.39log(R_{\mathrm{Edd},\mathrm{cat}})+1.3;& M={10}^9{M}_{\odot }\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill log(k_{5100,\mathrm{Netzer}})=0.43log(R_{\mathrm{Edd},\mathrm{cat}})+1.5;& M={10}^8{M}_{\odot }\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill log(k_{5100,\mathrm{Netzer}})=0.49log(R_{\mathrm{Edd},\mathrm{cat}})+2.1;& M={10}^7{M}_{\odot }.\hfill \end{array}$$(A.4)

To simplify the analysis, we took the relation for M = 10⁸M⊙ for the slope and scaled the normalisation by the mass. Taking the values of k_{5100, Netzer} at a single value of R_Edd, cat = 0.1, we obtained the relation between the bolometri correction and mass: log(k_{5100, Netzer}) = − 0.38log(M/M₈)+1.1. Putting both together, we obtained the approximate relation

$$\begin{array}{c}\hfill log(k_{5100,\mathrm{Netzer}})=0.4log(R_{\mathrm{Edd},\mathrm{cat}})-0.4log(M/M_8)+1.5,\end{array}$$(A.5)

and therefore the relation between R_Edd, model and R_Edd, cat is

$$\begin{array}{cc}\hfill log(R_{\mathrm{Edd},\mathrm{cat}})& =\frac{log(R_{\mathrm{Edd},\mathrm{model}})+0.4log(M/M_8)-1.5+log(k_{5100})}{1.4}\hfill \\ \hfill log(R_{\mathrm{Edd},\mathrm{cat}})& =0.7log(R_{\mathrm{Edd},\mathrm{model}})+0.3log(M/M_8)-0.38\hfill \end{array}$$

Recalling the dependence of f_b on these parameters and replacing R_Edd, cat gives

$$\begin{array}{cc}\hfill log{f}_{\mathrm{b}}& =logB+C\times log(M/M_8)+D\times log(R_{\mathrm{Edd},\mathrm{cat}}/0.1)\hfill \\ \hfill log{f}_{\mathrm{b}}& =logB+C\times log(M/M_8)+D\times 0.7log(R_{\mathrm{Edd},\mathrm{model}})\hfill \\ \hfill & +D\times (0.3log(M/M_8)-0.38-1)\hfill \\ \hfill log{f}_{\mathrm{b}}& =logB+(C+0.3D)\times log(M/M_8)+0.7D\times log{R}_{\mathrm{Edd},\mathrm{model}}\hfill \\ \hfill & -1.38D\hfill \\ \hfill log{f}_{\mathrm{b}}& =logB-0.655\times log(M/M_8)-0.245\times log{R}_{\mathrm{Edd},\mathrm{model}}\hfill \\ \hfill & +0.48,\hfill \end{array}$$

where in the last equation we have replaced C and D by their best-fitting values, -0.55 and -0.35, respectively. Therefore, including this correction of the accretion rate, from the constant bolometric correction used in the catalogue to the value predicted for the measured L₅₁₀₀ and mass by the thin disc model, results in an even stronger dependence of f_b on M and a weaker but still negative dependence of f_b on R_Edd.

Appendix B: Figures of merit as a function of the power spectral model parameters

Below we show the dependence on the figure of merit, fom/N, as a function of the four fitted parameters, for the five better-fitting models.

Fig. B.1. Figure of merit values as a function of the bend frequency for M = 108 M⊙ and REdd = 0.1, B, for the better-fitting models with different power law slopes (αL and αH, as noted in the title of each panel). Blue dots show the minimum, fom/N, for each value of B, with all other parameters free, corresponding to the second iteration, with a finer grid in the parameters. The horizontal black line marks the lowest fom/N achieved and the grey line marks the level of min(fom)/N + Δ, from which the error bounds were estimated. The value of Δ was obtained for each model, as is described in Sect. 6.1.

Fig. B.2. As Fig.B.1 but for parameter C, which determines the dependence of the bend frequency on M, that is, fb∝M C.

Fig. B.3. As Fig.B.1 but for parameter D, which determines the dependence of the bend frequency on REdd, that is, $f_{\text{b}}\propto R_{\text{Edd}}^D$.

Fig. B.4. As Fig.B.1 but for parameter F, which determines the dependence of the normalization of the power spectrum on REdd, that is, A∝ $R_{\text{Edd}}^F$.

Fig. B.5. Minimum fom for given values of D, as a function of F. Both parameters relate to the dependence of the power spectrum on REdd and are therefore partially degenerate, which produces the flat profiles of fom vs D. It is clear, however, that extreme values of D produce worse overall fits for any values of F.

All Tables

All Figures

Fig. 1. Low-resolution power spectra for different values of REdd within a narrow range in black hole mass (left) and different black hole masses but similar values of REdd (right). The markers represent the median variance (i.e., power × frequency) of all quasars falling in each M − REdd 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 the location of the bend in the power spectrum, as described in Sect. 4) and the accretion rate modifies mostly the normalisation, although some change in the slope is apparent. Frequencies are in units of day−1.

In the text

Fig. 2. Low-resolution power spectra (plotted as power × frequency) for different mass bins in different panels and different accretion rate levels in different colours. The markers represent the median variance of all quasars falling in each M − REdd bin and temporal frequency, and the error bars represent the rms of the median for each bin obtained by bootstrapping the samples. Frequencies are in units of day−1.

In the text

	Fig. 3. Similar to Fig. 2 but arranged in separate panels for different REdd bins, and with different colours used to differentiate by mass. In most REdd 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.
In the text

Fig. 4. Power spectra (plotted as power × frequency) of all the sources, binned by mass and Eddington ratio. The variances were calculated for timescales of 30, 75, 150, and 300 days, and in this plot the corresponding frequencies, f = 1/t, were multiplied by (M/M8)0.5 to put them on a common x-axis scaled to the bend timescale approximately, as was measured by Burke et al. (2021). The colour represents the different bins in REdd. The solid black 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 108 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 to rise with frequency.) Frequencies are in units of day−1.

In the text

Fig. 5. The large dots, colour-coded by REdd, 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 REdd of each bin and the best-fitting scaling parameters, C = −0.55, D = −0.35, and 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.

In the text

	Fig. 6. Residuals of the model and data plotted in Fig. 5, colour-coded by mass. Frequencies are in units of day−1.
In the text

	Fig. 7. Distributions of values of fom/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 high-frequency slopes of their power spectra and these are written in the legend.
In the text

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 REdd. The top panel shows the timescales as a function of M, colour-coded by REdd, together with a solid line displaying a torb ∝ M0.65 relation. The bottom panel shows the same model points as a function of REdd, colour-coded by mass, and the solid line displays a $t_{\mathrm{orb}}\propto {R_{\mathrm{Edd}}}^{0.35}$ relation.

In the text

Fig. B.1. Figure of merit values as a function of the bend frequency for M = 108 M⊙ and REdd = 0.1, B, for the better-fitting models with different power law slopes (αL and αH, as noted in the title of each panel). Blue dots show the minimum, fom/N, for each value of B, with all other parameters free, corresponding to the second iteration, with a finer grid in the parameters. The horizontal black line marks the lowest fom/N achieved and the grey line marks the level of min(fom)/N + Δ, from which the error bounds were estimated. The value of Δ was obtained for each model, as is described in Sect. 6.1.

In the text

	Fig. B.2. As Fig.B.1 but for parameter C, which determines the dependence of the bend frequency on M, that is, fb∝M C.
In the text

	Fig. B.3. As Fig.B.1 but for parameter D, which determines the dependence of the bend frequency on REdd, that is, $f_{\text{b}}\propto R_{\text{Edd}}^D$.
In the text

	Fig. B.4. As Fig.B.1 but for parameter F, which determines the dependence of the normalization of the power spectrum on REdd, that is, A∝ $R_{\text{Edd}}^F$.
In the text

	Fig. B.5. Minimum fom for given values of D, as a function of F. Both parameters relate to the dependence of the power spectrum on REdd and are therefore partially degenerate, which produces the flat profiles of fom vs D. It is clear, however, that extreme values of D produce worse overall fits for any values of F.
In the text