© The Authors 2024

1. Introduction

Variability is one of the most striking features of quasars and has been known since they were discovered (Matthews & Sandage 1963). The last 60 yr of observations have shown that it is observed at all wavelengths, both in the continuum and the emission lines, with timescales ranging from minutes to years and a completely stochastic and aperiodic behaviour (e.g. Ulrich et al. 1997; Giveon et al. 1999; Hawkins 2002; Vanden Berk et al. 2004; Wilhite et al. 2005, 2008). In effect, X-ray and UV/optical power spectral densities (PSDs, or power spectra) extracted from light curves of quasars show no periodic features and a typical red-noise trend, usually modelled with one or more power laws (e.g. Collier & Peterson 2001; Uttley et al. 2002; McHardy et al. 2004).

Quasars are the most luminous among non-obscured active galactic nuclei (AGN). It is generally accepted that most of their UV/optical light comes from thermal emission attributed to an accretion disc, while X-rays are produced by inverse-Compton scattered photons from a hot electron corona surrounding the central supermassive black hole (SMBH; see Netzer 2015 and references therein). However, the exact interplay between the two phenomena, as well as the origin of the variability itself, is still debated. Some authors support the idea of a thermal origin with propagation of instabilities through the accretion disc and/or changes in the accretion rate (e.g. Neustadt & Kochanek 2022 and references therein). Other evidences favour the idea of optical variability due to X-ray reprocessing (e.g. Kammoun et al. 2021; Panagiotou et al. 2022), but a combination of multiple effects is also not excluded.

The last two decades have seen an increase in both the quality and quantity of multi-epoch UV/optical data from large ground-based surveys, such as the Optical Gravitational Lensing Experiment (OGLE; Kelly et al. 2009; Kozlowski et al. 2010; Zu et al. 2013), the Sloan Digital Sky Survey (SDSS) Stripe-82 region (Sesar et al. 2006; MacLeod et al. 2010, 2012; Guo et al. 2017; Yu et al. 2022), the Palomar-QUEST Survey (Bauer et al. 2009), the Panoramic Survey Telescope And Rapid Response System survey (PanSTARRS-1; Morganson et al. 2014; Simm et al. 2016), LaSilla-Quest (Cartier et al. 2015; Sánchez-Sáez et al. 2018), the VST SUDARE-VOICE (Falocco et al. 2015; De Cicco et al. 2015, 2019, 2022; Poulain et al. 2020), the Palomar Transient Factory (PTF/iPTF surveys; Caplar et al. 2017), the Hyper Suprime-Cam Subaru Strategic Program (Kimura et al. 2020), the Catalina Real-time Transient Survey (Rakshit & Stalin 2017; Graham et al. 2020; Tachibana et al. 2020; Laurenti et al. 2020), the NASA Asteroid Terrestrial-impact Last Alert System (ATLAS; Tang et al. 2023), and the Zwicky Transient Facility (ZTF; López-Navas et al. 2023; Arévalo et al. 2023). Space missions such as Kepler and the Transiting Exoplanet Survey Satellite (TESS) also provided data with a high cadence (e.g. Mushotzky et al. 2011; Edelson et al. 2014; Kasliwal et al. 2015; Treiber et al. 2023).

Many efforts in combining data from different surveys enabled the variability on long timescales to be constrained. Rumbaugh et al. (2018) joined SDSS and Dark Energy Survey (DES) light curves, Suberlak et al. (2021) used SDSS Stripe-82 and PanSTARRS-1 data, while Stone et al. (2022) studied a sample of quasars observed with SDSS, PanSTARRS-1, DES, and a dedicated monitoring with the DECam on the CTIO-4m Blanco telescope, over 20 yr. The near future will be even more promising. The Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2009; Ivezić et al. 2019), expected to start in 2025, will observe the entire southern sky every 3–4 nights for 10 yr, across the ugrizy filters. This new facility will detect tens of millions of AGN (also using variability as a selection tool, see Savić et al. 2023), allowing a thorough analysis of light curves, including fast variability (e.g. Raiteri et al. 2022) and precise time-delay measurements (Kovačević et al. 2022; Czerny et al. 2023). Tests to define the final cadence requirements are still ongoing at the time of writing (e.g. Brandt et al. 2018; Kovačević et al. 2021; Sheng et al. 2022).

All surveys mentioned above allow us to study the statistical properties of optical variability and its correlations with the physical parameters of quasars. Ensemble analyses have shown that rest-frame variability amplitude increases with time (i.e. the red-noise trend of power spectra). Variability is also observed to be anti-correlated with luminosity (and/or the accretion rate) at all timescales, and many works also find an anti-correlation with rest-frame wavelength. Dependencies on the black hole mass, M_BH, are less certain, as well as the overall shape of the PSD (i.e. spectral indices and characteristic timescales). A partial list comparing the observed correlations between variability and quasar properties is shown in Table 4 of Suberlak et al. (2021), along with the different methods adopted to measure variability.

Optical light curves are often modelled using a damped random walk (DRW; Kelly et al. 2009), a first-order continuous-time autoregressive moving-average (CARMA) process. DRW is a statistical model which describes AGN variability via an exponential covariance matrix with two parameters: a variability amplitude, σ, and a characteristic damping timescale, τ. Such a model predicts a power-law PSD with a spectral index of −2, flattening to 0 (i.e. white noise) for times ≫τ. However, DRW processes are sometimes studied in time-domain via the structure function (SF), defined as the root mean square (rms) magnitude difference as a function of the time difference Δt between observation pairs (Kozłowski 2016a and references therein). DRW (or associated SF parameters) can be then studied as a function of physical properties of quasars (e.g. MacLeod et al. 2010; Zu et al. 2013; Kasliwal et al. 2015; Suberlak et al. 2021). Although there are known biases in the recovery of reliable damping timescales from observed light curves, as well as uncertainties on the SF itself due to sampling effects (e.g. Emmanoulopoulos et al. 2010; Kozłowski 2017), many attempts have been made to relate the DRW τ to any physical timescale (see e.g. Burke et al. 2021). Regardless of the proper estimate of τ, variability amplitude predicted by the DRW model seems to be consistent with light curves length of the order of a few years, but deviations have been observed on both longer and shorter timescales (e.g. Mushotzky et al. 2011; Guo et al. 2017). Possible steps forward may invoke the usage of higher-order CARMA processes, such as the damped harmonic oscillator (DHO; Yu et al. 2022), indirect PSD estimation through CARMA modelling (Kelly et al. 2014), or unsupervised machine learning analysis of time series (Tachibana et al. 2020).

In this work, we propose a model-independent variability analysis of ∼9000 spectroscopically confirmed quasars from the SDSS Stripe-82 region through direct PSD measurement. We show how to account for irregular sampling and take advantage of the large dataset to execute ensemble analysis and constrain the PSD shape as a function of redshift, luminosity, black hole mass, accretion rate, and wavelength. We provide scaling relations for the PSD on timescales of ∼[250, 1500] days (rest-frame) and present a framework which could be easily applied to other large multi-epoch surveys. We will discuss the PSD shape on shorter and longer timescales in a future work (Petrecca et al., in prep.).

This paper is organized as follows. In Sect. 2 we present the data used for our work. Section 3 describes the variability measurement method through PSD, while Sect. 4 is dedicated to the analysis of the correlations between variability and physical properties. Finally, we summarize our results and discuss possible interpretation and future perspectives in Sect. 5.

2. The sample

Our analysis requires a large sample of AGN with very well-known properties and well-sampled light curves. Among all the surveys listed in the Introduction, the Stripe-82 region of the SDSS (York et al. 2000) suits our purposes perfectly. This area consists of 290 deg² with right ascension from 20:00h to 4:00h and declination from −1.26° to +1.26°, imaged multiple times by the SDSS from 2000 to 2008. There are ∼60 observations in 5 filters u, g, r, i, z (Fukugita et al. 1996; Gunn et al. 1998; Smith et al. 2002) grouped in yearly seasons about 2–3 months long, albeit most of the data are in the last 6 seasons. Images in the different filters are nearly simultaneous, separated by 71.72 s of drift scan time.

Our sample of AGN consists of 9275 well-characterized, spectroscopically confirmed quasars selected by MacLeod et al. (2010, 2012), who provide tables with light curves¹. We retrieved virial black hole masses and bolometric luminosities for 9148 quasars from the SDSS DR7 quasar catalogue by Shen et al. (2011). Masses are estimated from emission line widths (H_β for z < 0.7, Mg II for 0.7 ≤ z < 1.9, and C IV for z ≥ 1.9). The nominal mean uncertainty on fiducial mass estimates is estimated to be ∼0.2 dex, although systematics from reverberation mapping could be as high as ∼0.5 (Krolik 2001). Despite such uncertainties, the wide range of available masses still allows us to study whether variability depends on mass, as explained in Sect. 4. Bolometric luminosities are derived from the continuum using bolometric corrections derived by Richards et al. (2006) from the composite spectral energy distribution (BC₅₁₀₀ = 9.26 for z < 0.7, BC₃₀₀₀ = 5.15 for 0.7 ≤ z < 1.9, and BC₁₃₅₀ = 3.81 for z ≥ 1.9). Average uncertainties on bolometric luminosities are of the order of 0.1 dex. We discard from the sample 125 sources with no valid estimate of black hole mass. Figure 1 shows the distribution of the quasar magnitudes in the ugriz bands. We restricted the sample to magnitudes 17.5 < mag_AB < 21.5 in all bands, as there are few quasars with mag_AB < 17.5, while at mag_AB > 21.5 there are no stars available for the necessary computation of the Poisson noise (see Sect. 3).

Fig. 1. Distribution of observed quasar magnitudes in the five SDSS filters. Dashed lines enclose the magnitude range of sources used in this analysis.

The sample, after the magnitude cuts, consists of 8119 quasars with known black hole mass, M_BH, bolometric luminosity, L_bol, and redshift, z. Their distribution is shown in Fig. 2 (left, middle and right panels, respectively). We cover a broad range of masses (from ∼10^7.5 up to ∼10¹⁰ M_⊙), and luminosities (10⁴⁵–10⁴⁷ ergs s⁻¹), which are typical of the quasars. Most of the sources have redshifts 0.5 ≤ z ≤ 2.5, although a tail at larger redshifts is also apparent in the figure. The large number of quasars over a broad range of M_BH and L_bol values will give us the opportunity to study if and how the variability properties of quasars (quantified by the PSD analysis) depend on each of the two parameters. The distribution of a large range of redshift will also allow us to investigate if the quasar optical variability properties evolve with time.

Fig. 2. Distribution of black hole mass, bolometric luminosity, and redshift for the quasars in our sample.

In addition, the availability of light curves across five bands will enable us to compare variability properties over different rest-frame continuum wavelengths. Figure 3 shows the distribution of the rest-frame wavelength of each light curve in all filters, of all quasars, computed as λ_rest = λ_obs/(1 + z), where λ_obs,u = 3520 Å, λ_obs,g = 4800 Å, λ_obs,r = 6250 Å, λ_obs,i = 7690 Å, and λ_obs,z = 9110 Å are the SDSS photometric system central wavelengths. The five left-most vertical dashed lines in the figure indicate the wavelength where the u, g, r, i and z rest-frame wavelength distribution peak is located (λ_rest, i = 1300, 1800, 2300, 3000, and 3400 Å, where i = 1, 2, 3, 4 and 5). The other two vertical lines indicate two additional rest-frame wavelengths, namely λ_rest, 6 = 4000 and λ_rest, 7 = 5000 Å. When we compute the power spectrum of quasars in various rest-frame bands (see Sect. 4), we consider all light curves, irrespective of the filter, whose rest-frame wavelength is within |Δλ|/λ_rest = 0.2, around the λ_rest wavelengths indicated by the vertical lines in Fig. 3. As for each quasar we have five distinct light curves in ugriz, we only use the one whose λ_rest is the closest to the centre of the selected wavelength bin. The number of quasars within the λ_rest, 6, and in particular within the λ_rest, 7 band, is smaller than the number of quasars in the other rest-frame wavelength bins. We considered those longer wavelengths, so that we can study the energy dependence of the power spectrum over of range of λ_rest, 7/λ_rest, 1 ∼ 4 in energy separation.

Fig. 3. Distribution of rest-frame wavelengths (λrest) of all light curves, for all quasars in the sample. The vertical dashed lines indicate the centres of the λrest bins at which we computed the PSD of quasars.

As this work makes use of broad-band photometry data, we expect some contamination from the broad line region (BLR) to the accretion disc continuum (mostly broad lines, blended iron multiplets, and the Balmer pseudo-continuum). Although spectroscopic studies of AGN from the SDSS show contamination levels up to ∼20% in the most severe cases, the mean contribution over a broad wavelength region is typically less than 10% (Vanden Berk et al. 2001; Calderone et al. 2017). Besides, the light-crossing time for the BLR is one order of magnitude greater than for the UV/optical disc, on average (Cackett et al. 2021). For this reason, we do not expect a major contribution from line variability to our analysis. Host-galaxy contamination is also negligible, as it starts to be significant for low-luminosity quasars with z ≲ 0.5 and at λ_rest above 5000 Å (Shen et al. 2011), while we mostly sample higher-luminosity and higher-redshift quasars at bluer wavelengths.

3. The variability analysis method

3.1. The SDSS light curves

Stripe-82 was observed for almost 10 yr, with a denser sampling cadence as the survey progressed, yielding very similar light curves in terms of observed sampling for all quasars in our sample. As an example, Fig. 4 shows the g-band light curves of two quasars, which have similar M_BH, L_bol and redshift. In the later observing seasons, sources were observed approximately 10 times in ∼3 months. SDSS magnitudes have been corrected for galactic extinction (coefficients from Schlegel et al. 1998) and converted into fluxes in units of nJy, assuming an AB photometric system with a zero point of 3631 Jy and following a standard procedure for error propagation.

Fig. 4. Light curves of two quasars in g-band with a similar BH mass (8.8 < log(MBH) < 9.0), luminosity (46.1 < log(Lbol) < 46.3), and redshift (1.1 < z < 1.3). Red points indicate the binned light curve.

We computed the (straight) mean flux per observing season (red dots in Fig. 4). The observing seasons are equally spaced by 1 yr in the observer frame. We verified that the deviations w.r.t. the mean epoch of all observations are of the order of just a few sigma. Errors on the binned fluxes are computed using the standard error on the mean $\sigma /\sqrt{n}$, where σ is the standard deviation and n is the number of points in each season. In this way we obtain evenly sampled light curves, with at least seven points. This is the case for almost all quasars in our sample, given the similarity of Stripe-82 observations in terms of cadence for all sources.

The second season data point at around ∼52200 MJD is missing in Fig. 4 as well as in many quasar light curves. In addition, there are a few more missing points among the last six observations for 77 sources. For that reason, we rejected them and kept N_Q = 8042 sources in the sample, for which we can compute binned light curves, with no missing points in the last six observations, like the light curves shown in Fig. 4.

3.2. Power spectrum derivation

We used the last six seasons of the binned light curves, where points are evenly sampled, and we computed the periodogram of each one, that is,

$$\begin{array}{c}\hfill I_N(f_j)=\frac{2\mathrm{\Delta }t}{N{\overline{x}}^2}[{\displaystyle \sum _{i=1}^N}[x(t_i)-\overline{x}]e^{-i2\pi f_j{t}_i}{]}^2,\end{array}$$(1)

which we accept as an estimate of the intrinsic power spectrum. In the equation above, N = 6 is the number of points in the binned light curves used for the analysis, x(t_i) are the binned flux values, and $\overline{x}$ is the mean flux. We computed the rest-frame power spectrum of each source, that is we assume Δt = 365/(1 + z) is the size of the time window, and all times t_i are divided by (1 + z). Frequencies are defined as f_j = j/(NΔt), with j = 1, 2 and 3 (by using six points in the binned light curves, we end up with three frequencies).

The periodogram is the most common estimator of the PSD function of a variable process. The PSD is the Fourier transform of the autocovariance function of the process, (which holds all information about the important variability properties of the process, including the presence of intrinsic timescales) and it is frequently used in many research fields, including astronomy (e.g. Press 1978; Priestley 1981; van der Klis 1988; and references therein). Knowledge of the intrinsic PSD of a variable source can be very important in constraining various physical (or statistical) models that have been proposed to interpret the observed variations of an AGN in the UV/optical bands. Although power spectral analysis is a well known and used tool for time series interpretation, it may suffer from biases depending on the sampling pattern and duration; thus we further validated the robustness of our method via simulations, as shown in Appendix B.

Observed PSDs should also be affected by the experimental Poisson noise. In principle, we could predict the resulting Poisson noise level knowing the photometric accuracy of the data. This is not true in our case due to various reasons. For example, as we mentioned above, we need to retrieve the original uncertainties transforming back the magnitudes to fluxes, but this process is usually based on a first-order approximation. In addition, there could be residual uncertainties (e.g. calibration) that are not properly taken into account when estimating the error on magnitudes. We therefore followed a different approach, and we computed the Poisson noise level using the light curves of non-variable stars that have the same magnitudes as the quasars in our sample, as suggested by, e.g. Kozłowski (2016b), and as we explain in detail in Appendix A.

Figure 5 show the rest-frame g-band PSDs, using the binned light curves plotted in Fig. 4. Since we consider only the last six points in the binned light curves, we compute the power spectrum in 3 frequencies, over half a decade in frequency space. Solid lines indicate the best-fit straight line to the PSDs (in the log–log space). The line appears to be a reasonable fit to the data, but we cannot compute errors on the best-fit parameters and we cannot judge the goodness of fit, since we do not know the error of each PSD measurement, as it depends on the (unknown) intrinsic PSD (e.g. Papadakis & Lawrence 1993). However, we can do much better than just estimating the PSD of each source in the sample, over three frequencies, by taking advantage of the large number of quasars in the sample. For example, we can group the quasars into narrow black hole (BH) mass, luminosity and redshift bins, compute the (noise-corrected) periodogram of each one of them, and then average all the periodograms in order to compute the ensemble power spectrum of all quasars in each [M_BH, L_bol, z] bin.

Fig. 5. Noise-corrected PSD of the two light curves shown in Fig. 4. The solid lines indicate the best-fit line to the data.

We show the logarithm of the (noise-corrected) periodogram of all quasars which belong to the BH mass, luminosity and redshift bin with: [8.8 < log(M_BH) < 9.0], [46.1 < log(L_bol) < 46.3], and [1.2 < z < 1.4] in Fig. 6. The periodograms are calculated in the rest-frame of each source and, since we consider a redshift bin with a small width of 0.2, the difference in rest-frame frequencies is small. We therefore compute the mean of all the periodograms in each frequency bin, and the standard error on the mean, and we accept this as representative of the average power spectrum of quasars with the chosen BH mass and luminosity. As long as there are at least 50 objects in each bin (and hence 50 periodogram estimates at each frequency), the mean periodogram should approximately follow a Gaussian distribution (Papadakis & Lawrence 1993). In practice, we compute the logarithm of the mean periodogram, and its error (following the standard error propagation formula), as it is more convenient to fit straight lines to sampled power-spectra, when the intrinsic power spectrum follows a power-law form (as is usually the case with the quasar power spectra).

Fig. 6. Logarithm of the noise-corrected periodogram of all quasars with [8.8 < log(MBH) < 9.0], [46.1 < log(Lbol) < 46.3], and [1.2 < z < 1.4]. Small blue dots indicate single PSD estimates, while black points show the logarithm of the mean PSD (with its error), and the red line is the best-fit linear model to the ensemble PSD.

We fitted the mean PSD in Fig. 6 with a straight line, in log–log space, i.e.

$$\begin{array}{c}\hfill {log}_{10}[\mathrm{PSD}(\nu )]={log}_{10}[{\mathrm{PSD}}_{\mathrm{amp}}]+{\mathrm{PSD}}_{\mathrm{slope}}[log(\nu )+2.6],\end{array}$$(2)

where PSD_amp is the power spectrum amplitude at ν = 10^−2.6 day⁻¹, i.e. ∼400 days. We chose to normalize the PSD at this frequency, as it is always within the range of sampled frequencies for our sample. In this way the error of the best-fit parameters is minimized. The solid red line in Fig. 6 shows the best-fit line to the mean PSD (χ² = 0.05 for 1 degree of freedom/d.o.f.).

4. Power spectral analysis results

As we already mentioned, our main objective is to study the optical power spectrum of quasars, using the light curves of 8042 sources in Stripe-82 in five bands ugriz. The main physical parameters of AGN are the BH mass, M_BH, and the accretion rate (derivable from the bolometric luminosity through the Eddington ratio, λ_Edd)², which determine the emitted energy spectrum. Thus, irrespective of the mechanism that drives the optical variability, it is plausible that the variability properties may also depend on these physical parameters or, perhaps, remain constant in all quasars. In addition, given the fact that quasars are powered by accretion of gas onto the BH, we expect the energy released through the accretion process, as well as the disc temperature, to increase inwards. This implies that the outer parts of the disc will be able to emit only longer wavelength photons. As the variability mechanism may depend on the local disc properties, we may expect the variability properties to also depend on radius, and hence, on the energy of the emitted photons (or their wavelength, i.e. λ_rest).

Based on the discussion above, we grouped the quasars into bins of BH mass and luminosity-accretion rate, we computed their average PSD (as explained in Sect. 3.2) in various rest-frame wavelengths, we fitted it with the model defined by Eq. (2), we determined PSD_amp, and PSD_slope, and we investigated whether PSD_amp and PSD_slope depend on M_BH, λ_Edd, and λ_rest, or if they are simply the same in all quasars, irrespective of energy and the quasar physical properties. Other relevant properties of AGN variability might be the line-of-sight orientation, the radio-loudness, the presence of a jet, and the BH spin, but these are much harder to measure. However, our sample consists only of quasars, i.e. with likely similar orientation, while significantly different variability amplitudes are only expected for extreme radio-loud sources (< 1% of the sample, see Table 2 in MacLeod et al. 2010).

We define BH mass bins, of size 0.2 dex, in the range 8.2 < log(M_BH) < 9.6, and L_bol bins of the same size, in the range 45.3 < log(L_bol) < 47.0. The lower/upper limits as well as the width of the bins are determined by the fact that, on the one hand, we need large bins with as many objects as possible but, on the other hand, objects in each bin should also have (almost) the same M_BH and L_bol. The first constraint requires large bin sizes, while the second implies the choice of small bin sizes. Given the uncertainty on the M_BH measurements (from the width of emission lines measured in optical spectra and systematics from reverberation mapping studies) and on estimating L_bol (from bolometric correction factors), a choice of a bin size of 0.2 may include some objects with similar BH mass and luminosity (and hence λ_Edd) in adjacent bins. Nonetheless, the wide range of masses and luminosities allows us to define many different bins and study any PSD trend with physical properties of quasars. We also define rest-frame wavelength (i.e. rest-frame energy) bins. The centres of those bins are indicated by the vertical lines plotted in Fig. 3. As we explained in Sect. 2, the width of each wavelength bin is such that, |Δλ_rest|/λ_rest = 0.2.

We compute the PSD of each single quasar, and then the mean PSD of all quasars in each [M_BH, λ_Edd(or L_bol), λ_rest] bin (as we explained in Sect. 3.2,) and we fitted it in the log–log space using the straight line model defined by Eq. (2). The number of quasars is larger than 50 in all bins (to ensure proper χ² statistics when assessing the goodness of fit to the PSDs). We found that a simple line can fit the mean quasar PSDs, in all bins. We discuss below how the best-fit PSD amplitude and slope depend on the various parameters. But before doing this, we first investigate the possibility that the variability properties may depend on redshift, which would imply that the variability mechanism may evolve with time.

4.1. Dependence of PSD on redshift

In order to investigate whether the PSD varies with redshift, we considered quasars with [8.8 < log(M_BH) < 9.0] and [46.1 < log(L_bol) < 46.3]. This BH mass and luminosity bin is the one with the largest number of quasars (see the peak of the BH mass and luminosity distribution in Fig. 1) and has a slightly larger bin size of 0.3 to further maximize the number of sources at different redshift. For each λ_rest, we divided the quasars in this [M_BH, L_bol] bin into various redshift bins, and we computed the mean PSD in each one of them (as described in Sect. 3.2). The results in the case when we use light curves with λ_rest = 3000 Å are shown in Fig. 7. The redshift bins were chosen such that each bin included at least 50 objects (to ensure Gaussian statistics, as we explain in Sect. 3.2). The average redshift in each bin is listed in the legend in Fig. 7, and their widths range from Δz = 0.05 to ∼0.2 − 0.3, in a few cases.

Fig. 7. Ensemble PSD of quasars with 8.8 < log(MBH) < 9.1 and 46.0 < log(Lbol) < 46.3 at λrest = 3000 Å in ten redshift bins (different symbols and colours). The average redshift for each bin is reported in the legend.

We fitted all the PSDs plotted in Fig. 7 at all redshifts, from z ∼ 0.9 − 2, with the model defined by Eq. (2). The solid line shows the resulting best-fit with a straight line. The best-fit χ² value is χ² = 47.4 for 28 d.o.f., and the resulting p_null³ is 0.012. Formally speaking, this means that we cannot reject (at the 99% confidence level) the null hypothesis that all points are drawn from the same PSD, and hence we can assume that the PSD is independent of redshift. In any case, even if there were some residual dependence on redshift, it would not be large enough to explain other strong trends which we report below.

Although we just report the results for λ_rest = 3000 Å as an example, we computed power spectra using light curves at all the rest frame wavelengths defined in Sect. 2 for the same M_BH and L_bol bin defined above, obtaining the same results. Furthermore, we repeated the same exercise considering quasars in smaller/larger BH mass/luminosity bins, so that we could sample smaller and larger redshift ranges, from z ∼ 0.5 to 2.5. In all cases, we found that the power spectra of quasars at different redshift bins are consistent with the same power-law model. This implies that we do not find significant evidence for the dependence of the quasar optical variability power spectrum on redshift. For the rest of the analysis, we consider all quasars with redshift from 0.5 to 2.5, when we compute their mean PSD in various mass and luminosity bins.

4.2. Dependence of PSD on L_bol

We proceed to investigate the dependence of the power spectrum on the bolometric luminosity, L_bol. Figure 8 shows the mean PSD of all quasars with [8.8 < log(M_BH) < 9.0] (i.e. the BH mass bin with the largest number of quasars, see left panel in Fig. 2), in six luminosity bins from log(L_bol) = 45.3 up to 46.5. The luminosity bin size is Δlog(L_bol) = 0.2, and we have used the λ_rest = 3000 Å light curves to compute the PSDs. The solid lines show the best-fit model lines (as defined by Eq. (2)). The model fits the data rather well. Since there are just three PSD points in each luminosity bin, there are some issues with the resulting error on the best-fit model parameter values. In some cases, the three PSD points just happen to be located exactly on the best-fit line (i.e. the red points in specific case of Fig. 8, where log(L_bol) = 45.8), while in other cases the points are spread around the best fit (e.g. the purple points in the same panel, with log(L_bol) = 45.4). Consequently, the error on the best-fit parameters varies considerably. Given the small number of data points in each PSD, the error from the best model-fits to the individual PSDs may not be representative of the true error of the best-fit parameters. For that reason, in all cases below, we consider the mean of all the best-fit parameter errors as a more representative error of the respective best-fit parameter.

Fig. 8. Ensemble PSDs of quasars with 8.8 < log(MBH) < 9.0 and all redshifts between 0.5 and 2.5 in different log(Lbol) bins. The legend reports the average luminosity in each bin of width 0.2.

Figure 8 shows that the power spectrum depends strongly on the source luminosity. The power spectra at different luminosity bins appear to have similar slopes but different amplitudes (we discuss possible mild dependencies of the slopes later in the text). It appears that the PSD amplitude increases with decreasing luminosity. This result is in agreement with a trend that has been known for a long time in many quasars, at all wavebands: more luminous sources are less variable (e.g. Barr & Mushotzky 1986; Cristiani et al. 1996; Papadakis & Lawrence 1993).

4.3. Dependence of PSD on M_BH

Similarly, we considered all quasars in the luminosity bin 46.1 < log(L_bol) < 46.3 (i.e. the L_bol bin with the largest number of quasars, see middle panel in Fig. 2) and we computed the mean PSD in narrow BH mass bins. The mass bin size is Δlog(M_BH) = 0.2, and we used the λ_rest = 3000 Å light curves to compute the PSDs, as before. The resulting PSDs are plotted in Fig. 9. As before with the PSDs in different luminosity bins, the PSDs in the different BH mass bins have similar slopes, but their amplitude appear to depend on M_BH, although not as strongly as the dependence on L_bol.

Fig. 9. Average PSDs of quasars with 46.1 < log(Lbol) < 46.3 and all redshifts between 0.5 and 2.5 in different log(MBH) bins. The legend reports the average mass in each bin of width 0.2.

Figure 9 shows that the variability amplitude increases with increasing BH mass. This is an unexpected result, as it is difficult to understand how a more massive BH will be more variable than a smaller one. However, since all quasars in the various mass bins have the same L_bol, they should have different accretion rates. Assuming that λ_Edd ∝ L_bol/L_Edd ∝ L_bol/M_BH, then larger BH mass objects in this case should correspond to objects with smaller accretion rates. Therefore, the apparent trend of PSD_amp increasing with increasing M_BH may in effect be due to the fact that power spectrum amplitude increases with decreasing accretion rate.

In fact, the same trend may be apparent in Fig. 8 as well. As luminosity decreases for quasars with the same M_BH, the accretion rate will also decrease. Therefore, the trend of PSD_amp increase with decreasing L_bol may be due to the fact that, in reality, PSD_amp increases with decreasing λ_Edd. But, when comparing Figs. 8 and 9, it seems that the PSD amplitude depends more strongly on L_bol than on M_BH. This is despite the fact that L_bol and M_BH vary by the same factor of 10 over the various bins we have considered in Figs. 8 and 9. These black hole and luminosity variations correspond to the same amplitude variations of λ_Edd. And yet, the PSD_amp clearly does not vary by the same amount in both cases. This result indicates that PSD_amp cannot depend only on λ_Edd, but a joint dependence with M_BH remains possible.

4.4. Modelling the joint dependence on λ_Edd and M_BH

4.4.1. PSD amplitude

In order to investigate the intrinsic dependence of the PSD amplitude on both λ_Edd and M_BH, we considered quasars in seven different BH mass bins (from 8.2 to 9.6, with ΔM_BH = 0.2). We grouped quasars in each BH mass bin in as many luminosity bins as possible (with ΔL_bol = 0.2), as long as the number of quasars in each bin was at least 50. There are in total 40 different bins. We then computed the mean PSD in each [M_BH–L_bol] bin, and we fitted them with the model defined by Eq. (2). The model fits all PSDs well and all the 40 best-fit PSD_amp parameters are plotted in Fig. 10 as a function of log(λ_Edd). The error on each point is the mean error of all the corresponding 40 best-fit parameters. As we commented above, this error should be more representative of the parameter’s uncertainty, than the error of each individual PSD parameter, as computed from the model fit to just three PSD points.

Fig. 10. PSD amplitudes for quasars in different MBH bins (points with different symbols and colours) as a function of λEdd. The effect of anti-correlation with luminosity (accretion rate) is still evident, as well as a dependence on MBH.

Points with the same symbols and colours, which fall roughly along consecutive diagonal lines, indicate the best-fit PSD_amp values when we fitted the PSD of quasars which have the same BH mass, but different L_bol (and hence λ_Edd). If PSD_amp depended only on λ_Edd, then all values should lay along the same diagonal line. This is obviously not the case. The results indicate that PSD_amp decreases both with increasing λ_Edd as well as with increasing M_BH.

We fitted the PSD_amp vs. log(λ_Edd) data in each BH mass bin (i.e. the same colour points in Fig. 10) with a line of the form:

$$\begin{array}{c}\hfill {log}_{10}[{\mathrm{PSD}}_{\mathrm{amp}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}})]=A_{\mathrm{amp}}(M_{\mathrm{BH}})+B_{\mathrm{amp}}(M_{\mathrm{BH}})log(\frac{{\lambda }_{\mathrm{Edd}}}{0.1}).\end{array}$$(3)

The best-fit lines from Eq. (3) are shown in Fig. 10. They are almost parallel, which suggests that the slope of the dependence of the power spectrum amplitude on λ_Edd is the same for all quasars, irrespective of their BH mass.

We plot the best-fit A_amp and B_amp values as a function of log(M_BH) in Fig. 11 (left and right panel, respectively). A_amp in the left panel is the PSD amplitude at log(ν) = 2.6 (days⁻¹), for quasars with log(λ_Edd) = 0.1. Clearly, this amplitude decreases with increasing BH mass. We fitted the A_amp vs. log(M_BH) data with a line of the form,

$$\begin{array}{c}\hfill A_{\mathrm{amp}}(M_{\mathrm{BH}})={\alpha }_{\mathrm{am}}+{\beta }_{\mathrm{am}}[log(M_{\mathrm{BH}})-8.9].\end{array}$$(4)

Fig. 11. Normalization, i.e. Aamp in Eq. (3) (left panel), and slope, i.e. Bamp in the same equation (right panel), plotted as a function of MBH. The parameters plotted in these panels determine how the PSD amplitude (plotted in Fig. 11) depends on BH mass and accretion rate.

This linear model provides an excellent fit to the data, as is shown by the solid line in the left (best-fit statistics reported in the same panel). The best-fit parameters are α_am = 0.45 ± 0.01, and β_am = −0.55 ± 0.03. On the other hand, B_amp does not appear to depend strongly on BH mass, in agreement with what we said above (see text below Eq. (3)). A constant line, that is

$$\begin{array}{c}\hfill B_{\mathrm{amp}}(M_{\mathrm{BH}})={\delta }_{\mathrm{am}},\end{array}$$(5)

where δ_am = −0.71 ± 0.03 is the mean of the B_amp values, provides an acceptable fit to the data (solid line in the right panel in Fig. 11; best-fit statistics reported in the same panel).

4.4.2. PSD slope

Similar trends, albeit of a smaller amplitude and higher uncertainty, can be found by repeating the same analysis with PSD_slope. Figure 12 shows the best-fit slopes of the average PSD of quasars in the 40 different bins of M_BH and L_bol defined above (colours and symbols as in Fig. 10). As before, if only λ_Edd determined PSD_slope, then all points in this plot should be located along the same line, but this does not appear to be the case. We fitted the PSD_slope vs log(λ_Edd) data in each BH mass bin (i.e. the same colour points in Fig. 12) with a line of the form:

$$\begin{array}{c}\hfill {\mathrm{PSD}}_{\mathrm{slope}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}})=A_{\mathrm{slope}}(M_{\mathrm{BH}})+B_{\mathrm{slope}}(M_{\mathrm{BH}})log(\frac{{\lambda }_{\mathrm{Edd}}}{0.1}).\end{array}$$(6)

Fig. 12. PSD slopes for quasars in different MBH bins (points with different symbols and colours) as a function of λEdd.

We therefore follow the same approach used for PSD_amp to analyse A_slope and B_slope as well. The best-fit A_slope and B_slope values as function of log(M_BH) are plotted in Fig. 13 (left and right panel, respectively). A_slope in the left panel is the PSD slope for quasars which have different BH mass and a fixed accretion rate log(λ_Edd) = 0.1. Clearly, this slope depends on BH mass: the PSD becomes steeper with increasing BH mass. We fitted the A_slope vs. log(M_BH) data with a line of the form:

$$\begin{array}{c}\hfill A_{\mathrm{slope}}(M_{\mathrm{BH}})={\alpha }_{\mathrm{sl}}+{\beta }_{\mathrm{sl}}[log(M_{\mathrm{BH}})-8.9].\end{array}$$(7)

Fig. 13. Normalization, i.e. Aslope in Eq. (6) (left panel), and slope, i.e. Bslope in the same equation (right panel), plotted as a function of MBH. The parameters plotted in these panels determine how the PSD slope (plotted in Fig. 13) depends on BH mass and accretion rate.

This line fits the data plotted in the left panel of Fig. 13 very well (fit statistics are reported on the panel). The best-fit parameters are α_sl = −1.17 ± 0.03, and β_sl = −0.5 ± 0.1. On the other hand, B_slope does not appear to depend on BH mass: the dependence of the PSD slope on λ_Edd appears to be independent of the quasar BH mass. Indeed, a constant line, that is

$$\begin{array}{c}\hfill B_{\mathrm{slope}}(M_{\mathrm{BH}})={\delta }_{\mathrm{sl}},\end{array}$$(8)

where δ_sl = −0.41  ±  0.05 is the mean of the B_slope values, provides an acceptable fit to the data plotted in the lower panel of Fig. 13 (best-fit statistics reported in the same panel).

Taking into account Eqs. (4), (5), (7), and (8), then Eqs. (3) and (6) can be re-written as follows:

$$\begin{array}{cc}& {log}_{10}[{\mathrm{PSD}}_{\mathrm{amp}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}})]={\alpha }_{\mathrm{am}}\hfill \\ \hfill & +{\beta }_{\mathrm{am}}[log\left(\frac{M_{\mathrm{BH}}}{{10}^{8.9}}\right)]+{\delta }_{\mathrm{am}}[log\left(\frac{{\lambda }_{\mathrm{Edd}}}{0.1}\right)],\hfill \end{array}$$(9)

and,

$$\begin{array}{c}\hfill {\mathrm{PSD}}_{\mathrm{slope}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}})={\alpha }_{\mathrm{sl}}+{\beta }_{\mathrm{sl}}[log(\frac{M_{\mathrm{BH}}}{{10}^{8.9}})]+{\delta }_{\mathrm{sl}}[log(\frac{{\lambda }_{\mathrm{Edd}}}{0.1})].\end{array}$$(10)

Equations (9) and (10), together with the best-fit parameter values listed above, represent our results regarding the power spectrum dependence on BH mass and λ_Edd, for all quasars in Stripe-82. However, all these results are derived using light curves at distinct filters, as long as their rest-frame wavelength lies within 3000 ± 300 Å. We discuss below the dependence of PSD on the rest frame wavelength (i.e. energy) as well.

4.5. Dependence of PSD on λ_rest

To investigate the PSD dependence on λ_rest, we repeated the same procedure as in Sect. 4.4, using light curves of various filters as long as their rest-frame wavelength would be within the λ_rest bins with |Δλ|/λ_rest = 0.2 indicated by the vertical lines in Fig. 3 (except from, λ_rest = 3000 Å, of course). So, for each waveband:

• We divided all quasars in M_BH bins, and each one of them in as many as possible L_bol bins, as long as there are at least 50 quasars in each [M_BH, L_bol] bin;

• We computed the mean PSD of all quasars in each bin, and determined PSD_amp(M_BH, L_bol), and PSD_slope(M_BH, L_bol);

• We created PSD_amp vs. log(λ_Edd) and PSD_slope vs. log(λ_Edd) plots, and we fitted the PSD_amp vs. log(λ_Edd) and the PSD_slope vs. log(λ_Edd) data (for the same BH mass) using Eqs. (3) and (6);

• We fitted A_amp, B_amp, A_slope and B_slope versus M_BH, and we calculated α_am, β_am, δ_am, and α_sl, β_sl, δ_sl. The respective values for each λ_rest are listed in Table 1.

  Table 1.

  Best-fit coefficients of Eqs. (9) and (10) for different λ_rest defined in Sect. 2.

Figure 14 shows a plot of α_am(λ_rest), β_am(λ_rest) and δ_am(λ_rest) versus λ_rest for the PSD_amp (left panels), as well as α_sl(λ_rest), β_sl(λ_rest) and δ_sl(λ_rest) versus λ_rest for PSD_slope (right panels). Solid lines in each panel represent the best linear fits, obtained by taking into account errors on both axes (where the error bar on the x axis is due to the width of the wavelength bin). In particular, β_am(λ_rest) and δ_am(λ_rest), as well as β_sl(λ_rest) and δ_sl(λ_rest) in the middle and lower panels, are acceptably fitted by their weighted means:

$$\begin{array}{c}\hfill {\beta }_{\mathrm{am}}({\lambda }_{\mathrm{rest}})=-0.52(\pm 0.01),{\delta }_{\mathrm{am}}({\lambda }_{\mathrm{rest}})=-0.71(\pm 0.02),\end{array}$$(11)

$$\begin{array}{c}\hfill {\beta }_{\mathrm{sl}}({\lambda }_{\mathrm{rest}})=-0.48(\pm 0.06),{\delta }_{\mathrm{sl}}({\lambda }_{\mathrm{rest}})=-0.40(\pm 0.04).\end{array}$$(12)

Fig. 14. PSD dependence on the rest-frame wavelength. Left panels show the PSDamp coefficients αam, βam, and δam from Eq. (9) plotted as a function of λrest. Right panels show the same plots for the PSDslope coefficients in Eq. (10). The solid lines indicate the best-fit lines to the data, the dashed line in the upper-right figure is a weighted mean (see text for details).

This result shows that the dependence of PSD_amp and PSD_slope on λ_Edd is the same at each probed wavelengths.

On the other hand, the plot in the upper-left panel shows that α_am, defining PSD_amp (see Eq. (9)), strongly depends on wavelength. In particular, it decreases with increasing rest-frame wavelength as

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{am}}({\lambda }_{\mathrm{rest}})=-0.55(\pm 0.03)-1.8(\pm 0.2)\times {10}^{-4}(\frac{{\lambda }_{\mathrm{rest}}}{3000\mathrm{\AA }}).\end{array}$$(13)

The dependence of the PSD slope, described by α_sl, is less straightforward. The weighted mean $\overline{{\alpha }_{\mathrm{sl}}}=-1.10(\pm 0.02)$ is not a good fit to the data plotted in the upper-right panel, (χ² = 18.2 for 6 d.o.f., p_null = 0.005). A linear fit returns

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{sl}}({\lambda }_{\mathrm{rest}})=-1.12(\pm 0.02)-9(\pm 3)\times {10}^{-5}(\frac{{\lambda }_{\mathrm{rest}}}{3000\mathrm{\AA }}),\end{array}$$(14)

where the slope is significant just at the 3σ level while χ² = 3.9 for 5 d.o.f., p_null = 0.56. A line improves the quality of the fit by Δχ² = 12.6 for an extra degree of freedom, which is significant according to the F-test (F-statistic = 18.3, p_null = 0.008).

5. Discussion and conclusions

We studied the continuum UV/optical variability of SDSS Stripe-82 quasars using ensemble power spectral density analysis. The final sample consists of 8042 sources with five nearly simultaneous light curves in u, g, r, i, z filters. Each object was observed in at least six yearly seasons, allowing us to sample power spectra over three timescales, in the rest-frame frequency range of log ν ∼ [ − 3.2, −2.4] day⁻¹, that is, on timescales of ∼[250, 1500] days. The available light curves have been extensively studied in the past but, to the best of our knowledge, this is the first time that they have been used to measure directly the power spectrum of quasars, and investigate how it depends on all major AGN physical parameters, at different UV/optical spectral bands. Although we estimate the power spectrum of individual objects in just three frequencies, the large number of quasars in the sample, with parameters distributed over a broad range, allowed us to analyse how variability properties, described by PSD_amp and PSD_slope, depend on quasars physical properties.

We divided the sample into bins of [M_BH, L_bol, z] containing at least 50 sources, and studied the variability properties at seven different rest-frame wavelengths λ_rest, of width |Δλ|/λ_rest = 0.2 (see Fig. 3).

1. The ensemble power spectrum in all [M_BH, L_bol, z] bins is well fitted by a simple power-law of the form:

   $$\begin{array}{c}\hfill \mathrm{PSD}(\nu )={\mathrm{PSD}}_{\mathrm{amp}}{(\nu /{\nu }_0)}^{{\mathrm{PSD}}_{\mathrm{slope}}},\end{array}$$(15)

   over the investigated range of frequencies and wavelengths.

2. The power spectrum shows no significant evidence for a dependence on redshift in the range 0.5 ≤ z ≤ 2.5 at any M_BH, L_bol, and λ_rest confirming that quasar variability properties do not evolve with time. Combining all redshifts together, we are able to study the dependencies of the PSD parameters, namely the amplitude and slope, on BH mass (8.2 ≤ log(M_BH)≤9.6) and accretion rate (−1.6 ≤ log(λ_Edd)≤ − 0.2).

3. At each rest-frame wavelength, both PSD_amp and PSD_slope show dependencies with M_BH and λ_Edd, independently (see Eqs. (9) and (10)). Such relations define a plane connecting the quasar variability properties (measured through the PSD), with BH mass and accretion rate, as shown in Fig. 15. Both the amplitude and the slope are anti-correlated with mass and accretion rate, i.e. power spectra become steeper and their amplitude decreases for larger black hole masses and higher accreting objects.

4. We found that the power spectrum also depends on the rest frame wavelength. This is shown in Fig. 16 where we plot the ensemble PSD of quasars of all BH masses and accretion rates, but at different λ_rest (PSDs for the various M_BH and λ_Edd were normalized to log(M_BH) = 8.9 and log(λ_Edd) = 0.1 using Eqs. (9) and (10)). Clearly, the PSD amplitude decreases significantly with increasing wavelength. We find some indication that PSD_slope may also depend on λ_rest, with power spectra becoming steeper at longer wavelengths. However, the significance for this effect is not as strong as for the other relation, and any PSD slope variations should be of relatively small amplitude (see upper panels of Fig. 14).

Fig. 15. Variability planes for the PSD amplitude and the PSD slope of the average power spectrum of the quasars at λrest = 3000 Å, computed using Eqs. (9) and (10) (left and right panels, respectively).

Fig. 16. Ensemble PSDs of quasars of all MBH and λEdd, with different λrest (indicated with different colours). Power spectra of quasars with various MBH and λEdd are normalized to log(MBH) = 8.9 and log(λEdd) = 0.1 using Eqs. (9) and (10). Solid lines show the best-fit linear models to the data (best-fit slopes are in agreement with the best-fit results plotted in the top panel in Fig. 14).

Our final results regarding PSD_amp and PSD_slope can be summarized by the following equations:

$$\begin{array}{cc}& {log}_{10}[{\mathrm{PSD}}_{\mathrm{amp}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}},{\lambda }_{\mathrm{rest}})]=-0.55(\pm 0.03)\hfill \\ \hfill & -1.8(\pm 0.2)\times {10}^{-4}(\frac{{\lambda }_{\mathrm{rest}}}{3000\mathrm{\AA }})-0.52(\pm 0.01)[{log}_{10}(\frac{M_{\mathrm{BH}}}{{10}^{8.9}})]\hfill \\ \hfill & -0.71(\pm 0.02)[{log}_{10}(\frac{{\lambda }_{\mathrm{Edd}}}{0.1})],\hfill \end{array}$$(16)

and

$$\begin{array}{cc}& {\mathrm{PSD}}_{\mathrm{slope}}(M_{\mathrm{BH}},{\lambda }_{\mathrm{Edd}},{\lambda }_{\mathrm{rest}})=-1.12(\pm 0.02)\hfill \\ \hfill & -9(\pm 3)\times {10}^{-5}(\frac{{\lambda }_{\mathrm{rest}}}{3000\mathrm{\AA }})-0.48(\pm 0.06)[{log}_{10}(\frac{M_{\mathrm{BH}}}{{10}^{8.9}})]\hfill \\ \hfill & -0.40(\pm 0.04)[{log}_{10}(\frac{{\lambda }_{\mathrm{Edd}}}{0.1})].\hfill \end{array}$$(17)

In using these results to constrain physical models of quasar variability, it must be stressed that our relations characterize the quasar PSD properties on timescales of 2–6 yr (observer frame), and were derived studying quasars with luminosity (45.3 ≤ log(L_bol)≤46.5), BH mass (8.2 ≤ log(M_BH)≤9.6) and accretion rate (−1.6 ≤ log(λ_Edd)≤ − 0.2).

5.1. Comparison with previous studies

We analyse the same sample of SDSS quasars as MacLeod et al. (2010). However, while we used a model-independent approach and measured power spectra directly from the data, they assumed a DRW model and found its best-fit parameters, namely the decorrelation time τ and the long-term rms variability amplitude SF_∞. In general, our results are in agreement with MacLeod et al. (2010) for what concerns the trends of the variability parameters (the PSD amplitude and slope, in our case) on AGN luminosity, BH mass, accretion rate and rest-frame wavelength. Like MacLeod et al. (2010), we also find that these parameters do not depend on redshift.

On the other hand, we find PSD slopes around −0.8/−1.2 (see left panel of Fig. 13), depending on BH mass, while the DRW model implies that the PSD slope is −2, irrespective of the BH mass, above the damping frequency. This discrepancy could imply a problem with the DRW assumption, as other works also found flatter slopes from ensemble studies (e.g. Vanden Berk et al. 2004). However, as we have stressed before, our results are based on the PSD analysis in a particular frequency range. It is possible that this frequency range is close to the range where the DRW power spectrum bends from a slope of −2 to zero, which could explain the flatter PSD slopes. In order to test this possibility, we used the scaling relations provided by (MacLeod et al. 2010, their Eq. (7) and Table 1) to derive the DRW parameters for quasars with log(λ_Edd) = 0.1 and log(M_BH) = 8.3, 8.9, and 9.5, at 3000 Å. We then determined the PSD resulting from the DRW model for these parameters, following Eq. (2) of Kozlowski et al. (2010):

$$\begin{array}{c}\hfill P(\nu )=\frac{2{\widehat{\sigma }}^2{\tau }^2}{1+{(2\pi \tau \nu )}^2},\end{array}$$(18)

where $\widehat{\sigma }={\mathrm{SF}}_{\infty }/\sqrt{\tau }$. Figure 17 shows the resulting DRW models (dashed lines) and the PSDs we have measured for quasars with the same BH mass and accretion rate (points and lines with the same colours correspond to models and measured PSDs of quasars with the same parameters).

Fig. 17. Power spectra of quasars at 3000 Å with log(λEdd) = 0.1 and three different mass bins: log(MBH) = 8.3, 8.9, and 9.5 (different colours). Points show PSD values measured in this work, while dashed lines are the best-fit DRW models from MacLeod et al. (2010).

Model predictions do not seem to be in good agreement with our power spectrum measurements as both the relative amplitudes and the slopes are different. The flatter slopes might be partly explained considering aliasing effects on the power spectra; however in order to be consistent with DRW models we need to have the DRW break exactly in the sampled frequency range (see Appendix B). In fact, even MacLeod et al. (2010) were not able to distinguish between a spectral index of −2 or −1 with the limited Stripe-82 sampling, albeit they assumed a DRW to derive scaling relations. Furthermore, the damping timescales predicted by MacLeod et al. (2010) are likely underestimated due to the finite length of light curves (see discussion in Kozłowski 2017 and Suberlak et al. 2021). If this was the case, the discrepancy between DRW and our results would increase, as model slopes would be closer to −2 and aliasing effects on our power spectra would be minimal. We therefore conclude that DRW are not fully consistent with our results, although putting definitive constraints on the intrinsic PSD shape (i.e. whether and where it is bending) with Stripe-82 data is difficult without making further assumptions.

Recently, Arévalo et al. (2024) modelled the quasar variability using data from the ZTF archive, with a range of physical properties similar to those analysed in our work. They found that the PSD is well fitted by a broken power law, where both the amplitude and the break frequency depend on BH mass and accretion rate. The best-fit model predicts power indices of −3 (at high frequencies) and −1 (at low frequencies), while the DRW model resulted in a significantly worse fit. Figure 18 shows a comparison between the PSD of the SDSS quasars from our analysis and their model predictions. Data points and dashed lines are as defined in Fig. 17. The overall absolute normalization of the Arévalo et al. (2024) models has been arbitrarily shifted to broadly match our PSD estimates, while the relative differences are constrained by the scaling relations derived in Arévalo et al. (2024). We also point out that these relations are defined at λ_rest = 2900 Å, while our PSD estimates are at λ_rest = 3000 Å.

Fig. 18. Power spectra of quasars at 3000 Å with log(λEdd) = 0.1 and three different mass bins: log(MBH) = 8.3, 8.9, and 9.5 (different colours). Points show PSD values measured in this work, while dashed lines are the best-fit models from Arévalo et al. (2024), at 2900 Å.

The Arévalo et al. (2024) model predictions also do not seem to reproduce well our results, albeit the discrepancy seems to be smaller than with the DRW models from MacLeod et al. (2010). Our flatter slopes seem to suggest a higher frequency break although, again, the limited frequency range explored in our analysis does not allow to robustly constraint the position of the break. A more thorough analysis on both longer and shorter timescales will be subject to a future work.

5.2. Addressing the meaning of the scaling relations

The scaling relations presented in this work describe the dependence between the main physical parameters of quasars and their variability properties, as determined by a model independent power spectrum analysis. As we argued above, such relations are important, as they should allow us to test different theoretical models. Testing specific models is beyond the scope of the present work. However, we can use the scaling relations and try to understand, in broad terms, possible physical reasons behind them.

5.2.1. The dependence of PSD amplitude and slope on BH mass

The anti-correlation of PSD_amp and PSD_slope with BH mass could be due to the fact that the power spectrum is the same in all AGN, but the intrinsic timescales are proportional to M_BH. This may be natural, as a larger BH mass implies a larger size for the accretion disc itself (the gravitational radius scales proportionally to M_BH). After all, the accretion disc characteristic timescales are proportional to M_BH. In this case, it is possible that the relations we found may be due to the fact that the observed power spectra sample different intrinsic physical timescales in each source. Similar considerations have also been proposed in other recent works using different data, such as Arévalo et al. (2024) and Tang et al. (2023).

In order to test this idea, we transformed the observed light curve of each object to an intrinsic light curve by dividing the observed times by the light travel time of the gravitational radius of the central BH, that is

$$\begin{array}{c}\hfill T_g=\frac{G{M}_{\mathrm{BH}}}{c^3},\end{array}$$(19)

where G is the universal constant of gravitation, and c is the speed of light. As our power spectra are computed using the periodogram defined in Eq. (1), this normalization implies that both the sampled frequencies and the power spectrum amplitude are respectively multiplied and divided by T_g. We applied this normalization to seven groups of quasars with 8.2 < log(M_BH) < 9.5 and log(λ_Edd) = −1.0, using light curves at 3000 Å rest-frame. The result is shown in Fig. 19.

Fig. 19. Ensemble PSDs of quasars with λrest = 3000 Å, log(λEdd) = − 1.0 and various MBH (points with different colours on the plot). Black points indicate the universal PSD of all the quasars, computed when light curves are normalized with the gravitational timescale, Tg. Black dashed and solid lines show the best-fit single and bending power law models to the black points, respectively. Coloured lines show the same best-fit models re-scaled in frequency by 1/Tg, for each of the BH mass listed in the bottom left part of the plot (see Sect. 5.2 for details).

Points with different colours in this figure show the original (i.e. non-scaled) power spectra of the different BH mass groups (note that the y-axis in this figure shows log₁₀[PSD(ν)×ν]). Black points show the scaled power spectra, as we explained above. The alignment of all the PSD estimates after the normalization with T_g is quite good. The fact that all the points become almost aligned after the PSD normalization with T_g indicates that the assumption of a power spectrum which has the same shape in all quasars, irrespective of their BH mass, while timescales depend linearly on M_BH is reasonable.

The dashed black line in Fig. 19 shows the best linear fit to the black points. We also fitted the power spectrum with a bending power law (e.g. McHardy et al. 2004):

$$\begin{array}{c}\hfill P(\nu )=N{\nu }^{-1}[1+(\frac{\nu }{{\nu }_b}{)}^{\alpha -1}{]}^{-1},\end{array}$$(20)

as shown by the solid black line. This is a model that fits well the X-ray PSDs of nearby Seyferts. It corresponds to a power law of slope −1 at low frequencies, which steepens to a slope of α at frequencies higher than ν_b. The best-fit slope of the linear fit is −1.34(±0.03), while the best-fit high frequency slope of the bending power law model is −1.6(±0.1), and log₁₀(ν_b) = − 2.90(±0.01). The high frequency slope is rather flat, and it is for this reason that both the linear and the bending power-law models fit the data equally well. Statistically speaking, neither of the two fits are acceptable (χ² = 46.5/19 d.o.f. and χ² = 40.6/18 d.o.f. for the linear and the bending power-law fits, respectively). However, as the figure shows, both models fit the normalized PSD rather well, with no systematic trends in the best-fit residuals.

The solid and dashed coloured lines in Fig. 19 show the best-fit power law and bending power fits to the normalized PSD (i.e. the black solid and dashed lines) scaled back to the PSD of the various mass groups plotted with the same colours. To re-scale the PSD models we have kept the best-fit parameters the same, and we re-scale the frequencies by 1/T_g, that is

$$\begin{array}{c}\hfill {\nu }_{M_{\mathrm{BH}}}={\nu }_{\mathrm{norm}}\frac{1}{G{M}_{\mathrm{BH}}/c^3},\end{array}$$(21)

where ν_norm are the frequencies of the universal quasar PSD (shown with the black points in Fig. 19), and ν_MBH the frequencies of the PSD of quasar of mass M_BH. The re-scaled PSDs agree well with the observed PSDs, as expected. Again, this supports the existence of a universal PSD shape for all quasars, where frequencies scale with BH mass, while normalization and slope(s) are fixed (at any given wavelength and accretion rate).

5.2.2. The dependence of PSD on λ_Edd, and λ_rest

Even if the hypothesis of a universal PSD which scales with BH mass in frequency could explain the dependence of the power spectrum on M_BH, it is difficult to understand the dependence on accretion rate. As the accretion rate increases, so does the surface area of the disc that emits at a given wavelength. If the UV/optical variability is due to variations that take place in the inner disc only (i.e. at distances within a given radius same in all objects in R_g units), then it is natural to expect that the variability amplitude should decrease as the disc emitting area increases (i.e. for larger accretion rates). At the same time, the power spectrum amplitude is expected to decrease with increasing λ_Edd as observed, in the case of X-ray reverberating discs. This is because the amplitude of the disc transfer function decreases with increasing accretion rate, at all wavelengths (see e.g. Fig. 10 in Panagiotou et al. 2022).

The dependence on λ_rest may be simpler to understand. In the case of X-ray reverberating discs, the amplitude of the disc transfer function decreases with increasing wavelengths, which implies that the PSD amplitude should also decrease (see e.g. all the figures in Appendix B of Panagiotou et al. 2022). In general, the disc-emitting area increases with increasing wavelength. Again, as we discussed above, it may be natural to expect that the variability amplitude would decrease as the disc emitting area increases (i.e. for longer wavelengths).

The discussion above does not propose a specific theoretical model for the observed UV/optical variability of quasars. It is meant to highlight a few of our results and demonstrate some potential, general constraints placed on various theoretical ideas that have been proposed in the past. We believe that, as long as the observed variations are associated with physical mechanisms that operate in the innermost region of quasars, then our results, i.e. the dependence of the PSD on BH mass, accretion rate, and rest frame energy, can constrain, strongly, any physical model for the UV/optical variability of quasars.

Acknowledgments

We thank the anonymous referee for the useful comments and suggestions. MP, VP and DD acknowledge the financial contribution from PRIN-MIUR 2022 and from the Timedomes grant within the “INAF 2023 Finanziamento della Ricerca Fondamentale”. VP also acknowledges Erasmus+ learning mobility founds 2022/2023 and the Insitute of Astronomy at FORTH, Greece, for hospitality. DD also acknowledges PON R&I 2021, CUP E65F21002880003. FEB acknowledges support from ANID-Chile BASAL CATA FB210003, FONDECYT Regular 1200495 and 1241005, and Millennium Science Initiative Program – ICN12_009. This work has been partially supported by ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU. The research has made use of the following Python software packages: Matplotlib (Hunter 2007), Pandas (van der Walt & Millman 2010), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020).

References

Appendix A: Measurement of Poisson noise in the PSD

We considered a sample of ∼5 × 10⁵ non-variable stars in Stripe-82 (Ivezić et al. 2007) to model the experimental Poisson noise in the quasars PSD. For each one of the five bands ugriz, we divided the sample of stars in bins of magnitude of width 0.2 and measured their PSD with the same prescription that we used to measure the quasars PSDs (i.e. we binned the light curves, computed the PSD at three frequencies for each star and then we constructed the ensemble PSD for the selected magnitude bin - see Sect. 3.2 for details).

Since we selected non-variable stars, we expect the ensemble PSD to be flat in all magnitude bins and spectral bands. We calculated the mean of the ensemble PSDs in each bin and we assume that this mean value is an accurate estimate of the Poisson noise contribution in each spectral band and magnitude range. Figure A.1 shows the mean ensemble PSD in each magnitude bin plotted as a function of magnitude for each filter. The solid lines show the best-fit of a third-order polynomial function to the data. As the Poisson noise should depend on the object’s flux only, we used the best-fit lines to predict the Poisson noise level of each quasar in our sample, according to its magnitude in each filter. Such a noise value is then subtracted from the PSD estimates obtained with Eq. 1, so that:

$$\begin{array}{c}\hfill PS{D}_{\mathit{corrected}}=PS{D}_{\mathit{original}}-noise.\end{array}$$(A.1)

Fig. A.1. White noise contribution as a function of magnitude for non-variable stars in the five SDSS filters. Continuous lines show third-order polynomial best-fits to the data.

Appendix B: PSD of simulated light curves

As we mentioned in Sect. 3.2, the periodogram is the usual estimator of the intrinsic PSD of a variable process. Its statistical properties are known very well (e.g. Priestley 1981) and they have been studied extensively even in the case of red-noise like processes which are quite common in astronomy (e.g. Papadakis & Lawrence 1993 and references therein). We use light curves with six points to compute the periodogram in three frequencies, for each quasar in the sample. Nevertheless, irrespective of the small number of points in the light curves, the periodogram statistical properties are only marginally affected.

To demonstrate this, we simulated stochastic light curves with a cadence of 1 day and a length of 20 years, assuming a DRW model with damping timescales in the range [10² − 10⁵] days (i.e. we went from a broken power-law, to a single power law with a slope of –2). The procedure for simulating DRW light curves is described in MacLeod et al. (2010) and Kovačević et al. (2021). For each fixed DRW τ and σ, we created 50 simulated light curves and we applied the same sampling pattern as that of the observed SDSS Stripe-82 light curves. We then processed the simulated data with the same procedure applied to real data (i.e. we binned the time series in yearly seasons) and we computed the periodograms at three frequencies with the resulting light curves which consist of six points.

As an example, we show the resulting PSDs of simulations with τ = 10³ days in Fig. B.1. Small grey dots represent the periodograms of the 50 original light curves, while the red points indicate the mean periodogram at each frequency. It is clearly evident that the mean of the periodogram estimates is equal to the intrinsic power-law model at all frequencies (solid black line, Eq. 18). This result is expected, as the periodogram is an unbiased estimator of the power spectral density function. Blue dots (and squares) on the same figure indicate the 50 periodograms (and their mean) which are calculated using the binned light curves (which consist of only six points). Again, they follow the intrinsic model for the frequency range they are sampling, showing that our approach is robust.

Fig. B.1. PSD of 50 simulated DRW light curves with a cadence of 1 day and a length of 20 years (grey points). Red points show the mean PSD values, at each frequency, on the original data. Blue squares are the ensemble PSD estimates from light curves binned as in this work. The black solid line shows the intrinsic PSD shape.

The slight flattening in the continuous PSD at very high-frequencies (i.e. at the right-side of the figure) is likely due to aliasing, which is due to the fact that higher (non-sampled) frequencies fold back to lower frequencies. This might also be affecting the highest frequency point on the binned PSD, but for power-law-like power spectra with a slope steeper than −1 at high frequencies the effect is typically small. As power decreases logarithmically with frequency, the amount of power folded back to lower frequencies is small with respect to the observed value. The steeper the power-law slope, the lesser the effect, with the exact contribution depending on the intrinsic PSD shape and the sampled frequency range. Aliasing effects should also be suppressed as each one of the six points in the light curves is the result of binning the original light curves over discrete time intervals of a couple of months. This procedure significantly reduce the amount of power which is aliased back to the sample frequencies, as shown by van der Klis (1988).

Our simulations show that the worst case scenario should appear if we sample an intrinsic PSD with a slope flatter than −1 at high frequencies and a constant power at low frequencies (e.g. if the sampled frequencies are close the DRW break frequency). This could lead to underestimation in the recovered PSD slope up to ∼0.5. However, this systematic can be easily reduced when more frequencies are available (i.e. with better sampled light curves), or corrected once a particular PSD model is assumed in order to explain the observed PSD slopes.

All Tables

All Figures

	Fig. 1. Distribution of observed quasar magnitudes in the five SDSS filters. Dashed lines enclose the magnitude range of sources used in this analysis.
In the text

	Fig. 2. Distribution of black hole mass, bolometric luminosity, and redshift for the quasars in our sample.
In the text

	Fig. 3. Distribution of rest-frame wavelengths (λrest) of all light curves, for all quasars in the sample. The vertical dashed lines indicate the centres of the λrest bins at which we computed the PSD of quasars.
In the text

	Fig. 4. Light curves of two quasars in g-band with a similar BH mass (8.8 < log(MBH) < 9.0), luminosity (46.1 < log(Lbol) < 46.3), and redshift (1.1 < z < 1.3). Red points indicate the binned light curve.
In the text

	Fig. 5. Noise-corrected PSD of the two light curves shown in Fig. 4. The solid lines indicate the best-fit line to the data.
In the text

	Fig. 6. Logarithm of the noise-corrected periodogram of all quasars with [8.8 < log(MBH) < 9.0], [46.1 < log(Lbol) < 46.3], and [1.2 < z < 1.4]. Small blue dots indicate single PSD estimates, while black points show the logarithm of the mean PSD (with its error), and the red line is the best-fit linear model to the ensemble PSD.
In the text

	Fig. 7. Ensemble PSD of quasars with 8.8 < log(MBH) < 9.1 and 46.0 < log(Lbol) < 46.3 at λrest = 3000 Å in ten redshift bins (different symbols and colours). The average redshift for each bin is reported in the legend.
In the text

	Fig. 8. Ensemble PSDs of quasars with 8.8 < log(MBH) < 9.0 and all redshifts between 0.5 and 2.5 in different log(Lbol) bins. The legend reports the average luminosity in each bin of width 0.2.
In the text

	Fig. 9. Average PSDs of quasars with 46.1 < log(Lbol) < 46.3 and all redshifts between 0.5 and 2.5 in different log(MBH) bins. The legend reports the average mass in each bin of width 0.2.
In the text

	Fig. 10. PSD amplitudes for quasars in different MBH bins (points with different symbols and colours) as a function of λEdd. The effect of anti-correlation with luminosity (accretion rate) is still evident, as well as a dependence on MBH.
In the text

	Fig. 11. Normalization, i.e. Aamp in Eq. (3) (left panel), and slope, i.e. Bamp in the same equation (right panel), plotted as a function of MBH. The parameters plotted in these panels determine how the PSD amplitude (plotted in Fig. 11) depends on BH mass and accretion rate.
In the text

	Fig. 12. PSD slopes for quasars in different MBH bins (points with different symbols and colours) as a function of λEdd.
In the text

	Fig. 13. Normalization, i.e. Aslope in Eq. (6) (left panel), and slope, i.e. Bslope in the same equation (right panel), plotted as a function of MBH. The parameters plotted in these panels determine how the PSD slope (plotted in Fig. 13) depends on BH mass and accretion rate.
In the text

	Fig. 14. PSD dependence on the rest-frame wavelength. Left panels show the PSDamp coefficients αam, βam, and δam from Eq. (9) plotted as a function of λrest. Right panels show the same plots for the PSDslope coefficients in Eq. (10). The solid lines indicate the best-fit lines to the data, the dashed line in the upper-right figure is a weighted mean (see text for details).
In the text

	Fig. 15. Variability planes for the PSD amplitude and the PSD slope of the average power spectrum of the quasars at λrest = 3000 Å, computed using Eqs. (9) and (10) (left and right panels, respectively).
In the text

	Fig. 16. Ensemble PSDs of quasars of all MBH and λEdd, with different λrest (indicated with different colours). Power spectra of quasars with various MBH and λEdd are normalized to log(MBH) = 8.9 and log(λEdd) = 0.1 using Eqs. (9) and (10). Solid lines show the best-fit linear models to the data (best-fit slopes are in agreement with the best-fit results plotted in the top panel in Fig. 14).
In the text

	Fig. 17. Power spectra of quasars at 3000 Å with log(λEdd) = 0.1 and three different mass bins: log(MBH) = 8.3, 8.9, and 9.5 (different colours). Points show PSD values measured in this work, while dashed lines are the best-fit DRW models from MacLeod et al. (2010).
In the text

	Fig. 18. Power spectra of quasars at 3000 Å with log(λEdd) = 0.1 and three different mass bins: log(MBH) = 8.3, 8.9, and 9.5 (different colours). Points show PSD values measured in this work, while dashed lines are the best-fit models from Arévalo et al. (2024), at 2900 Å.
In the text

Fig. 19. Ensemble PSDs of quasars with λrest = 3000 Å, log(λEdd) = − 1.0 and various MBH (points with different colours on the plot). Black points indicate the universal PSD of all the quasars, computed when light curves are normalized with the gravitational timescale, Tg. Black dashed and solid lines show the best-fit single and bending power law models to the black points, respectively. Coloured lines show the same best-fit models re-scaled in frequency by 1/Tg, for each of the BH mass listed in the bottom left part of the plot (see Sect. 5.2 for details).

In the text

	Fig. A.1. White noise contribution as a function of magnitude for non-variable stars in the five SDSS filters. Continuous lines show third-order polynomial best-fits to the data.
In the text

	Fig. B.1. PSD of 50 simulated DRW light curves with a cadence of 1 day and a length of 20 years (grey points). Red points show the mean PSD values, at each frequency, on the original data. Blue squares are the ensemble PSD estimates from light curves binned as in this work. The black solid line shows the intrinsic PSD shape.
In the text