Ensemble X-ray variability of optically selected QSOs: dependence on black hole mass and Eddington ratio

Although flux variability is one of the defining properties of accretion flows onto supermassive black holes, its dependence on physical parameters such as the mass of the compact object and the Eddington ratio remain under discussion. In this paper we address this issue using the structure function statistic to measure the variability at X-ray wavelengths of a sample of optically selected QSOs with available black hole masses and Eddington ratios. We present a new Bayesian methodology for estimating the structure function tailored to the Poisson nature of the X-ray data. This is applied to 15,548 SDSS DRQ16 QSOs with repeat observations in the XMM-Newton archive and/or the SRG/eROSITA All Sky Survey. The X-ray structure function monotonically increases to time intervals of about 10-15 years, suggesting a preference for scenarios in which instabilities of the accretion disk contribute to the X-ray variability on long timescales. Additionally, there is evidence that the amplitude of the stochastic X-ray flux variations rises with decreasing black hole mass and Eddington ratio. This finding imposes stringent constraints on empirical models of Active Galactic Nuclei variability derived from local samples, emphasizing the significance of high-redshift population studies for comprehending the stochastic flux variations in active black holes.


INTRODUCTION
The accretion of material onto supermassive black holes (SMBHs) at the centres of galaxies is characterised by a high level of stochasticity, which imprints temporal flux variations of the emitted radiation over a wide range of time scales.This fundamental property of the accretion process provides an important observational tool to explore the innermost regions of active SMBHs at spatial scales that are challenging to resolve by other means.This is because flux variations can be translated to light travel times and hence, provide information on the dimensions, spatial structure and dynamics of the accretion flow onto SMBHs (e.g.Uttley et al. 2014;Cackett et al. 2021).This potential has led to systematic multiwavelength studies to better understand the statistical properties and nature of the stochastic flux variations of accretion events onto the SMBHs, i.e.Active Galactic Nuclei (AGN).A central theme in these studies is the dependence of the variability on the fundamental physical parameters of the accreting system, e.g.black hole mass and Eddington ratio.Scaling relations between these physical quantities and the statistical properties of the observed temporal variations would allow the study of accreting systems over ★ E-mail: age@noa.gr a wide range of black hole masses and accretion rates under the same theoretical paradigm.
A popular method to explore the statistical properties of the stochastic AGN flux variability is the Power Spectral Density (PSD), which describes the distribution of the variance of a light curve in Fourier frequencies.X-ray monitoring campaigns of nearby Seyferts over long periods of time reveal PSDs that can be typically described by a broken power-law with slopes changing from about -2 above a break frequency to about -1 below it, similar to X-ray stellar binaries (e.g.Uttley et al. 2002a;Papadakis, I. E. et al. 2002;Papadakis et al. 2003;Markowitz et al. 2003;McHardy et al. 2004;González-Martín & Vaughan 2012).The position of the PSD break frequency is shown to scale with both the mass of the black hole and the Eddington ratio.More importantly, this trend can also be extended to X-ray stellar binaries with black holes many orders of magnitude less massive than AGN (e.g.McHardy et al. 2006).Unfortunately, the direct estimation of the PSD is limited to a small number of nearby bright AGN, for which long and uninterrupted campaigns, mostly from the Rossi X-ray Timing Explorer (RXTE), are available.One approach to overcome this limitation and enable the statistical characterisation of the variability of large samples of AGN is by estimating the normalised excess variance of light curves (Nandra et al. 1997).This quantity is mathematically related to the integral of the PSD over the time-scales sampled by a light curve.It can be measured even in the case of low signal-to-noise observations or not well-sampled light curves (e.g.Allevato et al. 2013).The rich archives of the XMM-Newton and NuSTAR X-ray telescopes has enabled the determination of the normalised excess variance at different timescales for large samples of low redshift AGN (e.g.Ponti et al. 2012;Akylas et al. 2022;Tortosa et al. 2023).These studies reveal correlations of the excess variance, and hence the PSD, on black-hole mass and possibly Eddington ratio, in broad agreement with earlier results based on the PSD analysis of a few local AGN.
Parallel to studies of the light curves of individual objects, there have been efforts to characterise the mean variability properties of high redshift AGN populations detected in extragalactic X-ray survey fields to constrain their PSDs (Paolillo et al. 2004;Papadakis et al. 2008;Paolillo et al. 2017Paolillo et al. , 2023)).This approach hinges on the fact that the total integration time of at least some X-ray surveys has been gradually built up by numerous repeat observations carried out over the course of many years.The light curves of individual sources in such surveys carry limited information but when averaged together they provide useful constraints on the PSD of the population.This approach shows that distant AGN have variability properties that scale with black hole mass and Eddington ratio in a similar fashion as their local counterparts (Georgakakis et al. 2021;Paolillo et al. 2023).However, the exact dependence of the variability on the physical parameters of the accreting systems is under investigation.The Georgakakis et al. (2021) modeling favors a bending power-law PSD with a break frequency that scale with black hole mass (Papadakis 2004) and an amplitude that anti-correlates with Eddington ratio (Ponti et al. 2012).The Paolillo et al. (2023) analysis is consistent with a bent power-law PSD model with constant amplitude and a brake frequency that depends on black hole mass and possibly Eddington ratio McHardy et al. (2006).
In this paper we adopt the ensemble variability approach to further explore the PSD of distant optically selected QSOs and get a handle on scaling relations with black hole mass and Eddington ratio.The variability statistic we chose to use is the structure function (Kozłowski 2016(Kozłowski , 2017)).It provides a measure of the flux variance at a given time scale based on two-epoch photometric measurements for a large sample of sources (e.g. de Vries et al. 2005;MacLeod et al. 2012;Vagnetti et al. 2016;Middei et al. 2017).We combine the SDSS DRQ16 QSO sample with XMM-Newton archival data and new observations from the eROSITA (extended ROentgen Survey with an Imaging Telescope Array, Predehl et al. 2021) instrument onboard the SRG (Spectrum-Roentgen-Gamma) spacecraft to measure the X-ray structure function to timescales of about 10 years at the rest-frame of the sources.Our analysis is similar to previous studies by Vagnetti et al. (2016) and Middei et al. (2017) but our larger sample size enables grouping sources by black-hole mass and Eddington ratio to investigate systematic variations of the structure function with these parameters.Moreover, we also present a new Bayesian algorithm for the estimation the structure function that robustly accounts for X-ray upper limits by modeling the Poisson nature of the X-ray data.Finally the interpretation of our results is supported by the light curve simulation tools presented by Sartori et al. (2018Sartori et al. ( , 2019)).
The analysis presented in this work is similar to the recently accepted paper by Prokhorenko et al. (2024).They also explore the structure function of SDSS QSOs by combining XMM-Newton observations with the eROSITA data in the part of the sky with Galactic longitudes 0 <  < 180 deg, in which the proprietary rights lie with the Russian SRG/eROSITA team.There are however, a number of differences between our work and that of Prokhorenko et al. (2024).Our QSO sample is selected from the latest SDSS DRQ16 catalogue (Lyke et al. 2020, see Section 2) instead of the DRQ14 one (Pâris et al. 2018) used by Prokhorenko et al. (2024).The structure function calculation in our study is based on 2-epoch X-ray observations of SDSS QSOs, in which at least one epoch is from XMM-Newton and the second is either from XMM-Newton (i.e.repeat observations) or eROSITA.Instead Prokhorenko et al. (2024) use only XMM-Newton/eROSITA pairs in their calculations.As a result of the points above the sample used here numbers approximately 15,500 unique QSOS (see Section 2) as opposed to 2344 used by Prokhorenko et al. (2024).We also re-iterate that our analysis uses on a new Bayesian methodology for calculating the structure function (see Section 3) that is adapted to the Poisson nature of the X-ray observations.

THE SAMPLE
The AGN sample used in this work is selected from the Sloan Digital Sky Survey (SDSS) quasar catalogue data release 16 (DR16Q Lyke et al. 2020).It provides photometric and spectroscopic information for 750,414 spectroscopically confirmed QSOs targeted by the SDSS.Black hole masses, bolometric luminosities and Eddington ratios for the DRQ16 QSOs are taken from Wu & Shen (2022).
In this work we explore the X-ray variability properties of the DR16Q QSOs via their structure function (see Section 3).For this calculation X-ray flux estimates for individual sources at two distinct epochs are required.Two independent datasets are used for the determination of the structure function.The first is based on repeat observations in the XMM-Newton archive.The second is combining XMM-Newton data with the German eROSITA All Sky Survey Data Release 1 (eRASS-DE DR1; Merloni et al. 2024).The latter dataset consists of observations carried out in the first six months of the SRG/eROSITA all-sky survey (eRASS1) whose proprietary rights lie with the German eROSITA consortium (eROSITA-DE).The statistical methodology adopted in this work to infer the X-ray structure function of DR16Q QSOs does not rely only on X-ray detected sources either in the XMM-Newton archive or the eRASS-DE DR1.Instead we perform aperture X-ray photometry at the positions of all DR16Q QSOs and then feed this information into the Bayesian model described in Section 3 to determine the structure function of the population at fixed rest-frame timescales.In the case of the XMM-Newton the RapidXMM database (Ruiz et al. 2022) is used to query aperture photometry products in the 0.2-2 keV spectral band (X-ray photons, background level, mean exposure time) at the positions of DR16Q QSOs.The RapidXMM aperture size is fixed to 15 arcsec radius in the case of pointed XMM-Newton observations.This radius roughly corresponds to about 70% of the Encircled Energy Fraction (EEF) of the XMM-Newton PSF (Ruiz et al. 2022).For eRASS-DE DR1 the apetool task of the eROSITA Science Analysis Software System (eSASS, Brunner et al. 2022) is used to extract X-ray photon counts, estimate the background level and corresponding mean exposure time in the 0.2-2.3keV energy band at the positions of DR16Q QSOs that overlap with the eRASS-DE DR1 footprint.The adopted aperture size in this case corresponds to the 75% EEF of the eROSITA PSF at the position of interest.
There are a total of 12,337 unique SDSS QSOs with at least two independent observations in the XMM-Newton archive.This number excludes sources (total of 2116) that lie close to either CCD gaps of the EPIC (European Photon Imaging Camera; Turner et al. 2001) detectors or the edge of the field of view of the XMM-Newton.The identification of such sources is based on the RapidXMM flag-ging scheme (see Table 2 of Ruiz et al. 2022, flagbits 0 and 1).Contamination of the RapidXMM photometric products by photons associated with the wings of the PSF of nearby X-ray sources is an issue that can potentially bias the analysis.This is addressed by first cross-matching the optical positions of the 12,337 SDSS QSOs with the 4XMM-DR12 serendipitous X-ray source catalogue (Webb et al. 2020) within a radius of 6 arcsec.This yield a total of 9474 associations.For the density of the 4XMM-DR12 X-ray detections (630,347 unique sources over an area of 1283 deg 2 ) and the adopted matching radius the expected chance association rate is 0.4%.We then search for additional QSOs in the DRQ16 catalogue that lie close to (within 45 arcsec) 4XMM-DR12 X-ray sources but are not associated with them (i.e. they do not belong to the subsample of the 9474 4XMM-DR13/DRQ16 identifications).The X-ray aperture photometry of such QSOs may be contaminated by photons associated with the PSF wings of the nearby X-ray detections.There are 598 DRQ16 QSOs that lie within 45 arcsec off the positions of 4XMM-DR12 X-ray sources.These are removed from the analysis.This results to a sample of 11,739 SDSS QSOs with a least two independent XMM-Newton observations.
The overlap between eRASS-DE DR1 footprint and the XMM-Newon archive yields a total of 9075 unique DRQ16 QSOs with at least one pair of observations in the two datasets.This number excludes systems (total of 430) located close to CCD gaps of XMM-Newton selected as explained above.We further identify SDSS QSOs that are in the vicinity to X-ray detections in either the XMM-Newton data or the eRASS-DE DR1 catalogue but are not associated with them.In the case of the XMM-Newton, such sources are selected as explained above.In the case of the eROSITA data, we define exclusion regions of 45 arcsec around all eRASS1 X-ray sources with detection likelihood det_ml> 7 and reject any SDSS QSOs that overlap with them.In this exercise we also use the eRASS1 optical identification catalogue presented by Salvato et al. (2024 in prep.) to identify a total of 896 unique DR16 SDSS QSOs that are associated with eRASS-DE DR1 X-ray detections in the Legacy Survey (Dey et al. 2019) Data Release 10 footprint.These sources are allowed in the structure function calculation.The selections above result to 8499 unqiue SDSS QSOs with at least two independent XMM-Newton and eRASS-DE DR1 observations.The total number of unique DRQ16 QSOs with at least two independent XMM-Newton observations or independent XMM-Newton and eRASS-DE DR1 observations is 17,940.X-ray diffuse emission associated with the hot gas of clusters of galaxies may contaminate the aperture photometry at the positions of SDSS QSOs, particularly in the case of the broad PSF of eROSITA (≈ 30 arcsec HEW).We therefore remove SDSS QSOs that lie in the close vicinity of known X-ray selected clusters (see Section 5.2 of Merloni et al. 2024) culled from the ROSAT All-sky Survey (RXGCC, Xu et al. 2022), the XMM Cluster Survey (XCS, Mehrtens et al. 2012), the XMM CLuster Archive Super Survey (X-CLASS, Clerc et al. 2012), the XMM-XXL survey (XXL365, Adami et al. 2018), the eROSITA Final depth Equatorial Survey (eFEDS, Liu et al. 2022).Each of these clusters is assigned an exclusion radius that is empirically determined and is different for the XMM-Newton and eROSITA observations because of the different PSF of the two telescopes.The choice of the size of the exclusion radius is a tradeoff between large QSO sample and low level of contamination from X-ray photons associated with cluster emission.For the eROSITA observations the exclusion radius is defined as 0.5 ×  500 in the case of systems with measured  500 parameter, or 500 kpc in the case of clusters with no  500 estimates, or 1 arcmin for systems with unknown redshifts.For the XMM-Newton observations the exclusion radius is defined as 0.15 ×  500 in the case of systems with measured  500 parameter, or 150 kpc in the case of clusters with no  500 estimates, or 0.5 arcmin for systems with unknown redshifts.This step filters out a total of 1573 unique SDSS QSOs.
SDSS QSOs associated with blazars may show X-ray variability that is not stochastic.We therefore also remove 79 SDSS QSOs that lie within 6 arcsec off blazars taken from the 5th edition of the Roma-BZCAT catalogue (Massaro et al. 2015).Blazars are a subset of the class of radio loud QSOs, the X-ray flux of which is often dominated by emission from the radio jets.As a result the high energy temporal properties of such sources may not be directly related to the X-ray corona and the accretion flow onto the central supermassive black hole.We quantify the radio loudness of a source via the relation (Kellermann et al. 1989) where  5 GHz and  4400 Å are the monochromatic flux densities at 5 GHz and 4400 Å respectively.For the determination of the  4400 Å we first estimate the flux density at rest-frame 2500 Å using the apparent -band magnitude as explained in Section 5 of Richards et al. (2006).The flux density at rest-frame 2500 Å is converted to  4400 Å assuming a power-law spectral energy distribution of the form   ∝  −0.44 .For the estimation of  5 GHz we use two different large-area radio surveys.The FIRST (Faint Images of the Radio Sky at Twenty centimeters Becker et al. 1995) survey and the 2nd data release of LoTSS (LOFAR Two-metre Sky Survey, Shimwell et al. 2017Shimwell et al. , 2022)).The former covers about 10,000 deg 2 and most of the SDSS area to a median 1 rms sensitivity of 0.141 mJy beam −1 at 1.4 GHz.The LoTSS DR2 currently covers about 5,700 deg 2 at 144 MHz to a deeper flux density limit compared to FIRST with a median 1 rms sensitivity of 83 Jy beam −1 .The SDSS DRQ16 catalogue has already been matched to the positions of sources detected in FIRST.For QSOs without FIRST counterparts we estimate 3 upper limits to their 1.4 GHz flux density using the FIRST RMS noise maps 1 .These upper limits are given by the relation 0.25 + 3  , where   is the FIRST RMS noise at the source position and 0.25 is the CLEAN bias correction (White et al. 1997).The association of the SDSS DRQ16 QSOs with the LoTSS catalogues uses the sky positions of the optical counterparts of LoTSS sources presented by Hardcastle et al. (2023).For QSOs without LoTSS matches we adopt an approximate 3 upper limit to their 144 MHz flux densities of 0.3 mJy based on the median RMS noise of 83 Jy beam −1 .For both radio surveys a power-law radio spectral shape with index −0.8(i.e.  ∝  −0.8 ; Sánchez-Sáez et al. 2018) is assumed to convert the observed FIRST 1.4 GHz and/or LoTSS 144 MHz flux densities to rest-frame ones at 5 GHz.
The traditionally adopted cut for defining radio quiet AGN is   < 10 ( Kellermann et al. 1989).Recent observational evidence nevertheless suggests that the contribution of the radio jet emission to the X-ray flux of AGN is subdominant until   > ∼ 100 (e.g Zhu et al. 2020).In this work we therefore adopt the less conservative threshold of   = 30 to differentiate radio loud from radio quiet DRQ16 QSOs (Timlin et al. 2020).Figure 1 compares the radio loudness parameters of DRQ16 QSOs estimated using the FIRST and LoTSS flux densities.This figure demonstrates that the LoTSS radio loudness upper limits provide much stronger constraints on the radio properties of SDSS QSOs compared to the shallower FIRST survey.For example, there are 6083 unique DRQ16 QSOs that have both FIRST and LoTSS coverage.From this sub-sample a total of 3410 have FIRST radio loudness upper limits that lie above our adopted threshold   = 30.For these sources it is not possible to conclude, based on the FIRST data alone, if they are radio quiet.The deeper LoTSS data however, yield   upper limits or measurements that are below the   = 30 threshold for 98% of these sources (3350 out of the 3410).These numbers show that the DRQ16 sub-sample without FIRST associations is dominated by radio quiet QSOs.We therefore define as radio loud those DRQ16 sources that are associated with radio detections in either FIRST or LoTSS and have   > 30 in either of these surveys.These sources are exclude from the analysis.Our sample may contain residual radio loud QSOs associated with FIRST non-detections.Such sources however.represent less than 2% of the total sample and are unlikely to have an impact on our results and conclusions.
It is recognised that the redshift estimates of the SDSS DRQ16 QSO catalogue become uncertain toward fainter optical magnitudes and lower signal-to-noise optical spectra (e.g.Menzel et al. 2016).The DRQ16 catalogue zwarning flag produced by the SDSS spectral reduction pipeline and the sn_median_all parameter, which provides an estimate of the signal-to-noise ratio of individual spectra, can be used to identify potentially problematic sources.A conservative approach proposed by Menzel et al. (2016) to select highly reliable redshifts is to retain only sources with zwarning=0 and sn_median_all> 1.6.This approach removes about 30% of the sample, inlcuding a large fraction of sources with bona-fide and reliable QSO redshifts (based on visual inspection of a subsample).In this work we therefore choose not to apply a redshift quality selection.Instead, we repeat the analysis presented in the following section for the sample with zwarning=0 and sn_median_all> 1.6 and confirm that the main trends and results of Section 4 are recovered for that conservative sample.
Following the selection criteria described above, the final (baseline) sample used in this paper to determine the X-ray structure function of SDSS DRQ16 QSOs consists of a total of 15,548 unique QSOs and 110,829 pairs of X-ray observations at distinct epochs.Figure 2 shows the distribution of the baseline subsample on the Eddington ratio vs black hole mass plane.In Figure 3 we present the distribution of the sample on the bolometric luminosity vs redshift plane.Figure 4 shows the histogram of the rest-frame time differences, Δ, of the sample.This quantity is defined as the time span between the two epoch X-ray observations of a given source divided by the factor (1 + ), where  is the redshift of the source.Figure 4 shows that the sample is probing timescales that extend to about 1 decade at rest-frame.Also, the contribution of the XMM/eROSITA flux measurement pairs increases toward longer timescales.In our analysis we use rest-frame timescales above the limit 10 −2 yr (i.e.few days).For shorter time intervals the expected flux differences between the two-epoch observations are small and systematics (e.g.uncertainties in the estimation of the encircled energy fraction for the RapidXMM counts) become important.In any case Figure 4 shows that the statistical power of the sample is for time intervals ranging from months to several years.These are the timescales of interest to this work.

STRUCTURE FUNCTION ESTIMATION
Multiple definitions of the structure function (SF) exist in the literature depending on whether the data are in flux or apparent magnitude units (Kozłowski 2017).In this work we define the structure function at a given time scale  as where   (),   ( + ) are the X-ray fluxes of a source measured at epochs separated by , i.e. at times  and  + .For a population of AGN the equation above estimates the variance of the logarithmic flux differences at the timescale .In practice observations are noisy and therefore any contribution from random flux errors must be accounted for in the SF estimation.This is usually done by modifying Equation 2 as (e.g.see Middei et al. 2017;Kozłowski 2017) where  2  () represents the uncertainty associated with the flux measurement log   ().X-ray observations are typically described by the Poisson distribution and therefore the shot noise term above depends on the instantaneous source flux and the instrumental background level.In the case of bright sources the errors may be approximated to an acceptable level of accuracy by the normal distribution.In the case of low photon counts however, such as the eRASS1 observations or shallow XMM-Newton data, this approximation breaks down.For example, many of the SDSS DRQ16 QSOs are not formally detected in the eRASS1 (see Section 2) resulting in upper flux limits, which may nevertheless provide useful constraints on the SF.For the calculation of the SF we therefore develop a custom Bayesian methodology that properly accounts for the Poisson nature of the data and allows the  consistent treatment of both formal X-ray detections and X-ray upper limits.At the core of the methodology is aperture photometry as described below.
Suppose a source  with X-ray observations in two distinct epochs.In the first dataset the source has  ,1 counts within an aperture with EEF    ,1 .The background level is assumed to be  ,1 and the exposure time at the source position is  ,1 .The probability that the source has flux  ,1 can be expressed by the Poisson probability Where  1 signifies the observations at hand (e.g. ,1 ,  ,1 etc) and the Poisson expectation value  ,1 is ,1 is the energy to flux conversion factor and depends on the source's spectral model and the characteristics of the X-ray telescope/detector that observed it.In the second X-ray observation carried out at a later epoch by the same or different telescope the source  has  ,2 counts within an aperture of EEF    ,2 , the background level is  ,2 and the exposure time is  ,2 .The probability that the source's flux has changed by Δ log   = Δ log(  ,1 /  ,2 ) = log  ,1 − log  ,2 to  ,2 is also given by the Poisson formula Here  2 signifies the second set of observations at hand.The flux variation is expressed by the ratio  ,1 /  ,2 , which links to the definition of the structure function in Equation 2. The Poisson expectation value  ,2 is Where  ,1 is the energy to flux conversion factor for the second observation.The likelihood of the two epoch observations of source  is the product of the Poisson probabilities of Equations 4, 6 and the probability of a flux ratio variation  ,2 /  ,1 .For the latter, it is assumed that the variations of the logarithmic flux ratio, log  ,2 /  ,1 , are drawn from a normal distribution with mean of zero and variance that depends on the (rest-frame) time interval Δ between the two observations.As a consequence of Equation 2, the variance of the adopted normal distribution is by definition the structure function at timescale Δ, which is to be constrained from the observations.Combining all the components above it is possible to write down the likelihood for the 2-epoch observations of source  as The last term describes the normal distribution with parameters  (mean) and  (scatter).In the case of  number of sources with two epoch observations separated by the same time interval Δ the likelihood of the population is The likelihood above is used to perform Bayesian inference on , which represents the structure function of the sample at a fixed time interval.The Hamiltonian Markov Chain Monte Carlo code Stan 2 is used to sample the likelihood of Equation 9 and produce parameter posterior distributions.The free parameters of the model are the  of the normal distribution in Equation 9 and the first epoch logarithmic flux log  ,1 of individual pairs.During the sampling process we adopt a flat prior on  between 10 −5 and 2. The log  ,1 is assigned a broad normal prior with mean −14 (i.e.flux of 10 −14 erg s −1 cm −2 ) and (logarithmic) scatter parameter of 10.Appendix A tests and validates the Bayesian SF estimation method described above using simulated light curves (Sartori et al. 2018(Sartori et al. , 2019) ) and X-ray photometric products that mimic those of the XMM-Newton and eRASS1-DE DR1 observations of SDSS DRQ16 QSOs.

The X-ray structure function of the full sample
We first estimate the structure function of the full sample of DRQ16 QSOs described in Section 2. The results are plotted as a function restframe time scale, Δ, in Figure 5.The SF monotonically increases with increasing ΔT in agreement with previous studies (e.g.Vagnetti et al. 2016;Middei et al. 2017).Assuming that the structure function is described by a power-law of the form  ∝ Δ  , we estimate an 2 https://mc-stan.orgindex  = 0.12 ± 0.01.This value translates to a Power Spectrum Density (PSD) power law index of about -1.2 (e.g.Bauer et al. 2009).
It is interesting to explore how these results compare with the predictions of empirical models for the PSD of AGN derived from detailed monitoring campaigns of local Seyferts.Because there is no analytic relation between PSD and SF, the generation of model predictions requires the intermediate step of generating simulated light curves for a given input PSD.The corresponding structure function at different rest-frame timescales can then be estimated numerically from these light curves.We use the method described in Emmanoulopoulos et al. (2013) as implemented in Sartori et al. (2018Sartori et al. ( , 2019) ) to generate light curve simulations.In this exercise we adopt a PSD model that follows a bending power-law form (McHardy et al. 2004;González-Martín & Vaughan 2012) where  is the frequency.The amplitude, , and break frequency,   , are suggested to depend on the physical parameters of the accreting system, e.g. the mass of the black hole and its Eddington ratio.Different observational studies of local AGN have proposed over the years different parameterisations for  and   .In our analysis we use three popular models put forward in the literature.The first one, which we refer to as Model 1, assumes a constant amplitude  = 2 •   •  (  ) = 0.02 and break frequency,   , that depends on the black hole mass of the accreting system as This parametrisation is motivated by the observations of Papadakis (2004) and González-Martín & Vaughan (2012).Simulated light curves are generated (Emmanoulopoulos et al. 2013;Sartori et al. 2018Sartori et al. , 2019) ) for two black hole masses log   / ⊙ = 8, 9 that bracket the mode of the distribution shown in Figure 2. The structure function that is numerically determined from these light curves is overplotted in Figure 5.It is clear that the PSD Model 1 predictions are inconsistent with the measured structure function of SDSS DR16 QSOs.Modifications to this model are therefore required to account for this inconsistency.One possibility is to consider different dependences of the PSD break frequency and and amplitude on black hole mass and Eddington ratio.We test two variants of Model 1.The first adopts the constant PSD amplitude of Model 1 and further assumes that the break frequency is a function of both the black-hole mass and the accretion rate as proposed by McHardy et al. (2006).This dependence is expressed in terms of the AGN bolometric luminosity ( bol ) as We refer to this PSD parametrisation as Model 2. We also consider a parameterisation in which the assumption of a constant PSD amplitude is relaxed.This is motivated by observations of local Seyferts suggesting that the amplitude of the PSD scales with Eddington ratio,   , as (Ponti et al. 2012) In this particular parameterisation, which is referred to as Model 3, the break frequency depends on black hole mass as in Equation 11.
Figure 5 shows the structure function of Models 2 and 3 for black hole masses log   / ⊙ = 8, 9 and Eddington ratio   = 0.1 that approximately corresponds to the mode of the distribution in Figure 2.For this choice of physical parameters the Model 2, 3 structure function predictions are in better agreement with the observations, particularly for log   / ⊙ = 8.We emphasise that the model curves plotted in Figure 5 are not fits to the structure function data points.Instead they help visualise the expectations of empirical models derived from extensive observational studies of nearby active black holes.In any case, the ensemble X-ray variability properties of optically selected QSOs at cosmological distances provide strong evidence that the PSD of AGN depends on both black hole mass and Eddington ratio.It is therefore interesting to further explore such a dependence by directly measuring the structure function of DRQ16 QSOs grouped by black hole mass and Eddington ratio.

Black hole mass dependence of the X-ray structure function
We select DRQ16 QSOs in a narrow Eddington ratio interval, −1.25 < log   < −0.5 that brackets the mode of the distribution in Figure 2.This sub-sample allows us to explore variations of the structure function with black hole mass at nearly fixed Eddington ratio.The choice of the Eddington ratio interval above is a trade-off between sufficiently large sample size and the desire to minimize the impact of Eddington ratio on the measured structure function.
We further split this sub-sample into two nearly equal size black hole mass bins at the threshold log   / ⊙ = 8.75.The structure function is estimated independently for the two sub-samples and the results are shown in Figure 6.DRQ16 QSOs with lower black-hole masses have a systematically higher structure function compared to more massive ones.This difference is not sensitive to the choice of black hole mass bins.Splitting for example the sample into three nearly equal size groups of black hole mass yields similar results, i.e. the structure function is systematically lower for QSOs with the most massive black holes (see Fig. C1 and discussion in Appendix C).Moreover, the observed trend is not dominated by a few high signal-to-noise ratio sources in the sample but is instead representative of the full population.We test this by re-running the analysis after removing observations in which individual QSOs have more than 100 X-ray counts (36% of pairs removed) or 250 X-ray counts (12% of the pairs removed).It is found that the main trend shown in Figure 6 holds for these smaller sub-samples, i.e. a systematically higher SF for lower black hole masses.
Although the structure functions of both sub-samples plotted in Figure 6 increase on average with increasing time-scale, there is scatter among neighboring data points that is larger than their errorbars.Moreover, there is one time-scale bin in Figure 6 for which the general trend of the SF with   is inverted, i.e. the high black hole mass sub-sample is more variable.Uncertainties in the determination of black-hole masses (≈ 0.5 dex, Wu & Shen 2022) may weaken or even invert any intrinsic correlation between structure function amplitude and black hole mass.
Additionally, although we attempt to control for variations of parameters that are important for the determination of the SF, e.g.black-hole mass, Eddington ratio or rest-frame time-scale, by selecting QSOs in narrow parameter slices, the distribution of these parameters within individual time-scale and black-hole bins may be different, thereby leading to increased scatter.We demonstrate this point in Figure 7 that shows for each time-scale bin of the data points in Figure 6 the distribution of key parameters, such as Δ,   ,   and redshift for the two QSO sub-samples selected by black hole mass.It is clear from this figure that e.g. the mean   of the two sub-samples is similar across different time-scale bins, but the underlying distributions span a range of values and in some cases show subtle differences.This is further demonstrated by applying the K-S (Kolmogorov-Smirnov) test to compare the histograms shown in each panel of Figure 7 and estimate the probability of the null hypothesis that they are drawn from the same parent population.For certain   panels the null hypothesis can be rejected with a probability of about 3 × 10 −3 , i.e at approximately the 3 confidence level.We caution against false positives in the case of hypothesis testing with multiple comparisons and the need to consider corrections to the estimated -values to account for this effect (e.g.Bonferroni adjustment).
The K-S test shows that the lowest null hypothesis probabilities occur for the redshift panels of Figure 7, with the high black-hole mass sub-samples being typically skewed to higher redshifts as a result of observational selection effects (e.g.flux limit, larger cosmological volumes).Variations of the mean redshift of a given black-hole mass subsample across different rest-frame time-scale bins are also visible in Figure 7.These effects translate to different rest-frame energies at which the variability of individual sources is estimated.Although the dependence of the variability amplitude on rest-frame energy is second order (e.g.Vagnetti et al. 2016;Ponti et al. 2012), it may contribute to the scatter among different data points in Figure 6.
Also plotted in Figure 6 are the structure function predictions of PSD Models 2 and 3 for different black holes masses and Eddington ratios that correspond to the averages of the two observational subsamples.We re-iterate that these curves are not fits to the data but instead serve to guide the reader on the expectations based on our understanding of the variability of AGN from monitoring campaigns of local samples.Qualitatively these models are consistent with the observations, i.e. the corresponding structure functions at fixed timescale decrease with increasing black hole mass.Quantitatively however, the model structure functions are systematically lower than the observations.This suggests that small modifications may be necessary to the scaling factors that govern the relation between PSD amplitude and/or break frequency with black hole mass or Eddington ratio, i.e.Equations 11, 12, 13.

Eddington ratio dependence of the X-ray structure function
Next, we select DRQ16 QSOs in a narrow black-hole mass interval, 8.5 < log   / ⊙ < 9.0 that brackets the mode of the distribution in Figure 2.This sub-sample is used to explore variations of the structure function with Eddington ratio at nearly fixed black hole mass.We define two Eddington ratio bins split at log   = −1.At this cut the two sub-samples have nearly equal number of 2-epoch X-ray observation pairs.The structure function is estimated for the two sub-samples and the results are plotted in Figure 8.There are four time scale bins where the structure function of low Eddington ratio DRQ16 QSOs is higher compared to that of high Eddington ratio sources.For three time scale bins the structure functions of the two samples are consistent within the uncertainties.Only at the highest time-scales probed by the current sample (≈ 10 yrs) is the structure function of the high Eddington ratio sources larger than that of low Eddington ratio systems.Although the differences between the two sub-samples in Figure 8 are not as clear as in the case of black-holes mass bins (i.e. Figure 6), there is evidence that low Eddington ratio DRQ16 QSOs tend to be more variable than high Eddington ratio systems, at least for time-scales between about 0.1-5 yrs.Splitting the sample into three nearly equal size groups of Eddington ratio yields qualitatively similar conclusions, i.e. the structure function tends to be higher for DRQ16 QSOs with low Eddington ratios (see  Middei et al. (2017) for SDSS DRQ12 QSOs (Pâris et al. 2017) using XMM-Newton and ROSAT observations.Also shown are the structure function predictions of the PSD models 1, 2 and 3 parametrised by Equations 10-13.Fig. C2 and discussion in Appendix C).The same trend is also recovered after removing sources with a high number of counts, e.g.> 100 (40% of pairs removed).We reiterate that uncertainties in the estimation of black hole masses propagate into the Eddington ratio determination and could wash out intrinsic correlations between amplitude of variability and Eddington ratio.
Figure 9 plots for each time-scale bin of Figure 8 and for each of the two QSO sub-samples selected by Eddington ratio the corresponding histograms of the quantities Δ,   ,   and redshift.We observe similar trends with black-hole mass, redshft and rest-frame time-scale as those already discussed in the case of Figure 7. Therefore part of the scatter observed in Figure 6 among neighboring data points is related to differences in the distribution of the parameters above for the various sub-samples.Also plotted in Figure 8 are the structure function predictions of PSD Models 2 and 3 for the mean Eddington ratios and black hole masses of the two observational subsamples.In Model 2 the structure function decreases with decreasing Eddington ratio, which is at odds with the observational results in Figure 8. Model 3 instead, is qualitatively consistent with the observations predicting a higher structure function amplitude for lower Eddington ratios.Quantitatively, although this model is in reasonable agreement with the structure function measurements of the low Eddington-ratio sub-sample (mean log   = −1.3), it underestimates the structure function of the high Eddington-ratio sub-sample (mean log   = −0.7).This suggests that modifications need to be applied to Model 3 to make it quantitatively consistent with the measurements in Figure 8. Exploring such changes is beyond the scope of this paper.

SF ( T)
Model 2: M BH = 10 8.5 M ; log Edd = 0.9 Model 2: M BH = 10 9.0 M ; log Edd = 0.9 Model 3: M BH = 10 8.5 M ; log Edd = 0.9 Model 3: M BH = 10 9.0 M ; log Edd = 0.9 logM BH = 8.0 8.75 logM BH = 8.75 9.5 Figure 6.Structure function versus rest-frame time-scale for the SDSS DRQ16 QSO subsample with Eddington ratios in the interval −1.25 < log   < −0.5, which is further split by black hole mass.The orange circles correspond to SDSS DRQ16 QSOs with log   / ⊙ > 8.75.Blue circles are for QSOs with black holes lower than the limit above, i.e. log   / ⊙ < 8.75.The horizontal error bars show the extent of the restframe time interval.Data points are plotted at the mean Δ  of a given sub-sample.Also shown are the SF predictions for the PSD Models 2 (solid lines; Eq. 10, 12) and 3 (dashed lines; Eq. 10, 11, 13).The blue curves correspond to a black hole mass log   / ⊙ = 8.5, while the orange curves assume log   / ⊙ = 9.The adopted Eddington ratio for all curves is log   = −0.9.

DISCUSSION
We present a new Bayesian methodology for measuring the X-ray structure function of large AGN samples and apply it to SDSS DRQ16 QSOs with repeat observations from XMM-Newton and/or eROSITA.This dataset enables measurement of the X-ray structure function to decadal timescales as well as investigations on the dependence of the AGN long-term flux variability on black hole mass and Eddington ratio.
An interesting result from our analysis is that we do not find evidence for a turnover in the structure function of SDSS DR1Q16 QSOs out to the longest time-scales (≈ 10 yr rest-frame) probed by the current sample.Instead Figure 5 shows a monotonic increase of the structure function toward longer time intervals with a (logarithmic) slope of about 0.12.This is consistent with the results of Middei et al. (2017) on the X-ray structure function of SDSS DRQ12 QSOs based on combined ROSAT and XMM-Newton observations (see also Vagnetti et al. 2016).Studies of the variability properties of QSOs in the UV/optical part of the electromagnetic spectrum also suggest a monotonically increasing structure function to decadal rest-frame time intervals with logarithmic slopes similar to those estimated by our analysis (de Vries et al. 2005, but see MacLeod et al. 2012 for possible evidence for flattening).These findings have implications on the size of accretion disk and the origin of the long-term variability at X-ray wavelengths.In particular, there is an increasing body of evidence suggesting that the UV/optical variability of  Orange circles are for QSOs with Eddington ratio larger than the limit above, i.e. log   > −1.The horizontal error bars show the extent of the restframe time interval.Data points are plotted at the mean Δ  of a given sub-sample.Also shown are the SF predictions for the PSD Models 2 (solid lines; Eq. 10, 12) and 3 (dashed lines;Eq. 10,11,13).The blue curves correspond to an Eddington ratio log   = −1.3,while the orange curves assume log   = −0.7.The adopted black hole mass for all the curves is log / ⊙ = 9.
QSOs has two components that dominate at different temporal scales (e.g.Czerny 2006;Breedt et al. 2009;Hernández Santisteban et al. 2020;Yao et al. 2023;Secunda et al. 2023).At short time intervals the reprocessing of the stochastic X-ray flux variations associated with changes in the X-ray emitting corona are believed to drive the UV/optical variability (e.g.Kammoun et al. 2023).On longer timescales magneto-rotational turbulence (Balbus & Hawley 1991) or turbulence generated by convection as a result of opacity variations of the disk material (Jiang & Blaes 2020) introduce accretion disk instabilities that are thought to propagate inwards on the viscous timescale thereby introducing an additional long-term modulation of the light curves of AGN.Such long term variations may also affect the observed X-ray flux variability (e.g.Papadakis et al. 2001;Arévalo & Uttley 2006;Yao et al. 2023;Secunda et al. 2023) and lead to the monotonically increasing structure function of Figure 5.
Next we turn to the dependence of the stochastic flux variations of AGN on the physical parameters of the accreting system such as black hole mass and Eddington ratio.Such dependences provide a handle on the physical processes that drive the observed variability (e.g.Tang et al. 2023) and could enable the unification of active black hole systems with diverse masses and Eddington ratios under the same theoretical paradigm (e.g.McHardy et al. 2006).Figures 6 and 8 indicate that the amplitude of the X-ray structure function of SDSS DRQ16 QSOs decreases with increasing black hole mass and toward higher Eddington ratios.This is consistent with the recent X-ray structure function calculations presented by Prokhorenko et al. (2024) based on SDSS DRQ14 QSOs (Pâris et al. 2018) with XMM-Newton/eROSITA overlapping observations in the part of the sky with proprietary rights belonging to the Russian SRG/eROSITA team.Our results are also consistent with results on the dependence of the UV/optical variability of SDSS QSOs on the physical properties of the accreting system (e.g.Arévalo et al. 2023).Similar trends, but on order of magnitude shorter timescales (typically fraction of a day), are also recovered for the X-ray variability of local Seyferts based on either the normalised excess variance (Ponti et al. 2012; -2.0 -1.8 -1.6 6.6e-01 6.6e-01 [ 2.0, 1.6] 8.5 9.0 6.8e-01 6.8e-01 -2 -1 0 0 1 2 3 4   et al. 2004, 2006).
It is also interesting to explore different empirical scaling relations that can potentially describe the observed dependence of the structure function on physical parameters.The fundamental quantity in this exercise is the PSD of the observed variations that describes the variability amplitude as a function of Fourier frequency.We first assume that the PSD is described by the bending power-law of Eq. 10 with parameters (i.e.break frequency, amplitude) that scale with black hole mass and/or Eddington ratio.We then consider different parametrisations for these relations derived from detailed studies of the temporal properties of local Seyeferts.It is shown that PSD models with constant amplitude and a break frequency that is only a function of black hole mass (see Eq. 11, Papadakis, I. E. et al. 2002;Uttley et al. 2002b) are inconsistent with the X-ray variability of distant QSOs (see Fig. 5).Moreover, PSD models in which   scales with both black hole mass and Eddington ratio (see Eq. 12, McHardy et al. 2004McHardy et al. , 2006) ) predict a dependence of the structure function on Eddington ratio that is at odds with the observations in Fig. 8.The PSD model that better describes the X-ray temporal properties of DRQ16 QSOs, at least qualitatively, is the one proposed by Ponti et al. (2012) with an amplitude that is inversely proportional to Eddington ratio (Eq.13) and a break frequency that decreases with increasing black hole mass as in Eq. 11.
Similar conclusions are reached from studies of the ensemble variability of a different sample of high redshift active supermassive black holes, X-ray selected AGN in the Chandra Deep Field South (Luo et al. 2017).The mean normalised excess variance of these sources is found to decrease with increasing X-ray luminosity (Paolillo et al. 2017).Modelling suggests that this dependence cannot be reproduced by a PSD model with a constant amplitude and a break frequency given by either Equation 11 or Equation 12 (Georgakakis et al. 2021).Instead, a PSD amplitude that scales with Eddington ratio as in Equation 13 provides a better description of the ensemble variability of X-ray selected AGN in the Chandra Deep Field South (Georgakakis et al. 2021).Recently, Paolillo et al. ( 2023) combined archival observations of local Seyferts and higher redshift X-ray selected broad optical emission-line AGN detected in deep extragalactic surveys (COSMOS, Chandra Deep Field South) to measure the ensemble normalised excess variance of the population at different time scales and as a function of black hole mass.They adopt a PSD model with a constant amplitude (i.e.independent of redshift or Eddington ratio) and a break frequency that scales with black hole mass to demonstrate that it provides a satisfactory description of their variability data.Nevertheless, Paolillo et al. (2023) do not explore the consistency of their observations with more complex PSD model parameterisations that e.g.include a dependency of the amplitude on Eddington ratio.Further work is therefore needed to compare our results and conclusions with those of Paolillo et al. (2023).At any rate, our analysis demonstrates the power of ensemble X-ray flux variability studies of high redshift AGN samples as a means of constraining empirical PSD models.
Dedicated spectroscopic follow-up programmes of the eROSITA X-ray sources, e.g.SDSS-V Black Hole Mapper (Kollmeier et al. 2017;Almeida et al. 2023), 4MOST AGN Survey (Merloni 2019), in combination with the multi-epoch X-ray photometry from the successive eROSITA All Sky Surveys (eRASS1-5), will significantly increase the sample of AGN for ensemble variability investigations.These datasets can help refine and improve the observational trends presented in this work by both increasing the statistics and providing a wider black-hole mass and Eddington ratio baselines.

CONCLUSIONS
This paper revisits the X-ray structure function of SDSS QSO samples using a novel Bayesian methodology and combining XMM-Newton archival data with new eROSITA observations carried out in the first six months of the SRG/eROSITA all-sky survey (eRASS1) and whose proprietary rights lie with the German eROSITA consortium.
We find no evidence for a turnover in the X-ray structure function to the longest time-scales of about 10 yr rest-frame probed by our sample.This finding may have implications on the size of the accretion disk, under the assumption that the long-term variations of the accretion flow drive the variability at X-ray wavelengths on decadal timescales.
There is also evidence that the structure function of DRQ16 QSOs on timescales of few months to years increases with decreasing black hole mass and toward lower Eddington ratios.This dependence is best represented by PSD models with an amplitude that scales with Eddington ratio.
Our analysis highlights the importance of ensemble variability studies as a complement to the characterisation of the temporal properties of individual AGN.

APPENDIX A: VALIDATION OF THE BAYESIAN STRUCTURE FUNCTION ESTIMATION METHOD
We use simulated light curves with known statistical properties, i.e. power spectrum shape and normalisation, to validate the Bayesian approach of Section 3 for the estimation of the structure function of QSO populations.
The starting point of the validation exercise is a light curve generated by the code presented by Sartori et al. (2018Sartori et al. ( , 2019) ) based on the methods of Emmanoulopoulos et al. (2013).For a given rest-frame time-scale, Δ, we can then randomly choose  sample number of data point pairs on the light curve each of which are separated by the time interval Δ.These pairs represent the 2-epoch X-ray photometry of a population of QSOs with size  sample that can be used to determine their structure function.
Before drawing  sample pairs however, mean fluxes need to be assigned to the light curves that are representative of the SDSS DRQ16 QSO sample.Using the RapidXMM aperture photometry of the SDSS QSO sample we estimate 0.2-2 keV fluxes for individual sources following the approach explained in Appendix B and then build the distribution of the sample shown in Figure B1.We use this historgam to randomly draw  sample fluxes that are representative of the SDSS DRQ16 QSO sample.These are used to normalise the mean of the simulated light curves.The end product of this step are  sample light curves with the same power spectrum but different mean fluxes.It is these  sample renormalised light curves that we use to draw (from each) two random fluxes separated by the time interval Δ.These flux pairs are then converted to XMM-Newton EPIC (first epoch) and eROSITA (2nd epoch) photon counts or just XMM-Newton photon counts (both epochs).This conversion uses exposure times and backgrounds (XMM-Newton or eROSITA) randomly drawn from the RapidXMM or apetool X-ray aperture photometry products extracted at the positions of SDSS QSO as described in Section 2. This process yields simulated data sets that are similar to real observations that can be analysed using the Bayesian methodology of section 3 to infer the structure function of the population at the time-scale Δ.The inferred SF can then be compared with the known structure function of the input light curve.
Figure A1 shows the results of this experiment by comparing the inferred SF with the input one as a function of the sample size,  sample .Two separate sets of data points are shown in this figure.The first combines XMM-Newton EPIC (first epoch) with eROSITA-DE DR1 (2nd epoch) photon counts to infer the SF.The second assumes XMM-Newton photon counts for both epochs to derive the structure function.Moreover, different timescales Δ and/or PSD normalisations are used as input, which result in three different values of the input structure function that are to be recovered.These are shown with the three horizontal lines in Figure A1.This figure demonstrates that when increasing  sample the inferred SF converges to the input one.Additionally, the inferred SF is always consistent with the input one within the 1 uncertainty estimated from the posteriors returned by the Bayesian methodology of Section 3.These results test and validate the Bayesian structure function calculation developed in Section 3.

APPENDIX B: X-RAY FLUX ESTIMATION USING APERTURE PHOTOMETRY PRODUCTS
This Section describes the determination of X-ray fluxes for the SDSS QSOs based on the RapidXMM aperture photometry.Suppose a source on the detector  of a given XMM-Newton observation for which the RapidXMM database lists   total number of photons, a background level   and an exposure time   .The probability that this source has flux   is given by the Poisson probability distribution function where  = 3 =1   is the sum of the source's photons in each of the three EPIC detectors.The Poisson expectation value in the equation above can be written as where   is the count to flux conversion factor for the camera  and     is the encircled energy fraction of the RapidXMM aperture for the detector  (see Table 1 of Ruiz et al. 2022).The summation is over the three EPIC cameras, PN, MOS1 and MOS2.The count to flux conversion factor is estimated for each EPIC camera separately assuming a power-law X-ray spectrum with Γ = 1.9 absorbed by the appropriate Galactic hydrogen column density (Kalberla et al. 2005) in the direction of the SDSS QSO under consideration.Equation B1 can be solved numerically to determine the X-ray flux probability

Figure 1 .
Figure 1.Comparison of the radio loudness parameter (see Section 2) estimated using either the FIRST data (y-axis) or the LOFAR LoTSS survey (x-axis).All the data points are SDSS DR16 QSOs that lie in the overlap region of two radio samples.The black diagonal solid line shows the one-toone relation.The horizontal and vertical dashed lines show the radio loudness threshold   = 30 adopted in our work to separate radio loud from radio quiet sources.Filled blue circles represent QSOs detected in both radio surveys.The orange open triangles pointing downward are SDSS QSOs that are detected in LoTSS but not in FIRST.The green open triangles pointing to the left correspond to SDSS QSOs that are detected in FIRST but not in LoTSS.The pink crosses are SDSS QSOs with upper limit radio flux density measurements in both surveys.They are offset from the one-to-one relation because of the different depths of the FIRST and LoTSS radio surveys.

Figure 2 .Figure 3 .
Figure 2. Distribution of the baseline sample of SDSS DRQ16 QSOs (see Section 2) on the Eddington ratio vs black hole mass plane.The contour levels are chosen to enclose 34, 68, 95 and 99 per cent of the QSO population.The 1-dimensional projections of this distribution along the black hole mass and Eddington ratio axes are also shown by the histograms on the top and to the right of the main panel respectively.

Figure 4 .
Figure 4. Histogram of the rest-frame time scales probed by the baseline sample of SDSS DRQ16 QSOs (see Section 2).These timescales are estimated as the time difference between either XMM-Newton repeat observations or XMM-Newton and the eRASS-DE DR1 data at the rest-frame of each QSO.The blue hatched histogram corresponds to the full sample.The orange histogram is for the subsample that combines XMM-Newton observations with the eRASS-DE DR1 data.The contribution of the latter sample to the total increases with increasing timescale.

Figure 5 .
Figure 5. Structure function versus rest-frame time-scale.The blue circles are the results for SDSS DRQ16 QSOs based on the combined XMM-Newton and eRASS1-DE DR1 observations.The horizontal error bars show the extent of the rest-frame time interval.Data points are plotted at the mean Δ  of a given sub-sample.The observations are compared against the SF estimates presented by Middei et al. (2017) for SDSS DRQ12 QSOs (Pâris et al. 2017) using XMM-Newton and ROSAT observations.Also shown are the structure function predictions of the PSD models 1, 2 and 3 parametrised by Equations 10-13.

Figure 7 .
Figure7.The histograms in each panel show the distributions of each of the key parameters, log Δ (top row), log   (2nd row from top), log   (3rd row from top) and redshift (last row), for the two subsamples with Eddington ratio in the range −1.25 <   < −0.5 and black holes masses in the intervals log   / ⊙ = 8 − 8.75 (blue lines) and log   / ⊙ = 8.75 − 9.5 (orange lines).Each column of panels corresponds to one of the rest-frame times scales of the data points plotted in Fig.6as indicated at the top of the column.The vertical lines in each panel show the mean of the distribution of the corresponding parameter.The number in each panel corresponds to the Kolmogorow-Smirnow test p-value, i.e. the probability of the null hypothesis that the blue and orange histograms in each panel are drawn from the same parent distribution.

Figure 8 .
Figure 8. Structure function versus rest-frame time-scale for the subsample with 8.5 < log   / ⊙ < 9.0, which is further split by Eddington ratio.The blue circles correspond to SDSS DRQ16 QSOs with log   < −1.Orange circles are for QSOs with Eddington ratio larger than the limit above, i.e. log   > −1.The horizontal error bars show the extent of the restframe time interval.Data points are plotted at the mean Δ  of a given sub-sample.Also shown are the SF predictions for the PSD Models 2 (solid lines; Eq. 10, 12) and 3 (dashed lines; Eq. 10, 11, 13).The blue curves correspond to an Eddington ratio log   = −1.3,while the orange curves assume log   = −0.7.The adopted black hole mass for all the curves is log / ⊙ = 9.

Figure 9 .
Figure9.The histograms in each panel show the distributions of each of the key parameters, log Δ (top row), log   (2nd row from top), log   (3rd row from top) and redshift (last row), for the two subsamples with black holes masses in the range log   / ⊙ = 8.5 − 9.0 and Eddington ratios in the intervals −2.0 < log   < −1.0 (blue lines) and −1.0 < log   < 0.0 (orange lines).Each column of panels corresponds to one of the rest-frame times scales of the data points plotted in Fig.8as indicated at the top of the column.The vertical lines in each panel show the mean of the distribution of the corresponding parameter.The number in each panel corresponds to the Kolmogorow-Smirnow test p-value, i.e. the probability of the null hypothesis that the blue and orange histograms in each panel are drawn from the same parent distribution.

Figure C1 .
Figure C1.Structure function versus rest-frame time-scale for the SDSS QSO subsample with Eddington ratios in the range −1.25 < log   < −0.5 and black hole masses in the intervals log   / ⊙ = 8.0 − 8.5 (blue data point) and log   / ⊙ = 9.0 − 9.5 (orange data points).The vertical error bars correspond to the 68% confidence interval uncertainties in the determination of the structure function.The horizontal error bars show the extent of the rest-frame time interval.Data points are plotted at the mean Δ  of a given sub-sample.

Figure C2 .
Figure C2.Structure function versus rest-frame time-scale for the SDSS QSO subsample with black hole masses in the range 8.5 < log   / ⊙ < 9.0 and Eddington ratios in the intervals −2 < log   < −1.5 (blue data points) and −1 < log   < 0 (orange data points).The vertical error bars correspond to the 68% confidence interval uncertainties in the determination of the structure function.The horizontal error bars show the extent of the rest-frame time interval.Data points are plotted at the mean Δ  of a given sub-sample.