© The Authors 2026

1. Introduction

Active galactic nuclei (AGNs) are among the most energetic sources in the Universe, powered by supermassive black holes (BHs) actively accreting matter at the centre of galaxies (Netzer 2015). One of their striking features is variability, observed across the entire electromagnetic spectrum, both in the continuum and the emission lines, with timescales ranging from minutes to years (e.g. Ulrich et al. 1997; Giveon et al. 1999; Cackett et al. 2021). The study of AGN variability is becoming increasingly relevant, especially in the era of large time-domain surveys. Indeed, variability allows us to understand the geometry and physics of the accretion process and has proven to be a valuable selection tool (e.g. Butler & Bloom 2011; De Cicco et al. 2022, 2025; Savić et al. 2023). Most of the studies focus on continuum X-ray and ultraviolet/optical (UV/optical) variability measured from light curves, which show a stochastic behaviour and directly probe the accretion disc and the X-ray corona around the central source. Although there is debate about the interplay between these two AGN regions, the correlations between variability and AGN properties (e.g. BH mass, accretion rate, wavelength) and the origin of variability itself, there is general consensus that power spectral density functions (PSDs or power spectra) or equivalent measurements (e.g. structure function, SF) display a red-noise trend and are usually modelled with one or more power laws (see Paolillo & Papadakis 2025 and references therein for a recent review of continuum variability). Optical light curves are usually modelled with a damped random walk (DRW; Kelly et al. 2009). This is a first-order continuous-time autoregressive moving average (CARMA) process, which predicts a PSD with a spectral index of –2, flattening to 0 after a damping timescale, τ. Although the model seems to be consistent with light curves of the order of a few years, deviations have been observed on both shorter and longer timescales (e.g. Mushotzky et al. 2011; Guo et al. 2017), as well as possible biases in the proper recovery of DRW parameters due to sampling effects (e.g. Emmanoulopoulos et al. 2010; Kozłowski 2017). To improve our understanding of AGN continuum variability, we need large samples of light curves that probe multiple timescales and to explore new analysis techniques, such as fitting the data with higher-order CARMA models (e.g. Yu et al. 2022b), unsupervised machine learning analysis (e.g. Tachibana et al. 2020), or directly measuring the PSD with no a priori model assumptions.

In Petrecca et al. (2024), hereafter P24, we studied a sample of 8042 spectroscopically confirmed quasars (i.e. the most luminous among unobscured AGNs) from the SDSS Stripe-82 region, first selected and analysed by MacLeod et al. (2010, 2012). Unlike many of the previous works on optical light curves, P24 directly computed the power spectra via traditional Fourier techniques and studied relations between the ensemble PSD and quasar physical properties. They found that variability does not depend on redshift, while both the PSD amplitude and slope depend on the accretion rate, BH mass, and rest-frame wavelength (e.g. Eqs. 16 and 17 in P24). They also discussed the possibility of a universal PSD shape for all quasars, where frequencies scale with the BH mass, while normalisation and slopes are fixed (at any given wavelength and accretion rate). While similar ideas have also been proposed by other studies (e.g. Tang et al. 2023; Arévalo et al. 2024), limitations in cadence, temporal baseline, and survey coverage have hindered the ability to unambiguously determine the intrinsic PSD shape of all AGN and its dependencies.

The forthcoming Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) expected to start by the end of 2025 at the Vera C. Rubin Observatory will mitigate many of these problems. It will observe the entire southern sky every three-four nights for ten years across ugrizy filters and a 5σ point source depth of ∼24.7 mag_AB in the r band for a single visit (Bianco et al. 2022). In addition, five deep drilling fields (DDF¹) covering ∼60 deg² will be observed even more intensively in terms of depth and cadence. The LSST is expected to detect tens of millions of AGNs, allowing for a detailed analysis of light curves on multiple timescales, including fast and extreme variability (e.g. Raiteri et al. 2022; Komossa et al. 2024), and precise time-delay measurements (Kovačević et al. 2022; Czerny et al. 2023). However, challenges such as seasonal gaps, uneven sky coverage, and missing data due to other survey-related factors will inevitably persist. There is a general community effort based on simulations and archival data aimed at supplementing our understanding of what to expect from the LSST data, while preparing a machinery of software, techniques, and ancillary datasets ready to work with the first Rubin data releases. Our PSD framework, introduced in P24, is among these methods, but other aspects such as optimising light-curve pre-processing to extract more variability features also need to be inspected.

This work is a follow-up to the P24 paper. We present an ensemble PSD analysis of quasars combining data from multiple surveys. The starting sample consists of the SDSS Stripe-82 quasars from P24, which all feature information related to redshift, BH mass (M_BH), luminosity (L_bol), and accretion rate (λ_Edd = L_Bol/L_Edd). The combination of light curves from different facilities allows us to sample a larger frequency range with both shorter and longer timescales, thus having a better constraint on the PSD shape. Moreover, this methodology will be most likely applied in the future to LSST data to extend the observing baseline, but taking advantage of the better cadence and wavelength coverage similar to the SDSS data studied in P24. Finally, we also present a library of calibrated light curves covering ∼20 years, expanding on previous efforts by Rumbaugh et al. (2018), Suberlak et al. (2021), and Stone et al. (2022).

This paper is organised as follows. In Sect. 2, we present the data used for our work. Section 3 describes the PSD estimates from multiple surveys, while in Sect. 4 we discuss the broad-band PSD shape. We investigate how the power spectra depend on BH mass and accretion rate in Sect. 5 and we report our final considerations in Sect. 6.

2. Data

In an effort to build a legacy dataset of archival light curves in preparation for LSST, we started collecting all publicly available data complementing the SDSS Stripe-82 region already studied for the LSST AGN Science Collaboration Data Challenge (Yu et al. 2022a; Savić et al. 2023) by P24 and others (e.g. MacLeod et al. 2010; Kozłowski 2016). As our sample of quasars from Stripe-82 consists of bright point-like sources with mag < 21.5 (see Fig. 1 in P24), it is also possible to recover most of them with shallower and lower resolution surveys. The point-like nature of quasars enables an easier comparison of PSF photometry from different facilities, with minor issues related to the host-galaxy contribution.

The most recent, ongoing time-domain survey is the ZTF (Bellm et al. 2019; Graham et al. 2019), which shares many similarities with the future LSST, including alert production and difference image analysis (DIA), although at a lower depth. It started in March 2018, using the Samuel Oschin 48-inch Telescope at the Palomar Observatory to observe the northern equatorial sky (including Stripe-82) in g,r bands with a few nights cadence. Additional i band observations are also available (although with a more erratic cadence), as well as repeated g,r exposures in the same night. Currently, public data releases are updated every two months and include raw CCD-based images and calibrations in FITS format, reference and difference images from DIA, updated objects database, and single-epoch photometry (Masci et al. 2019). We used ZTF DR19 to cross-match the object table with the catalogue of SDSS Stripe-82 quasars. The match returned 8250 sources within 0.5 arcsec (mean separation 0.12 arcsec, standard deviation 0.06 arcsec), including all the 8042 quasars used by P24. We downloaded PSF photometry in g,r bands up to July 2023 (i.e. the last visits in DR19) for each object with an ad hoc python code querying the IRSA database².

Although the ZTF high-cadence data are ideal for studying variability on short timescales, the temporal baseline of the DR19 used in this work is ∼5.5 years and it does not add any information on timescales longer than those probed with the SDSS light curves (the case is the same when using the most recent data release at the time of writing, which is about seven years, but this could change in the future if the survey progresses and overlaps with the LSST). Furthermore, the close to nine year gap between the two surveys prevents us from simply combining the light curves to constrain even longer times. Filling the gap requires gathering archival data from other facilities, as proposed by Suberlak et al. (2021); see also Figure 2 of Paolillo & Papadakis (2025) for an updated list. However, the differences between the surveys in terms of sky coverage, photometric systems, sensitivity, and resolution make combining and homogenising the data a non-trivial task; thus, the selection of a good sample is critical to minimise any bias in the analysis. Among recent time-domain surveys that overlap with the Stripe-82 region, we have Catalina Real Time Survey (CRTS; Drake et al. 2009), Palomar Transient Factory (PTF/iPTF survey; Law et al. 2009), Panoramic Survey Telescope And Rapid Response System Survey (Pan-STARRS-1, or PS1 for short; Kaiser et al. 2010; Chambers et al. 2016), All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014), NASA Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018), and the Gaia space telescope (Gaia Collaboration 2016).

In an attempt to use all the available information, we explored the light curves of each survey and reached the same conclusion as Suberlak et al. (2021), namely, to exclude CRTS and PTF from the analysis. In fact, their shallower photometry and large photometric uncertainties affect the quality of the observed PSDs, whereas the temporal coverage is similar to the deeper PS1 observations (actually CRTS has a better sampling, but this is spoiled by the lower photometric quality and the need to reduce the sample to fewer good sources). Similar considerations also apply to ATLAS and ASAS-SN, whose limiting magnitude is shallower than the typical ∼20 mag_AB, r of Stripe-82 quasars (see Fig. 1 in P24). On the other hand, we included ZTF data in our analysis, as they provide a fundamental contribution to the high-frequency PSD, as shown in Sect. 3.1, and cover a different time period. As a further ground-based survey, we included PS1 as well, a wide-field imaging facility that contains measurements in five filters, grizy, for 30 000 square degrees of the sky north of declination –30° (including Stripe-82). The data set is composed of roughly ∼12 epochs for each filter, from May 2010 to March 2014, publicly available through the DR2³ (Flewelling et al. 2020). A spatial cross-match with the Stripe-82 quasar catalogue returns all the sources in common with ZTF, but there is still a significant gap of ∼4–5 years between the observations.

We also used light curves from the ESA’s Gaia satellite, designed for an accurate astrometric, photometric, and spectroscopic map of the Milky Way, therefore exhibiting extremely good resolution and photometric quality, albeit less deep and employing broad-band filters. Having started its operations in late 2014 and just concluding its nominal mission in early 2025, Gaia is the perfect survey to complement the previous observations. We cross-matched the Stripe-82 quasars with the Gaia DR3 variable AGN catalogue (Carnerero et al. 2023), covering 34 months from July 2014 to May 2017. We found a match for 4037 objects. This low number is due to many factors, including the lower magnitude limit (mag_G < 21) and the fact that Gaia only published light curves for sources detected as variable within the survey. The last limitation is the most severe, considering that AGN variability tends to be larger on longer timescales and a few years might not be enough for a complete selection (e.g. see discussion in De Cicco et al. 2019). However, by looking at the quasar physical properties from SDSS+PS1+ZTF and Gaia (Fig. 1), we see that they cover the same range of L_bol and M_BH with a similar distribution; thus, the analysis on the limited sample is not expected to introduce any severe bias. Future data releases could allow for the recovery of the light curves for the missing sources and the addition of more epochs, possibly overlapping with ZTF and improving the analysis.

Fig. 1. Distributions of bolometric luminosity (left panel) and BH mass (right panel) for the 8042 SDSS+PS1+ZTF quasars (red) and the subset of 4037 sources in common with Gaia DR3 (blue).

To summarise, we cross-matched data from SDSS (1998–2008), PS1 (2010–2014), Gaia (2014–2017), and ZTF (2018–2023) to study AGN variability and build a legacy dataset in preparation for the LSST. The ZTF’s high cadence allows for a detailed analysis of high-frequency PSDs, while the full combined light curves cover more than 20 years, probing lower frequencies. Our main objective is to compute the power spectrum of the quasars in the widest possible frequency range. To do so, the AGN flux from the various catalogues must be cross-calibrated. This is necessary for various reasons. First, the filters used by the SDSS, PS1 and ZTF surveys are not identical (in width and shape). Second, Gaia measurements are made using a single broad filter, which partially overlaps but is not similar to any of the filters used by the other surveys. Therefore, it is necessary to cross-calibrate all the data we use to a single waveband (or filter) before using them to calculate the PSD.

In Appendix A, we explain in detail the process we followed to calibrate all light curves from the various surveys to the SDSS r-filter flux, which is the band minimising the corrections. We show an example of such a light curve in Fig. 2, highlighting the contributions from the different surveys. The cross-calibrated light curves prepared for this work will be made available to the community.

Fig. 2. Example of r-band light curve for a Stripe-82 quasar combining SDSS (blue), PS1 (green), Gaia (orange), and ZTF (red) data. Colour corrections and calibrations are given in Appendix A.

3. Estimating the power spectrum.

As the variability of AGNs is a stochastic process, we would ideally like to calculate the power spectrum on the longest and best sampled light curves as possible for a significant number of sources. In fact, sampling a given time range multiple times for sources with similar properties is the only way to determine the typical variability behaviour and it is the main motivation behind the effort of building a dataset as the one described in Sect. 2. To study the variability, we used traditional Fourier techniques to estimate the PSD; namely, we used the periodogram as an estimator of the power spectrum, as done in P24. Fourier analysis is fast and reliable, and the periodogram is a well-studied robust estimator of the power spectrum of stationary random processes (e.g. Press 1978; Priestley 1981; van der Klis 1988; and references therein). We can also bin a large number of periodograms (at least 50; see e.g. Papadakis & Lawrence 1993) to calculate appropriate estimators to model fit the power spectra using ordinary χ² minimisation. However, traditional Fourier analysis requires evenly sampled light curves.

We tried combining the SDSS+PS1+Gaia+ZTF data using bins of various sizes, but were unable to construct homogeneous, long, and evenly sampled light curves (with a reasonable number of missing points), spanning the whole interval from the start of the SDSS observations to the end of the ZTF data points. This was mainly due to the large gap between SDSS and the Gaia data points and the erratic sampling of quasars by PS1. For that reason, we adopted a different strategy. We used the ZTF data for the determination of the high-frequency power spectrum of the sources, and we kept the SDSS PSD estimates of P24 at low frequencies. To improve the determination of the low-frequency PSD, we calculated additional low-frequency PSD estimates using the individual Gaia and the ZTF light curves appropriately binned, using large bin sizes. We explain in detail the calculation of the high- and low-frequency PSD in the following sub-sections.

3.1. High-frequency PSD from ZTF

ZTF light curves from DR19 are composed of five yearly seasons, each with a cadence of a few nights and about five-six months in length. To measure the PSD on very short timescales, we considered each season separately assuming that they are multiple realisations of the same stationary random process. First, we converted magnitudes into fluxes in units of nJy and obtained an evenly sampled time series for each seasonal segment of the light curve, in the same way as described in P24 (i.e. starting from magnitudes in the AB photometric system with a zero point of 3631 Jy and following a standard procedure for error propagation). We adopted a bin size of 12 days, which allows us to have at least three independent visits in most of the bins. For each bin, we calculated the mean flux and the standard error on the mean (i.e. $\sigma /\sqrt{n}$, where σ is the standard deviation and n is the number of points). As is often the case with high-cadence surveys, there are missing observations within each seasonal segment (e.g. due to bad weather, full moon, etc.). We used linear interpolation to fill the gaps. We rejected segments having more than three consecutive gaps (i.e. > 30 days of missing data) and when the interpolated values accounted for more than 15% of the total number of points in the segment. This choice maximises the amount of data used for the analysis (more than ∼80% of the seasonal segments for the quasars we considered), while minimising biases introduced from the interpolation procedure, as we show on simulated light curves in Appendix B. We also tested different rejection thresholds and higher-order interpolation methods, but linear interpolation appears to be the best option for the sort term light curves. An example of the ZTF light curve with the binned data plotted on top of the observations is shown in Fig. 3.

Fig. 3. Example light curve for a ZTF quasar in the r band. Grey points are the original flux measurements, while yellow stars indicate the binned-interpolated values turning each season into an evenly sampled time series to compute the high-frequency PSD in Sect. 3.1. The second half of the first season shows data points rejected from the analysis due to the quality cuts applied. Blue squares are yearly-binned points used for the low-frequency PSD in Sect. 3.2.

3.2. Low-frequency PSD from SDSS, ZTF, and Gaia

In P24, we calculated the power spectrum at three low frequencies from the SDSS data, using six-year-long light curves with a yearly bin size (i.e. we used the last six SDSS observing seasons from the typical light curve in Fig. 2 as the epochs before were sparser; see P24 for more details). The Gaia light curves are also sufficiently well sampled to be binned with a 180-day window. The resulting light curves do not have (almost any) missing points and can be used to calculate the periodogram at three new independent frequencies, located between the SDSS and ZTF frequencies. Two additional PSD estimates at frequencies within the range sampled by SDSS could also be obtained by using the full ZTF light curves when they are binned with a bin size of one year. Contrary to Gaia and ZTF, the PS1 light curves’ sampling is more erratic and makes it difficult to to construct evenly sampled light curves without adding a large number of interpolated points. For that reason, we did not use the PS1 light curves to calculate the power spectrum.

We calculated separate periodograms for each seasonal segment in ZTF light curves, for the full GAIA light curves, and for the full ZTF light curves, as in P24, namely,

$$\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)

where N is the number of points in the light curve, x(t_i) are the flux values and $\overline{x}$ is the mean flux. Defined in this way, the PSD estimates have units of 1 per unit frequency (i.e. days). The power spectrum is computed in the rest frame of each source, that is, all t_i’s in the above equation are divided by (1 + z), including the time window, Δt. Periodograms are defined in the usual set of frequencies, f_j = j/(NΔt), with j = 1, 2, …, N/2.

We could combine the periodograms from the seasonal segments of the ZTF light curves (there are typically five of them for each object), from the SDSS, the Gaia and the long-term binned ZTF light curves to construct a broad frequency band power spectrum for each quasar in our sample. However, periodogram estimates are distributed as χ² variables with 2 degrees of freedom, and their error depends on the intrinsic power spectrum, which is unknown (e.g. Papadakis & Lawrence 1993). Therefore, they are not suitable for fitting models to the PSD using χ² techniques. For that reason, we followed the same approach as with P24, and we calculated ensemble power spectra for many quasars which have similar BH mass, luminosity, and redshift (as a proxy of rest-frame wavelength) as we explain below.

4. The broad-band power spectrum of quasars

The ensemble PSD analysis requires averaging the periodograms of quasars with similar physical properties, as outlined in P24. The availability of information on black hole mass, luminosity, redshift, and accretion rate (from the SDSS DR7 quasar catalogue by Shen et al. (2011), see Sect. 2 of P24 for more details) allows us to explore how the PSD features (slope, normalisation, and bending/break frequencies) depend on these parameters. Compared to P24, the joint contribution to the PSD from SDSS, Gaia and long-term ZTF at low frequencies, and from the seasonal ZTF segments at high frequencies enable us to probe the quasar PSD over a wider frequency range, for any given bin of quasars with similar physical properties. P24 showed that quasar variability does not depend on redshift (at any given wavelength), but depends on M_BH, L_bol (or λ_Edd = L_Bol/L_Edd), and rest-frame wavelength. However, the availability of simultaneous SDSS observations in ugriz bands easily allowed P24 to select many bins of [M_BH, λ_Edd, z] and study how the PSD depends on them. Our sample of combined light curves has only the r band, forcing us to select quasars in narrow redshift bins to not introduce a potential bias in the ensemble PSD due to mixing of sources with a significant large range of rest-frame wavelengths. This, together with the smaller sample of 4037 quasars, reduces our ability to study the PSD dependence on a broad range of BH masses and Eddington ratios because the number of adjacent [M_BH, λ_Edd, z] bins with enough sources to constrain their ensemble PSD accurately is more limited. Nonetheless, for each bin [M_BH, λ_Edd, z], we have a wider frequency range to constrain the PSD shape, while P24 had only three frequencies with Δ(log ν)∼0.5 day⁻¹.

To probe the PSD shape, we initially focussed on a single subset of 69 quasars and the following physical parameters: 8.7 < log(M_bh < 9.0, 45.9 < log(L_bol) < 46.1, and 1.25 < z < 1.45 (hereafter the master sample). This range of values corresponds to the peak of the distribution of the respective parameter values (see Fig. 1); hence, the number of quasars is the highest in the narrowest bin and their ensemble power spectrum is ideal for studying its shape in detail. Given the mean redshift of the quasars in this sample, the observed light curves can probe intrinsic quasar variations in the UV band (i.e. around ∼2600 Å). Other sources of the larger sample are used in Sect. 5 to discuss the dependence of the PSD on BH mass and accretion rate.

The combined PSD for the master sample is shown in Fig. 4, with different colours and symbols for the different contributions. The points on the left side of the plot show the low-frequency PSD, calculated on timescales of the order of years, computed with the SDSS, Gaia, and long-term ZTF light curves. For each group of these light curves, we computed the 69 periodograms, we subtracted the Poisson noise, and calculated their mean and its error, in each frequency bin. For the SDSS and the long-term ZTF light curves the Poisson noise was directly measured by computing the variance of light curves from non-variable stars as a function of magnitude, as described in Appendix A of P24. For Gaia PSDs, where light curves from non-variable sources are not available, we computed the theoretical Poisson noise contribution as the mean error squared of the light curve points (which is representative of the light curve variance due to Poisson noise) divided by the frequency range over which we calculate the power spectrum and the light curve mean squared (since this is how we normalised the periodogram; see Eq. 1).

Fig. 4. Ensemble power spectrum for the master sample of the quasars. Different colours and symbols show PSD estimates calculated using various light curves, as labelled (see Sect. 4 for details).

On the right side of Fig. 4, blue circles show the high frequency part of the ensemble PSD, calculated from the seasonal segments of the 12-day binned ZTF light curves. These points sample the power spectrum on time scales of weeks to months. To compute the ensemble PSD at these frequencies, we merged all the PSDs calculated from the approximately five seasonal PSDs of all the quasars in the master sample and computed the mean PSD (and its error) over bins of 100 points. The Poisson noise contribution was not subtracted from the high-frequency PSD estimates; instead, it was fitted a posteriori by adding a constant to the models we fit to the ensemble PSD.

The power spectrum estimates of the SDSS and the long-term ZTF light curves are consistent within the errors, as shown by the respective data points in the left part of Fig. 4. This is also the case with the Gaia PSD estimates, despite the fact that Gaia measurements include contributions from a wide range of wavelengths (due to the width of its filter). According to P24, the UV/optical power spectra of the quasars depend on the wavelength, and this could affect the Gaia measurements, although Fig. 4 shows that this is not the case. The agreement of the ensemble power spectra estimated from light curves taken years apart and with different instruments implies that the rest-frame UV variability process remain constant over the timescales explored in this work (∼20 years). This is an important result and provides a strong indication that the long-term variability in AGN is stationary⁴ on timescales of several years, although we cannot exclude non-stationarity at much longer timescales. The second result is that the variability amplitude continues to increase with increasing time scales, without clear evidence of a complete flattening of the PSD (to zero slope) down to about 10⁻³ days⁻¹ (i.e. ∼3 years in rest-frame).

To quantitatively determine the PSD shape, we fit the ensemble PSD of the master sample with a single power law, a bending power law, and the PSD associated with the DRW model, as we explain below. The single power law (PL) is defined as

$$\begin{array}{c}\hfill P(\nu )=PS{D}_{\mathit{norm}}{(\nu /{\nu }_0)}^{\alpha }+noise,\end{array}$$(2)

where ν₀ = 0.002 days⁻¹ is the frequency at which we evaluate PSD normalisation, PSD_norm. The more general bending power law (BPL) model takes the form

$$\begin{array}{c}\hfill P(\nu )=PS{D}_{\mathrm{norm}}{\nu }^{{\alpha }_{\mathrm{low}}}[1+(\frac{\nu }{{\nu }_b}{)}^{{\alpha }_{\mathrm{low}}-{\alpha }_{\mathrm{high}}}{]}^{-1}+\mathrm{noise}.\end{array}$$(3)

This model was introduced by McHardy et al. (2004). It describes the broad-band (i.e. low- and high-frequency) power spectrum of AGN in X-rays well. In the above equation, ν_b is the bending/break frequency above which the PSD slope steepens from α_low at ν ≪ ν_b, to α_high at ν ≫ ν_b. PSD_norm is equal to 2PSD(ν_b)×ν_b, and equal to PSD(ν)×ν at frequencies ν ≪ ν_b, for the BPL model with α_low = −1. In this case, when the power spectrum is plotted in the PSD(ν)×ν representation, it appears flat below ν_b. Its value is equal to PSD_norm, and for that reason, PSD_norm is frequently called the ‘power spectrum amplitude’, in the case of BPL models with α_low = −1. The power spectrum of the DRW model is equivalent to the BPL model with α_low = 0 and α_high = −2.

We fit the PL, BPL, and DRW models to the ensemble PSD of the master sample plotted in Fig. 4. We kept α_low = −1 fixed during the model fitting process of the BPL. If left free to vary during the fits, α_low is unconstrained, with best-fit values close to –1, but large errors due to the fact that the ensemble PSD does not extend to low enough frequencies. The best-fit results for the master sample are listed in Table 1, while the best-fit models and residuals are plotted in Fig. 5.

Fig. 5. Top panel: Best-fit PL, DRW, and BPL models to the ensemble PSD of the master sample. Horizontal lines mark the noise floor level for each model. Bottom panels: Respective best-fit residuals plots; i.e. observed PSD–model+error.

All models provide statistically acceptable fits to the data, as p_null > 0.01 for all of them. However, the improvement of the BPL fit to the data is significant when compared with the PL best-fit. The difference in the best fit χ² of the two models is Δχ² = 11.9 for 1 extra degree of freedom. This is statistically significant according to the F-test (F statistic of 12.4, p_null = 0.002). This result is strengthened and confirmed by the poor fits the PL model provides to the ensemble PSDs of quasars across various BH mass and accretion rate bins that are analysed in Sect. 5. We therefore conclude that the intrinsic PSD in the UV band of quasars with log(M_BH) ∼8.85 and log(L_bol)∼46 shows a significant flattening at low frequencies, below ν_b ∼ 1 − 3 × 10⁻³ day⁻¹.

BPL provides the best fit to the PSD. The best-fit high-frequency slope is −2.60 ± 0.25 which implies a difference ∼2.4σ compared to the slope of ∼ − 2 predicted by the DRW model. However, the goodness-of-fit of the BPL model is not statistically significant when compared with the DRW best-fit to the ensemble PSD of the master sample.

5. The dependence of the quasar PSD on BH mass and accretion rate

As in P24, we also investigated whether the power spectrum varies with M_BH and/or λ_Edd. To this end, we created six additional samples of quasars. Three of the samples include quasars with the same accretion rate (–1.0 < log(λ_Edd) < –0.7) and redshift (1.15 < z < 1.45), but different BH mass masses, BHM1, BHM2, and BHM3, including quasars with 8.4 ≤ log(M_BH) < 8.7 (61 quasars), 8.7 ≤ log(M_BH) < 9.0 (113 quasars), and 9.0 ≤ log(M_BH) ≥ 9.30 (63 quasars), respectively. The other three samples include objects with the same BH mass (8.7 < log(M_BH) < 9.0) and redshift (1.15 < z < 1.45), but different accretion rates, AR1, AR2, and AR3, including quasars with –1.30 ≤ log(λ_Edd) < –1.00 (69 quasars), –1.00 ≤ log(λ_Edd) < –0.75 (96 quasars), and –0.75 ≤ log(λ_Edd) ≥ –0.40 (66 quasars). We note that both the BHM2 and the AR2 samples are similar to the master sample in terms of the range of BH mass and accretion rates of quasars in them. We cannot calculate the ensemble PSD of quasars in additional bins to increase the range of the BH mass or the accretion rate values because the number of sources in each bin becomes too small and their distribution will not be Gaussian (Papadakis & Lawrence 1993). We can increase the redshift bin size beyond 0.3, but this will result in considering quasars with rest-frame light curves in significantly different energy bands. Since the quasar variability depends on energy (see P24), the ensemble PSDs will be increasingly difficult to interpret.

Figures 6 and 7 show the ensemble power spectra (calculated as before for the master sample) for the quasars in the various BHM and AR samples. Clearly, the PSDs are not identical, which implies that their properties must vary with both the BH mass and the accretion rate.

Fig. 6. Top panel: Ensemble power spectra for the quasars in the BHM1, BHM2, and BHM3 samples (mean values of the BH mass in each sample are reported in the legend). Dashed, solid, and dotted lines indicate the the best-fit BPL models to the power spectra. Thinner lines of the same colour and style are noise-subtracted. Bottom panels: Best-fit residuals.

Fig. 7. Same as in Fig. 6, but for the AR1, AR2, and AR3 samples (mean values of the accretion rate for each sample are reported in the legend).

We fit the ensemble PSDs of the BHM and AR quasar samples with the DRW and BPL models. We did not consider the PL model, as it gives very poor fits to most of the data, which is in agreement with the result in the previous section. The best-fit results are listed in Tables 2 and 3. Statistically speaking, both the DRW and the BPL models fit the PSDs well (i.e. p_null is larger than 0.01 in all cases). The only exception is the BHM3 PSD, where p_null is low for both best-fit models. The respective residuals plot in Fig. 6 shows that the BPL model fits well the overall PSD shape, but the PSD is rather noisy, hence, the low p_nul.

BPL provides a better fit to all ensemble power spectra, for both the BHM and AR samples. Based on the assumption that there is an underlying universal model that has to fit the PSDs in all bins, the total χ² for the best-model fits in the case of the DRW model is 183.8 for 129 degrees of freedom. The respective numbers are χ_total² = 152.2 for 124 degrees of freedom in the case of the BPL model (excluding the best-fit models to the BHM3 PSD for both models). The goodness of the DRW fit to all PSDs is poor (p_null = 0.0012), while the BPL provides an acceptable fit to all PSDs (p_null = 0.02). Furthermore, the improvement in the goodness of fit in the case of the BPL model (Δχ² = 24.5/5 dof) is statistically significant according to the F-test (F-statistic 3.8, p_null = 0.0029). We therefore conclude that, in general, the quasar power spectra are not consistent with the predictions of the DRW model.

The ensemble PSDs plotted in Figs. 6 and 7 can be used to investigate the dependence of the PSD characteristics (i.e α_high, PSD_norm and ν_b) on the physical parameters of the quasars, that is, M_BH and λ_Edd. We find that α_high does not depend on either M_BH or λ_Edd. All the best-fit α_high values listed in Tables 2 and 3 remain consistent with the high-frequency slope of the ensemble PSD listed in Table 1. The weighted mean of α_high in the master, BHM1, BHM3, AR1, and AR3 samples (which include objects which do not appear in the other samples) is −2.7 ± 0.1 and this is our best estimate of the high-frequency slope in quasars (in the UV band).

All the best-fit values of PSD_norm listed in Table 2 are consistent with the best-fit PSD_norm of the master sample, implying that PSD_norm does not depend on M_BH. However, the best fit value of PSD_norm for the AR3 sample [log(λ_Edd)∼ − 0.64] is approximately half that of PSD_norm in the other two AR samples listed in Table 3, and the difference is highly significant. These results indicate a possible decrease in PSD amplitude at accretion rates greater than 10^−0.64 (i.e. ∼0.23 of the Eddington limit).

Finally, we also considered the dependence of the bending frequency, ν_b, on the BH mass and accretion rate. Regarding the latter parameter, we did not notice any trend between ν_b and λ_Edd. The best-fit values ν_b in Table 3 are consistent (within the errors) with the best-fit bending frequency in the PSD of the master sample. However, we noticed a trend of decreasing bending frequency with increasing BH mass, in Table 2. A straight line (in the log–log space) with a slope of −0.6 ± 0.1, fits the ν_b vs M_BH plot well. This implies an anti-correlation between the bending frequency and BH mass, although it is based on only three points.

However, if PSD_norm and α_high are not dependent on M_BH, then the relation ν_b∝ M${}_{\mathrm{BH}}^{-0.6}$ we find can explain the difference between the PSDs shown in Fig. 6. The PSD amplitude at low frequencies [in PSD(ν)×ν] could be the same in all objects, but because the bending frequency decreases with increasing BH mass, the PSD amplitude in the sampled frequency range appears to decrease with increasing BH mass. In the case of the PSDs in Fig. 7, the main difference is between the λ_Edd ≈0.23 PSD (green triangles) and the other two PSDs. The lower amplitude of this PSD could, in fact, be due to a difference in PSD_norm, which can decrease with increasing accretion rate. The ratio between the normalisation in the first and last bin is PSD_norm(AR3)/PSD_norm(AR1) ∼0.6 ± 0.1.

6. Discussion and conclusions

We computed the ensemble power spectrum of hundreds of quasars, which are part of the SDSS Stripe-82 sample of P24, in the rest-frame UV band and in various BH mass and accretion rate bins. We were able to estimate the power spectrum over two orders of magnitude in frequency, from ∼1/(10 days), to ∼1/(1000 days). This is a significant advance compared to our previous work (P24), where we were able to estimate the quasar power spectra in just three frequencies, over a narrower range. The advance was made possible because we considered data from two additional surveys in addition to SDSS (ZTF and Gaia), allowing us to extend the PSD at high frequencies and better constrain it at low frequencies. We estimated the power spectrum using traditional Fourier techniques. As a result, the PSD units are well defined, and the UV PSDs can be directly compared with AGN PSDs estimated in other energy bands (i.e. X-rays). However, cross-matching data from multiple surveys resulted in limiting the original P24 sample to half of the sources and a single photometric band (to limit cross-calibration biases in combining the light curves; see Appendix A). As UV/optical variability of AGN depends on energy (e.g. P24), we restricted the analysis to a narrow redshift range 1.25 < z < 1.45, resulting in PSDs at ∼2600Å rest-frame (i.e. the UV band). While P24 was able to exploit a larger sample and five filters ugriz to get many [M_BH, λ_Edd] bins at any given wavelength, this work focuses on a limited number of bins spanning ≲1 dex in each parameter. The number of sources in other bins was significantly smaller than the minimum necessary for the statistical properties of the ensemble PSD estimates to be appropriate for model fitting. Nonetheless, for each bin we have a wider frequency range to better constraint the intrinsic PSD shape. Our main results are summarised below.

1. We were able to accurately estimate the ensemble PSD of quasars in three BH mass bins with an average M_BH = 3.8 × 10⁸, 7.1 × 10⁸, and 1.3 × 10⁹ M_⊙, and with three accretion rate bins with an average λ_Edd = 0.08, 0.13, and 0.23.

2. We find that a power-law model does not fit the PSDs well. We detected bending frequencies, ν_b, in the ensemble PSDs in all BH mass and accretion rate bins. This is the first time that bending frequencies have been detected directly in the UV power spectra of a single sample of luminous quasars, without the need to combine PSDs calculated on narrower frequency ranges and then rescaling frequencies according to pre-defined relations between ν_b and BH mass.

3. A bending power-law model fits all PSDs well, with a mean slope at high frequencies, above ν_b, of −2.7 ± 0.1 depending neither on M_BH nor on λ_Edd. The slope at low frequencies is consistent with –1 in all ensemble PSDs (irrespective of the BH mass and/or accretion rate).

4. The DRW model does not fit the ensemble PSDs equally well. A low frequency slope of zero (as expected in the case of the DRW model) is not consistent with the data (as the best-fit results worsen considerably, when compared with the BPL best-fit models). Furthermore, the high-frequency BPL best-fit slope of −2.7 ± 0.1 is significantly different from the DRW prediction of −2.

5. The amplitude of the BPL model, PSD_norm, does not depend on M_BH, although it might decrease with increasing accretion rate. The ratio between the normalisation in the first and last bin, PSD_norm(AR3)/PSD_norm(AR1), is ∼0.6 ± 0.1. The possible existence of a trend is tentative and we need to study ensemble PSDs of quasars on a wider range of accretion rates to reach conclusive results.

6. The bending frequency decreases with BH mass, approximately as ${\nu }_b\propto M_{\mathit{BH}}^{-0.6\pm 0.1}$. We did not detect any dependence of ν_b on the accretion rate but, as in the case of PSD_norm, we need to fit the ensemble PSDs of quasars over a wider range of accretion rates to reach conclusive results.

Figure 8 shows a schematic representation of the quasar’s power spectrum in the UV band, based on our results. The power spectrum in this figure is plotted in the PSD(ν)×ν representation. The low-frequency PSD slope was fixed at –1. If we leave α_low free to vary during model fitting, the best-fit results are still consistent with –1, but the errors are large (not only in α_low but also in the other parameters). The dependence of the bending frequency on BH mass, as well as the possible dependence of the PSD amplitude on the accretion rate, are based on the results of fitting the ensemble PSDs over three bins covering less than 1 dex in BH mass and λ_Edd. Clearly, it is important to extend the study of the ensemble power spectrum of quasars to many more BH mass and accretion rate bins to better constraints the scaling relations. An important factor that determines the size of the sample and the number of quasars in narrow BH mass and accretion rate bins is the use of the Gaia light curves. The future Gaia data release will hopefully include more quasars and additional epochs overlapping with the beginning of ZTF. This, together with efforts in populating the gaps between the surveys and/or investigating sophisticated interpolation methods, could be used to (reliably) produce evenly spaced time series using as input light curves such as the one plotted in Fig. 2 and to study the PSD on an even larger frequency range. We discuss this topic further in this section. In the future, we also plan to explore other methods to estimate the ensemble power of quasars, which will require fewer objects per bin. Despite the cautious remarks above, Fig. 8 could be used to gain important information on the UV variability of quasars.

Fig. 8. Schematic representation of the best-fit bending power-law PSD model and its scaling with BH mass and accretion rate.

6.1. Comparison with X-rays and additional physical constraints

The best-fit PSD normalisation when we fit the ensemble PSD of the master sample with the BPL model is 0.011± 0.002 (day)⁻¹. Since a_low = −1 in Eq. 3, then PSD(ν)×ν will be equal to 0.011 at very low frequencies (i.e. when ν < < ν_b). Interestingly, this is also the case with the X-ray PSDs of AGNs. The results indicate that the X-ray PSD normalisation is roughly the same in all AGNs we have studied so far, and that PSD_{X − rays}(ν)×ν ∼ 0.01 − 0.02 at frequencies much lower than the break frequency. For example, Paolillo et al. (2023) assumed a BPL model (with a_low = −1) for the X-ray PSD and found that A = 0.016 ± 0.003 Hz⁻¹ for a large sample of AGNs, which also includes very bright quasars. This implies that PSD_{X − rays}(ν)×ν ∼ 0.016 at frequencies lower than ν_b. Furthermore, Papadakis & Binas-Valavanis (2024) studied the 14–195 keV band power spectrum of the brightest AGN in the BAT survey sample. They found that the PSD is consistent with a power-law model with a slope of −1 at low frequencies (i.e. at frequencies lower than the expected break frequency) and they also found that PSD_{X − ray}(ν)×ν = 0.014 ± 0.003 at these low frequencies. These results indicate that when normalised to the (square) of the mean flux, the PSD amplitude at frequencies much lower than the break frequency (i.e PSD(ν)×ν) is ∼0.01 − 0.015, both in X-rays and in the UV band. This is an interesting result. It may be a coincidence, but it is also consistent with X-ray reverberation. Papoutsis et al. (2025) showed that the observed UV/optical variability in SDSS quasars, from 1300 Å to 4000 Å, could be the result of X-ray reverberation, if the X-ray corona is powered by the accretion power, the black hole spin is lower than ∼0.7, and the height of the X-ray corona is in the range of 20–40R_g, where R_g = GM_BH/c² is the gravitational radius of the BH.

The best-fit bending frequencies in the UV power spectra can also be used to constrain various models. In the case of X-ray reverberation, we expect two breaks in the UV (and optical) PSD (see Papoutsis et al. 2025). The first is the bending frequency that also appears in the X-ray PSD. However, the bending timescale, T_b, UV = 1/ν_b, is ∼285 days for the master sample (see Table 1). The respective timescale in X-rays for a quasar with M_BH ∼7 × 10⁸ M_⊙ and λ_Edd ∼0.13, is T_{b, X − rays} ∼ 22 days, assuming the McHardy et al. (2006) relation. Therefore, the bending frequencies we detect cannot be associated with the bending frequency in the X-ray PSD.

The second bending frequency in the UV/optical PSDs should be due to the break frequency in the disc transfer function. Roughly speaking, this bending time scale is associated with the time it takes for the X-rays to propagate from the inner disc up to the outer radius of the part of the disc that emits in the UV. Interestingly, Panagiotou et al. (2022) predicted that this PSD break should be proportional to M${}_{\mathrm{BH}}^{-0.65}$ (see their Eq. A2), which is very similar to what we found. This result supports the X-ray reverberation hypothesis.

However, the bending timescales that we detect in the UV PSDs are rather low. For example, Papoutsis et al. (2025) showed that the bending frequency in the disc transfer function should be around ∼10⁻² day⁻¹ for an AGN with a mass of 8 × 10⁸ and an accretion rate of 0.1 of the Eddington limit (see Fig. 1 in their paper). However, the best fit ν_b in the combined PSD of the master sample is about two to three times smaller. We need PSDs for many more objects, preferably sampled at even lower frequencies, to perform tests with theoretical models.

The other possibility is that the observed optical/UV variations are due to intrinsic disc variations. Lyubarskii (1997) proposed a detailed model with respect to the variability of the flux of the disc. They showed that the resulting PSD will have a slope of –1 if there are fluctuations of the accretion rate that propagate inward, and their characteristic timescale is on the order of the local viscous timescale. In this case, the bending timescale should correspond to the viscous time at the inner disc radius. If that is equal to 6R_g (in the case of BH with spin zero), then the expected viscous time scale for a ∼7 × 10⁸ M_⊙ is ∼600 days, assuming that the radius vs disc height ratio is 10, and the viscosity parameter is 0.1 (Paolillo & Papadakis 2025). This is twice the timescales we observed in our work. Perhaps a different viscosity parameter and/or disc radius-to-height ratio could explain the observed discrepancy. Another possibility could be that the flux of the disc varies on the thermal time scale. For the PSD of the master sample (with M_BH ∼ 7 × 10⁸ M_⊙), the break time scale of ∼300 days in the PSD could be equal to the thermal time scale at ∼80R_g (assuming that the viscosity parameter is 0.1). It is not clear why the disc should be variable on the thermal time scale at this radius. In any case, a model should not only match the observed time scale with some theoretical timescale, but should also be able to explain the amplitude and slope of the PSD both at low and high frequencies.

More recently, Sun et al. (2020) proposed a Corona-Heated Accretion-disk Reprocessing (CHAR) model, which fits the data with steeper than DRW slopes at high-frequencies, and no strong evidence for a complete flattening with light curves over than 20-years long (Zhou et al. 2024). Another possibility is that variability at low and high frequencies are associated with two different components, accretion rate fluctuations and X-ray reprocessing, respectively (e.g. Yu et al. 2025).

6.2. Comparison with P24 and other past results

P24 estimated the ensemble PSD of SDSS quasars over many bins of BH mass, accretion rate and rest-frame wavelength, collectively covering a frequency range from 250 to 1500 days (rest-frame). However, as the SDSS light curves used for the analysis were six years long, they were limited to three frequencies per bin and unable to directly detect any bending time scale in the power spectrum. They found that the PSD over the SDSS range could be fitted well with a PL model and that the amplitude and slope of the best-fit PL models depend on M_BH and λ_Edd, as follows: log(PSD_amp, P24)∝ − 0.5log(M_BH)−0.7log(λ_Edd), and PSD_slope, P24 ∝ −0.5log(M_BH)−0.4log(λ_Edd). Our results, taking advantage of the additional frequencies, put better constraints on the intrinsic PSD shape and can be used to understand the P24 relations.

We plotted the P24 PSD parameters (i.e. log(PSD_amp, P24) and PSD_slope, P24 with the best-fit PSD parameters we determined (i.e. ν_b, and PSD_norm, BPL). We find that PSD_amp, P24 ∝ ν_b^0.7 and PSD_slope, P24 ∝ 0.7log(ν_b). Since ν_b∝ M${}_{\mathrm{BH}}^{-0.6\pm 0.1}$ (see Sect. 5), then we can easily explain the dependence of PSD_amp, P24 and PSD_slope, P24 on BH mass as the result of ν_b shifting across the narrow frequency range sampled by P24 (as also discussed in their Fig. 19). The dependence of PSD_amp, P24 and PSD_slope, P24 on the accretion rate could be due to a possible dependence of ν_b on λ_Edd. As discussed in Sect. 5, we do not detect a correlation between ν_b and λ_Edd, but this may be due to the limited range and number of accretion rate bins in which we could estimate ensemble PSDs. It is also possible that the dependence of PSD_amp, P24 and PSD_slope on the accretion rate may be due to the decrease of PSD_norm with increasing λ_Edd. However, we cannot confirm a correlation between PSD_amp, P24 and PSD_slope, P24 with PSD_norm, BPL in the limited range of λ_Edd we considered in this work.

The relation between break timescales and BH mass is found in many recent works, using independent methods and datasets (e.g. Tang et al. 2023; Arévalo et al. 2024; Yuk & Dai 2025), although there is no clear consensus on the exact dependence as the comparison between explicit PSD calculation, CARMA model fitting or SF analysis is not trivial. Burke et al. (2021) fit the light curves of ∼70 quasars with a DRW model, and found that the characteristic bending timescale is related to the BH mass as τ∝ M${}_{\mathrm{BH}}^{0.4\pm 0.05}$. This relation is flatter than the one we determined here (${\nu }_b\propto M_{\mathit{BH}}^{-0.6\pm 0.1}$). Furthermore, we used their Eq. (1) and computed the expected ν_b for quasars with the mean BH mass in the BHM1, BHM2, and BHM3 samples. We found that the Burke et al. (2021) predictions for ν_b are systematically ∼1.6 times smaller than our best-fit values. The (significant) discrepancy between our results and those of Burke et al. (2021) is probably due to the fact that these authors assumed that the light curves are well fitted by a DRW model and computed the associated timescale, while we fit the ensemble PSD with a general BPL.

A work similar to P24 was also presented by Arévalo et al. (2024), who used g-band light curves obtained by ZTF through their public data release DR14, and considered a large number of quasars in the small redshift range of 0.6 < z < 0.7, over bins of BH mass and accretion rate. Their redshifts imply a rest-frame wavelength of ∼2900Å, which is close to the one presented in this work. They initially assumed the Burke et al. (2021) correlation between timescale and BH mass to shift the frequency of individual PSDs over a broader frequency range. Then, they fit all PSDs with BPL models at various fixed slopes, allowing for a dependence on both BH mass and accretion rate. The final best-fit has low- and high-frequency slopes of –1 and –3, respectively, with ${\nu }_b\propto M_{\mathit{BH}}^{-0.67to-0.5}\times {\lambda }_{\mathrm{Edd}}^{-0.41to-0.12}$. While we are in agreement with the dependence of ν_b on BH mass, we do not find and trend with the accretion rate and our mean high-frequency slope of –2.7 is a bit flatter (although –2.5 is their second-best fit, and their procedure is fixing the slope a priori while to us is a free parameter). We also find lower PSD amplitudes, which may be due to the different technique to calculate the PSD, and a stronger dependence of PSD_norm on the accretion rate [PSD_norm(AR3)/PSD_norm(AR1) ∼0.6 ± 0.1, while according to Arévalo et al. (2024) it should be ∼0.3, using their Eqs. 3 and 4, the values listed in their Table 3, and the mean accretion rate and the BH mass of the quasars in the AR3 and AR1 bins].

We believe the main reason for the differences between the Arévalo et al. (2024) results and our work stems from the use of the model-dependent assumption adopted to shift PSDs. This clearly shows that it is important to estimate the power spectra of quasars over the broadest possible frequency range to detect the bending frequency and the low-frequency slope with high accuracy. At the same time, a proper characterisation of the PSD dependence on the AGN physical properties requires having a large sample of sources to be binned in many subsets. While analyses such as P24 and Arévalo et al. (2024) definitely satisfied the second requirement, the present work is aimed at constraining the PSD shape of fewer quasar samples with the most model-independent approach as possible. As the study of AGN variability becomes more and more detailed, the general agreement is that DRW can be used as a first order approximation, but it does not accurately describe the results. Several works report of high-frequency PSD slopes steeper than the –2 value of the DRW, consistent with our findings, along with possible evidences of multiple break frequencies (e.g. Mushotzky et al. 2011; Tachibana et al. 2020; Stone et al. 2022). Analogous conclusions have also been reported by Yu et al. (2025), who analysed a sample of combined AGN light curves similar to ours, showing that a damped harmonic oscillator (i.e. a higher order CARMA model) fits the data significantly better than a DRW.

6.3. Very-low-frequency PSD from combined light curves: Attempts at interpolation

The only way to determine the PSD at frequencies lower than the the range of this work is to combine light curves from as many different surveys as possible and compute a unique PSD from them. This requires a considerable effort in calibration and homogenisation, but also accounting for the gaps between the surveys with no overlapping region. In the light curves we initially considered, albeit photometry was cross-calibrated (Appendix A), the extent of such gaps and their number prevent us from using a simple linear interpolation to account for the missing points, as we did for the high-frequency ZTF data. Several advanced techniques have been proposed to deal with non-uniform time series (e.g. Lefkir et al. 2025), but the biggest challenge is to avoid (or limit) the addition of spurious variability and keep the analysis model-independent. Moreover, when combining data from different surveys, ensuring there are no residual offsets between them is fundamental and requires applying the same approach on non-variable objects (which we could not do with Gaia DR3).

A possibility is to model the light curve with a continuous function based on the measured points, and then use yearly-spaced epochs to measure the PSD. We tested a Radial Basis Function (RBF) interpolation method, which is commonly used to get a continuous approximation of a dataset from scattered and irregular points, such as uneven time series. More details about the implementation, the choice of parameters and its limitations are provided in the Appendix C. We focused on the master sample of 69 quasars defined in Sect. 4 and used the full r-band light curves, including the early SDSS epochs and the PS1 data, as in Fig. 2. The RBF-interpolation procedure produced evenly sampled light curves with a 356-day cadence, which we used to compute the PSD as described in Sect. 3. While at high frequencies the reconstructed PSD-shape from the interpolation is strongly dependent on the input parameters (i.e. kernel function, smoothing, polynomial degree, sampling cadence), at frequencies ν ≲ 2 × 10⁻³ day⁻¹ the outcomes are more robust for any reasonable input combination.

The results of the interpolated PSD are shown in Fig. 9, as purple points on top of the other independent estimates from the single surveys. There is a good agreement in the overlapping region, of around ν ∼ 10⁻³ day⁻¹. The key improvement is that the new light curves allow for the calculation of the PSD at frequencies that are about ten times lower than the lowest frequency in the PSDs we studied. The PSD at these frequencies appears to still increase, but at a lower rate. If we fit the new PSD with the same BPL model as before, we found best-fit parameters that are consistent within the errors with those of Fig. 5. However, the new, low-frequency points appear to have a higher normalisation than the other PSD points. This is a potential issue, as it is possible that creating an evenly sampled light curve via RBF interpolation could introduce additional variance to the light curves. To have a fully reliable interpolation method, thorough simulations of all possible outcomes are required. The tests we performed on simulated DRW light curves show that the method is potentially able to reconstruct the mean long-term variability trends, but the reliability of the reconstructed ensemble PSD depends on the position of the break and the variability amplitude (see Appendix C). This is a critical point for all the possible light curve fitting or interpolation techniques, which can potentially bias the results of the analysis. Given these issues, we did not measure the ensemble PSD on the combined light curves, but we did choose to present our preliminary attempts to create evenly sampled light curves, combining data from various surveys as they might be useful for future works. In fact, this will be the only way to probe the quasar long-term variability, until LSST survey light curves are delivered a few years from now. However, even following the start of LSST, such approaches will be crucial to investigating AGN variability on timescales of over a decade, as well as for high-redshift samples where the probed rest-frame timescales are much shorter than those of local AGNs.

Fig. 9. Ensemble power spectrum for the master sample of quasars, as in Fig. 4, with the addition of estimates from the RBF-interpolated light curve combining all the surveys (purple stars).

Acknowledgments

MP, VP, DD and MF acknowledge the financial contribution from PRIN-MIUR 2022 and from the Timedomes grant within the “INAF 2023 Finanziamento della Ricerca Fondamentale”. IEP acknowledges support from the Visiting Professor programme of the Federico II University. DD also acknowledges PON R&I 2021, CUP E65F21002880003. CMR and MIC acknowledge financial support from the INAF Fundamental Research Funding Call 2023 (PI: Raiteri). 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: Cross-calibration of light curves

The main challenges when comparing data from different surveys arise from differences in the the detector properties, photometric systems, and the overall data reduction process. The normalised filter transmission curves for SDSS, PS1, and ZTF gri bands are reported in Fig. 1 of Ngeow et al. (2019), whereas Fig. 3 of Jordi et al. (2010) shows the Gaia broad filters. It is clearly evident that while the filters for the ground-based surveys can be adapted via relatively small colour corrections, the Gaia wide field photometry requires more dedicated attention. Indeed, by following Jordi et al. (2010) and the DR3 documentation⁵, it is possible to use three wide bands to match the SDSS system with a 4th order polynomial expression:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill G-r^{\mathit{SDSS}}& =-0.09837+0.08592(G_{\mathit{BP}}-G_{\mathit{RP}})+\hfill \\ \hfill [0.5em]& 0.1907{(G_{\mathit{BP}}-G_{\mathit{RP}})}^2-0.1701{(G_{\mathit{BP}}-G_{\mathit{RP}})}^3\hfill \\ \hfill [0.5em]& +0.02263{(G_{\mathit{BP}}-G_{\mathit{RP}})}^4.\hfill \end{array}\end{array}$$(A.1)

A similar equation can be used to go from G to g^SDSS, but starting from the the same broad band photometry means that it would not add information to the data, as they would be nothing but shifted version of the same magnitude. Hence we can only use the broad filter from Gaia calibrated to either the g or r SDSS band.

The ZTF Data Products also come with filter-based zero-points (ZP_f) and colour-coefficients (c_f) for photometric calibration with the PS1 system (Masci et al. 2019). The colour correction is defined as polynomial equations:

$$\begin{array}{c}\hfill \begin{array}{c}{g}^{PS1}-g^{\mathit{ZTF}}=Z{P}_g+c_g(g^{PS1}-r^{PS1}),\hfill \\ [1em]r^{PS1}-r^{\mathit{ZTF}}=Z{P}_r+c_r(g^{PS1}-r^{PS1}),\hfill \end{array}\end{array}$$(A.2)

where the ZP and c coefficients are provided in ZTF DR19, for each epoch.

The final step in colour correction is to compare the SDSS and PS1 filters, in order to correct any offsets and standardise all data to a common photometric system. To achieve this, we used a sample of ∼5 × 10⁵ non-variable Stripe-82 stars from Ivezić et al. (2007) to create colour-magnitude diagrams with PS1 and SDSS PSF photometry. We found that r band photometry is nearly equivalent, at a 1% level up to 20.5 mag, thus requiring no correction (as also reported in Suberlak et al. 2021). The g band, instead, presents an offset that we can correct with a linear fit,

$$\begin{array}{c}\hfill g^{PS1}-g^{\mathit{SDSS}}=-0.04-0.11(g^{PS1}-r^{PS1}).\end{array}$$(A.3)

The colour differences and best-fit results are shown in Fig. A.1.

As an additional check, we verified the validity of the corrections from ZTF to PS1 system through Eqs. A.2. Furthermore, we used light curves of non-variable stars to search for any residual biases. After a 5-σ clipping in magnitude and errors to limit the impact of bad photometry (mostly due to a few erratic PS1 epochs), deviations with respect to the mean value remain of the order of a few sigma. However, it is worth mentioning that the ZTF median error on magnitude is ten times larger than that of SDSS, which is expected because of the lower sensitivity. The same test could not be done for Gaia, due to the lack of light curves for non-variable sources. This is a critical point if our aim is to use the full combined light curve to extract the PSD, but future Gaia data releases could help resolve this issue.

To minimise the uncertainties due to photometric transformations, we limited our analysis to the r band for the combined light curves from the four surveys (following Suberlak et al. 2021). This way, SDSS and PS1 require no modification, while ZTF and Gaia are corrected as explained above. While the photometric uncertainties on the single light curve might prevent an accurate determination of variability parameters, our binning procedure over discrete time intervals and ensemble PSD analysis should mitigate the impact of small fluctuations, as shown in P24.

An example of cross-calibrated r-band light curve is shown in Fig. 2. Magnitudes are converted into fluxes in units of nJy, following standard procedure for error propagation and using the zero points associated with each survey.

Fig. A.1. Colour-magnitude diagrams with PS1 and SDSS magnitudes of non-variable stars. The points are colour-coded by density, with yellow indicating regions of higher point density, and darker areas representing lower-density regions. The upper panel shows the linear correction applied to the g band, while the lower panel demonstrates that the r filters are nearly identical, requiring no correction.

Appendix B: PSD of simulated light curves

Fig. B.1. Ensemble power spectra of simulated DRW light curves with the same observing pattern as the ZTF data. Red circles and squares show the measured PSDs at high and low frequencies, respectively. Solid black line refers to the power spectrum associated to a DRW model with σ = 0.28, and log τ = 3 (left panel), log τ = 3 (middle panel), and log τ = 1 (right panel).

To assess potential biases introduced by the linear interpolation in ZTF light curves for the high-frequency PSD, we tested it on simulated data for which we know the intrinsic PSD shape. We used the same approach as P24 (see their Appendix B) to simulate DRW light curves with any given amplitude, σ, and damping timescale, τ, following the procedure described in MacLeod et al. (2010) and Kovačević et al. (2021). For each fixed DRW τ and σ, we created 69 simulated light curves and we applied the same sampling pattern as the ZTF data from the master sample. We then process the data to measure the PSD exactly as we do in Sect. 3.1, that is, we bin each seasonal segment with a 12 days window, apply linear interpolation for eventual gaps, and compute the ensemble PSD. The requirements for the interpolation of missing data are very strict: no more than three consecutive gaps and no more than 20% made up of interpolated points. We stress again that the only reason behind the need for interpolation is that Fourier analysis requires evenly sampled light curves with no zeros or gaps. Otherwise, these zeros will be interpreted as real numbers, introducing very large amplitude features at high frequencies, which will be artificial and seriously affect the estimation of the true, intrinsic PSD of AGN.

Figure B.1 shows examples of power spectra derived from simulated ZTF-like light curves, together with the intrinsic DRW model. To further test the reliability of our analysis, we show both the high-frequency PSD (subject to linear interpolation) and the two low-frequency points estimated as in Sect. 3.2. For any fixed DRW variability amplitude, we explored different damping timescales to test the accuracy of the reconstructed PSD against different break positions. The figure clearly shows that the recovered PSD is in agreement with the input in any configuration. The slight mismatch at very high-frequencies is probably due to aliasing (i.e. higher non-sampled frequencies folded back to lower frequencies). This is a well know and expected effect, hence, we believe this small offset is not due to interpolation effects. In any case, in the real light curves, the frequency range where it appears is already dominated by the Poisson noise PSD component, which has a much larger amplitude than this mismatch, and it is not affect it by it. Therefore, we do not believe that the linear interpolation scheme we have adopted in this work is affecting the reliability of our PSD analysis.

Appendix C: RBF interpolation

In Sect. 6.3, we describe our test of a radial basis function (RBF) interpolation method, using the Python SciPy class RBFInterpolator⁶. RBFs are scalar functions whose value at any point x can be expressed in terms of the distance from the centre, r = ||x − c||, where c is the centre of the RBF. For applications with time series, points are the individual observations, while spatial distances are indeed temporal. RBF interpolation consists of expressing a target function f(x) as a linear combination of RBFs, each centred at one of the data points, and polynomials (Wahba 1990; Fasshauer & Zhang 2007). A proper combination of such functions can allow us to gather information about missing epochs, by following the overall trend of the light curve. One of the advantages of RBF interpolation over more advanced methods, like Gaussian processes (GPs), is its fast and straightforward implementation, with minimal assumptions on the intrinsic variability model.

More into detail, for a vector of fluxes f measured at times t, the RBF interpolant function at the new times $\stackrel{\sim }{\mathbf{t}}$ can be written as

$$\begin{array}{c}\hfill \mathbf{f}(\stackrel{\sim }{\mathbf{t}})=K(\stackrel{\sim }{\mathbf{t}},\mathbf{t})\mathbf{a}+P(\stackrel{\sim }{\mathbf{t}})\mathbf{b},\end{array}$$(C.1)

where $K(\stackrel{\sim }{\mathbf{t}},\mathbf{t})$ is a matrix of RBFs with centres at t evaluated at the times $\stackrel{\sim }{\mathbf{t}}$, and $P(\stackrel{\sim }{\mathbf{t}})$ is a matrix of polynomials of any degree. The use of polynomials is optional, and in some cases they may help to preserve global trends or satisfy boundary conditions. The coefficients a and b are determined by solving the following system of linear equations:

$$\begin{array}{c}\hfill (K(\mathbf{t},\mathbf{t})+\lambda I)\mathbf{a}+P(\mathbf{t})\mathbf{b}=\mathbf{f},\end{array}$$(C.2)

and

$$\begin{array}{c}\hfill P{(\mathbf{t})}^T\mathbf{a}=0,\end{array}$$(C.3)

where λ is a non-negative smoothing parameter that controls how closely the interpolant fits the data (exact fit for λ = 0), and I is the identity matrix. In the system, Eq. C.3 in the system ensures that the RBF terms are independent of the polynomial terms. In other words, it separates the local interpolation from any global trends the polynomial might capture, such as constant shifts or other complex behaviour in the data, while the RBFs handle local variations.

To measure the PSD on long timescales, we tested the RBF interpolation on the combined light curves with different choices of parameters. The best results are found with a multiquadric RBF kernel,

$$\begin{array}{c}\hfill \varphi (r)=\sqrt{1+(\epsilon r{)}^2},\end{array}$$(C.4)

where r is the distance and ϵ is the shape parameter controlling the width of the kernel. We used the default shape given by the median pairwise distance between the input times. As a smoothing factor we use λ = 0.8, which balances the fit to the data while avoiding over-fitting the noise. Adding high-order polynomials does not improve the PSD, so we use the default implementation with degree = 0.

Once the interpolant is built, we can evaluate it on a uniform grid,

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathbf{t}}=[t_{min},t_{min}+\mathrm{\Delta }t,\cdots ,t_{max}],\end{array}$$(C.5)

with a user-defined cadence. We tested values of 36, 180, or 360 days, with shorter cadences resulting in better reconstructed light curves. However, we estimate the final PSD always with a 360 days window, and frequencies shorter than 10⁻³day⁻¹ are not affected by the RBF interpolation cadence. For each interpolated point, flux errors are calculated using a linear interpolation of the original errors. Additionally, we use a running sigma-clipping over a 36-day window on the light curve to reduce biases from bad photometry, and inject Gaussian noise to the interpolated points to account for measurement uncertainties.

An example of combined light curve, and the corresponding RBF interpolation with a 36-day cadence is shown in Fig. C.1. While the interpolation may be not reliable to reconstruct the short-terms variability (strongly depending on the RBF parameters and the distribution of points), it might be used to fill sporadic gaps on timescales longer than the typical length of an observing season. We used interpolated light curves as described here to measure the PSD at short frequencies shown in Fig. 9.

To verify the reliability of the method in reconstructing the power spectrum, we tested it on simulated light curves as described in Appendix B. For each fixed DRW τ and σ, we created 69 simulated light curves and we applied the same sampling pattern as that of the combined SDSS+PS1+Gaia+ZTF data from the master sample. We then interpolated the data with the RBF method and computed the power spectrum as we do for the real data. Figure C.2 shows examples of power spectra derived from interpolated light curves, together with the intrinsic DRW model. For a fixed variability amplitude, we explored damping timescales in the range [10¹ − 10⁴] days, to probe regions where the interpolated points are before, within or after the break. The low-frequency part of the reconstructed PSD is generally compatible with the model in the power-law and the early-break regimes (first four panels). However, in the last two cases where the intrinsic PSD flattens, the interpolated light curves show residual power in disagreement with the model. Moreover, also in the cases where they seem to be in agreement, spurious breaks appear in the interpolated power spectra. This clearly shows that our RBF-interpolation cannot be used as a universal interpolation or the reconstruct the unknown PSD shape of a stochastic process, as the outcome depends on the underlying process. For this reasons, we do not include the interpolated points in the analysis, albeit we show here our preliminary attempts and the drawbacks to the interested reader. More extended simulations with other RBF parameters and additional interpolation methods are required, to fully exploit the potential of combined light curves.

Fig. C.1. Example of quasar r-band light curve combining all the surveys (black points), with the RBF interpolated version over a 36-days window (green points).

Fig. C.2. Ensemble power spectra of simulated DRW light curves with the same observing pattern as combined SDSS+PS1+Gaia+ZTF data, interpolated via RBFs. Red points show the measured power spectral densities, while the black solid line refers to the power spectrum associated to the DRW model. Each panel has DRW light curves with different damping timescales, as listed on the figures, probe regions where the interpolation falls before, within or after the break.

All Tables

All Figures

	Fig. 1. Distributions of bolometric luminosity (left panel) and BH mass (right panel) for the 8042 SDSS+PS1+ZTF quasars (red) and the subset of 4037 sources in common with Gaia DR3 (blue).
In the text

	Fig. 2. Example of r-band light curve for a Stripe-82 quasar combining SDSS (blue), PS1 (green), Gaia (orange), and ZTF (red) data. Colour corrections and calibrations are given in Appendix A.
In the text

Fig. 3. Example light curve for a ZTF quasar in the r band. Grey points are the original flux measurements, while yellow stars indicate the binned-interpolated values turning each season into an evenly sampled time series to compute the high-frequency PSD in Sect. 3.1. The second half of the first season shows data points rejected from the analysis due to the quality cuts applied. Blue squares are yearly-binned points used for the low-frequency PSD in Sect. 3.2.

In the text

	Fig. 4. Ensemble power spectrum for the master sample of the quasars. Different colours and symbols show PSD estimates calculated using various light curves, as labelled (see Sect. 4 for details).
In the text

	Fig. 5. Top panel: Best-fit PL, DRW, and BPL models to the ensemble PSD of the master sample. Horizontal lines mark the noise floor level for each model. Bottom panels: Respective best-fit residuals plots; i.e. observed PSD–model+error.
In the text

	Fig. 6. Top panel: Ensemble power spectra for the quasars in the BHM1, BHM2, and BHM3 samples (mean values of the BH mass in each sample are reported in the legend). Dashed, solid, and dotted lines indicate the the best-fit BPL models to the power spectra. Thinner lines of the same colour and style are noise-subtracted. Bottom panels: Best-fit residuals.
In the text

	Fig. 7. Same as in Fig. 6, but for the AR1, AR2, and AR3 samples (mean values of the accretion rate for each sample are reported in the legend).
In the text

	Fig. 8. Schematic representation of the best-fit bending power-law PSD model and its scaling with BH mass and accretion rate.
In the text

	Fig. 9. Ensemble power spectrum for the master sample of quasars, as in Fig. 4, with the addition of estimates from the RBF-interpolated light curve combining all the surveys (purple stars).
In the text

Fig. A.1. Colour-magnitude diagrams with PS1 and SDSS magnitudes of non-variable stars. The points are colour-coded by density, with yellow indicating regions of higher point density, and darker areas representing lower-density regions. The upper panel shows the linear correction applied to the g band, while the lower panel demonstrates that the r filters are nearly identical, requiring no correction.

In the text

	Fig. B.1. Ensemble power spectra of simulated DRW light curves with the same observing pattern as the ZTF data. Red circles and squares show the measured PSDs at high and low frequencies, respectively. Solid black line refers to the power spectrum associated to a DRW model with σ = 0.28, and log τ = 3 (left panel), log τ = 3 (middle panel), and log τ = 1 (right panel).
In the text

	Fig. C.1. Example of quasar r-band light curve combining all the surveys (black points), with the RBF interpolated version over a 36-days window (green points).
In the text

Fig. C.2. Ensemble power spectra of simulated DRW light curves with the same observing pattern as combined SDSS+PS1+Gaia+ZTF data, interpolated via RBFs. Red points show the measured power spectral densities, while the black solid line refers to the power spectrum associated to the DRW model. Each panel has DRW light curves with different damping timescales, as listed on the figures, probe regions where the interpolation falls before, within or after the break.

In the text