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

A. Georgakakis,1, J. Buchner2,3, A. Ruiz1, T. Boller2, A. Akylas1, M.
Paolillo4,5,6, M. Salvato2,A. Merloni2, K. Nandra2, T. Dwelly2
1Institute for Astronomy & Astrophysics, National Observatory of Athens, V.
Paulou & I. Metaxa, 11532, Greece
2 Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, 85741
Garching, Germany
3Excellence Cluster Universe, Boltzmannstr. 2, D-85748, Garching, Germany
4Dipartimento di Fisica "Ettore Pancini", Università Federico II, via Cinthia,
80126 ,Napoli, Italy
5INAF - Osservatorio Astronomico di Capodimonte, via Moiariello 16, 80131,
Napoli, Italy
6Istituto Nazionale di Fisica Nucleare, Sezione di Napoli, I-80126 Napoli,
Italy E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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 yrs, 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.

###### keywords:

X-rays: general – galaxies: active – Galaxies: nuclei – quasars: general –
quasars: supermassive black holes

††pubyear: 2024††pagerange: Ensemble X-ray variability of optically selected
QSOs: dependence on black hole mass and Eddington ratio–B

##  1 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 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 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., 2017,
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 debate. The Georgakakis et al. (2021) modeling
favors a bending power-law PSD with a brake 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
points to a PSD with constant amplitude and a brake frequency that depends on
both Eddington ratio and black hole mass as in 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, 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. (2018); Sartori 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<l<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 nature of the X-ray
observations.

Throughout this work, we adopt a flux $\Lambda$-CDM (Cold Dark Matter)
cosmology with $\rm H_{0}=70km\,s^{-1}$ and $\Omega_{\Lambda}=0.7$.

##  2 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., submitted). 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.3 keV 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
flagging 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 $\rm 1283\,deg^{2}$) and the adopted
matching radius the expected chance association rate is 0.4%. We then remove
from the analysis a total of 598 QSOs from the full DRW16Q catalogue that lie
within 45 arcsec off the positions of these 9474 associations. The X-ray
aperture photometry of these 598 QSOs may be contaminated by photons
associated with the wings of the nearby X-ray detections. In summary a total
of 598 SDSS QSOs are rejected because they lie close to but are not physically
associated with X-ray detections. 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 ($\approx 30$ arcsec
HEW). We therefore remove SDSS QSOs that lie in the close vicinity of known
X-ray selected clusters (see Section 4.2.2. of Merloni et al., submitted)
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 trade-off 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\times R_{500}$ in
the case of systems with measured $R_{500}$ parameter, or 500 Mpc in the case
of clusters with no $R_{500}$ estimates, or 1 arcmin for systems with unknown
redshifts. For the XMM-Newton observations the exclusion radius is defined as
$0.15\times R_{500}$ in the case of systems with measured $R_{500}$ parameter,
or 150 Mpc in the case of clusters with no $R_{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)

$R_{F}=\frac{F_{\rm 5\,GHz}}{F_{\rm 4400\,\mathring{A}}}$ (1)

where $F_{\rm 5\,GHz}$ and $F_{\rm 4400\,\mathring{A}}$ are the monochromatic
flux densities at 5 GHz and 4400 Å respectively. For the determination of the
$F_{\rm 4400\,\mathring{A}}$ we first estimate the flux density at rest-frame
$\rm 2500\,\mathring{A}$ using the apparent $i$-band magnitude as explained in
Section 5 of Richards et al. (2006). The flux density at rest-frame $\rm
2500\,\mathring{A}$ is converted to $F_{\rm 4400\,\mathring{A}}$ assuming a
power-law spectral energy distribution of the form
$F_{\nu}\propto\nu^{-0.44}$. For the estimation of $F_{\rm 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., 2017, 2022). The former
covers about 10,000 deg2 and most of the SDSS area to a median 1$\sigma$ rms
sensitivity of 0.141$\rm\,mJy\,beam^{-1}$ at 1.4 GHz. The LoTSS DR2 currently
covers about 5,700 deg2 at 144 MHz to a deeper flux density limit compared to
FIRST with a median 1$\sigma$ rms sensitivity of 83 $\rm\mu 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\sigma$
upper limits to their 1.4 GHz flux density using the FIRST RMS noise
maps111sundog.stsci.edu/first/catalogs/readme.html#coverage. These upper
limits are given by the relation $0.25+3\sigma_{RMS}$, where $\sigma_{RMS}$ 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\sigma$ upper limit to their 144 MHz flux
densities of 0.3 mJy based on the median RMS noise of 83 $\rm\mu
Jy\,beam^{-1}$. For both radio surveys a power-law radio spectral shape with
index $-0.8$ (i.e. $F_{\nu}\propto\nu^{-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 $R_{F}<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 $R_{F}\ga 100$ (e.g Zhu et al., 2020). In this work we
therefore adopt the less conservative threshold of $R_{F}=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 $R_{L}=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 $R_{L}$ upper limits or measurements that
are below the $R_{L}=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 $R_{L}>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. The histogram of the rest-frame time differences between the
two epoch X-ray observations of that sample is shown in Figure 4. The sample
is probing timescales that extend to about 1 decade (rest-frame). Also, the
contribution of the XMM/eROSITA flux measurement pairs increases toward longer
timescales.

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-to-one relation. The
horizontal and vertical dashed lines show the radio loudness threshold
$R_{F}=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: Distribution of the
baseline sample of SDSS DRQ16 QSOs (see Section 2) on the Eddington ratio vs
black hole mass plane. The 1-dimensional projections of 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 the main panel respectively.
Figure 3: Distribution of the baseline sample of SDSS DRQ16 QSOs (see Section
2) on the bolometric luminosity vs redshift plane. The 1-dimensional
projections of of this distribution along the bolometric luminosity and
redshift axes are also shown by the histograms on the right and at the top of
the the main panel respectively. 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.

##  3 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 $\tau$ as

$SF^{2}(\tau)=\langle[\log f_{X}(t)-\log f_{X}(t+\tau)]^{2}\rangle,$ (2)

where $f_{X}(t)$, $f_{X}(t+\tau)$ are the X-ray fluxes of a source measured at
epochs separated by $\tau$, i.e. at times $t$ and $t+\tau$. For a population
of AGN the equation above estimates the variance of the logarithmic flux
differences at the timescale $\tau$. 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)

$SF^{2}(\tau)=\langle[\log f_{X}(t)-\log
f_{X}(t+\tau)]^{2}\rangle-\langle\sigma^{2}_{X}(t)+\sigma^{2}_{X}(t+\tau)\rangle,$
(3)

where $\sigma^{2}_{X}(t)$ represents the uncertainty associated with the flux
measurement $\log f_{X}(t)$. 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 $i$ with X-ray observations in two distinct epochs. In the
first dataset the source has $N_{i,1}$ counts within an aperture with EEF
$EEF_{i,1}$. The background level is assumed to be $B_{i,1}$ and the exposure
time at the source position is $t_{i,1}$. The probability that the source has
flux $f_{i,1}$ can be expressed by the Poisson probability

$P(f_{i,1}\,|\,O_{1})=\frac{e^{-\lambda_{i,1}}\times\,\lambda_{i,1}^{N_{i,1}}}{N_{i,1}!}.$
(4)

Where $O_{1}$ signifies the observations at hand (e.g. $N_{i,1}$, $B_{i,1}$
etc) and the Poisson expectation value $\lambda_{i,1}$ is

$\lambda_{i,1}=f_{i,1}\times ECF_{i,1}\times t_{i,1}\times
EEF_{i,1}\,+B_{i,1}.$ (5)

$ECF_{i,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 $i$ has
$N_{i,2}$ counts within an aperture of EEF $EEF_{i,2}$, the background level
is $B_{i,2}$ and the exposure time is $t_{i,2}$. The probability that the
source’s flux has changed by $\Delta\log
F_{i}=\Delta\log(f_{i,1}/f_{i,2})=\log f_{i,1}-\log f_{i,2}$ to $f_{i,2}$ is
also given by the Poisson formula

$P(f_{i,1}\cdot\frac{f_{i,2}\,}{f_{i,1}}\,|\,O_{2})\,=\,\frac{e^{-\lambda_{i,2}}\cdot\,\lambda_{i,2}^{N_{i,2}}}{N_{i,2}!}.$
(6)

Here $O_{2}$ signifies the second set of observations at hand. The flux
variation is expressed by the ratio $f_{i,1}/f_{i,2}$, which links to the
definition of the structure function in Equation 2. The Poisson expectation
value $\lambda_{i,2}$ is

$\lambda_{i,2}=\left(f_{i,1}\cdot\frac{f_{i,2}}{f_{i,1}}\right)\times
ECF_{i,2}\times t_{i,2}\times EEF_{i,2}\,+B_{i,2}.$ (7)

Where $ECF_{i,1}$ is the energy to flux conversion factor for the second
observation. The likelihood of the two epoch observations of source $i$ is the
product of the Poisson probabilities of Equations 4, 6 and the probability
that a flux ratio variation $f_{i,2}/f_{i,1}$. For the latter, it is assumed
that the variations of the logarithmic flux ratio, $\log f_{i,2}/f_{i,1}$, is
drawn from a normal distribution with mean of zero and variance that depends
on the (rest-frame) time interval $\Delta T_{i}$ between the two observations.
This variance is the structure function at timescale $\Delta T_{i}$, 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 $i$ as

$\mathcal{L}_{i}=P(f_{i,1}\,|\,O_{1})\cdot
P(f_{i,1}\cdot\frac{f_{i.2}}{f_{i,1}}\,|\,O_{2})\cdot\mathcal{N}\left(\log\frac{f_{i.2}}{f_{i,1}}\,|\,\mu,\sigma\right)$
(8)

The last term describes the normal distribution with parameters $\mu$ (mean)
and $\sigma$ (scatter). In the case of $N$ number of sources with two epoch
observations separated by the same time interval $\Delta T$ the likelihood of
the population is

$\begin{split}\mathcal{L}=\prod_{i=1}^{N}\mathcal{L}_{i}=\prod_{i=1}^{N}&P\left(\log
f_{i,1}\,|\,O_{1}\right)\,\cdot\\\ &P[\log
f_{i,1}+\log(f_{i,2}/f_{i,1})\,|\,O_{2}]\cdot\\\
&\mathcal{N}[\log(f_{i,2}/f_{i,1})\,|\,\mu,\sigma]\end{split}$ (9)

The likelihood above is used to perform Bayesian inference on $\sigma$, which
represents the structure function of the sample at a fixed time interval. The
Hamiltonian Markov Chain Monte Carlo code Stan222https://mc-stan.org is used
to sample the likelihood of Equation 9 and produce parameter posterior
distributions.

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.

##  4 Results

###  4.1 The X-ray tructure 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 rest-frame time
scale, $\Delta T$, in Figure 5. The SF monotonically increases with increasing
$\rm\Delta 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 $SF\propto\Delta T^{\gamma}$, we estimate an index
$\gamma=0.12\pm 0.01$. This value translates to a Power Spectrum Density (PSD)
power law index of about –1.2 (e.g. Bauer et al., 2009). Also shown in Figure
5 is the expected structure function of a bending power-law model PSD of the
form (McHardy et al., 2004; González-Martín & Vaughan, 2012)

$PSD(\nu)=A\,\nu^{-1}\,\left(1+\frac{\nu}{\nu_{b}}\right)^{-1},$ (10)

with a constant amplitude $A=2\cdot\nu_{b}\cdot PSD(\nu_{b})=0.02$ and break
frequency, $v_{b}$, that depends on the black hole mass of the accreting
system as

$\nu_{b}=\frac{580}{M_{\rm BH}/M_{\odot}}\;(s^{-1}).$ (11)

This parametrisation, which we will refer to as Model 1, is motivated by the
observations of Papadakis (2004) and González-Martín & Vaughan (2012) that
explore the variability properties of nearby Seyferts as a function of
physical parameters. Simulated light curves are generated from this PSD model
using the method described in Emmanoulopoulos et al. (2013) as implemented in
Sartori et al. (2018); Sartori et al. (2019). The simulated light curves are
then used to estimate the structure function that corresponds to the adopted
PSD model. In this exercise we assume two black hole masses $\log
M_{BH}/M_{\odot}=8,\;9$ that bracket the mode of the distribution shown in
Figure 2. It is clear that the structure function of these simulated light
curves is inconsistent with the measured one for SDSS DR16 QSOs. Modifications
to the PSD Model 1 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 as

$\nu_{b}=\frac{200}{86400}\cdot\frac{L_{\rm bol}}{10^{44}{\rm
erg\,s^{-1}}}\cdot\left(\frac{M_{\rm BH}}{10^{6}\,M_{\odot}}\right)^{-2}\;(\rm
s^{-1}),$ (12)

where $L_{\rm bol}$ is the bolometric luminosity and the black hole mass is
measured in solar units. 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,
$\lambda_{Edd}$, as (Ponti et al., 2012)

$A=2\cdot\nu_{b}\cdot PSD(\nu_{b})=3\times 10^{-2}\cdot\lambda_{Edd}^{-0.8}.$
(13)

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
M_{BH}/M_{\odot}=8,\;9$ and Eddington ratio $\lambda_{Edd}=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 M_{BH}/M_{\odot}=8$. The ensemble X-ray variability properties of
optically selected QSOs at cosmological distances therefore 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.

###  4.2 Black hole mass dependence of the X-ray structure function

We select DRQ16 QSOs in a narrow Eddington ratio interval,
$-1.25<\log\lambda_{Edd}<-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. We further
split this sub-sample into two nearly equal size black hole mass bins at the
threshold $\log M_{BH}/M_{\odot}=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. 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 $M_{BH}$ is
inverted, i.e. the high black hole mass sub-sample is more variable.
Uncertainties in the determination of black-hole masses ($\approx 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 $\Delta T$, $M_{BH}$, $\lambda_{Edd}$
and redshift for the two QSO sub-samples selected by black hole mass. It is
clear from this figure that e.g. the mean $\lambda_{Edd}$ 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. Figure 7 further shows that the high black-hole mass sub-sample
is 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. 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.

###  4.3 Eddington ratio dependence of the X-ray structure function

Next, we select DRQ16 QSOs in a narrow black-hole mass interval, $8.75<\log
M_{BH}/M_{\odot}<9.25$ 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\lambda_{Edd}=-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 ($\approx 10$ yrs) have high Eddington ratio sources a
higher structure function compared to low Eddington ratio ones. Although the
differences between the two sub-samples in Figure 8 are not as conclusive 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
($\log\lambda_{Edd}<-1.5$, $-1.5<\log\lambda_{Edd}<-1$ and
$\log\lambda_{Edd}>-1$) yields qualitatively similar conclusions, i.e. the
structure function tends to be higher for DRQ16 QSOs with low Eddington
ratios. 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 $\Delta T$,
$M_{BH}$, $\lambda_{Edd}$ 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\lambda_{Edd}=-1.3$), it
underestimates the structure function of the high Eddington-ratio sub-sample
(mean $\log\lambda_{Edd}=-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.

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 $\Delta
T_{rest}$ 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 6: Structure function versus rest-frame time-scale for the SDSS
DRQ16 QSO subsample with Eddington ratios in the interval
$-1.25<\log\lambda_{Edd}<-0.5$, which is further split by black hole mass. The
orange circles correspond to SDSS DRQ16 QSOs with $\log M_{BH}/M_{\odot}>8.7$.
Blue circles are for QSOs with black holes lower than the limit above, i.e.
$\log M_{BH}/M_{\odot}<8.7$. The horizontal error bars show the extent of the
rest-frame time interval. Data points are plotted at the mean $\Delta
T_{rest}$ 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 M_{BH}/M_{\odot}=8.5$, while
the orange curves assume $\log M_{BH}/M_{\odot}=9$. The adopted Eddington
ratio for all curves is $\log\lambda_{Edd}=-0.9$. Figure 7: The histograms in
each panel show the distributions of each of the key parameters, $\log\Delta
T$ (top row), $\log M_{BH}$ (2nd row from top), $\log\lambda_{Edd}$ (3rd row
from top) and redshift (last row), for the two subsamples with Eddington ratio
in the range $-1.25<\lambda_{Edd}<-0.5$ and black holes masses in the
intervals $\log M_{BH}/M_{\odot}=8-8.75$ (blue lines) and $\log
M_{BH}/M_{\odot}=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. 6 as
indicated at the top of the column. The vertical lines in each panel show the
mean of the distribution of the corresponding parameter. Figure 8: Structure
function versus rest-frame time-scale for the subsample with $8.75<\log
M_{BH}/M_{\odot}<9.25$, which is further split by Eddington ratio. The blue
circles correspond to SDSS DRQ16 QSOs with $\log\lambda_{Edd}<-1$. Orange
circles are for QSOs with Eddington ratio larger than the limit above, i.e.
$\log\lambda_{Edd}>-1$. The horizontal error bars show the extent of the rest-
frame time interval. Data points are plotted at the mean $\Delta T_{rest}$ 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\lambda_{Edd}=-1.3$, while the
orange curves assume $\log\lambda_{Edd}=-0.7$. The adopted black hole mass for
all the curves is $\log M/M_{\odot}=9$. Figure 9: The histograms in each panel
show the distributions of each of the key parameters, $\log\Delta T$ (top
row), $\log M_{BH}$ (2nd row from top), $\log\lambda_{Edd}$ (3rd row from top)
and redshift (last row), for the two subsamples with black holes masses in the
range $\log M_{BH}/M_{\odot}=8.75-9.25$ and Eddington ratios in the intervals
$-2.0<\log\lambda_{Edd}<-1.0$ (blue lines) and $-1.0<\log\lambda_{Edd}<0.0$
(orange lines). Each column of panels corresponds to one of the rest-frame
times scales of the data points plotted in Fig. 8 as indicated at the top of
the column. The vertical lines in each panel show the mean of the distribution
of the corresponding parameter.

##  5 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 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 ($\approx 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.

There is an increasing body of evidence suggesting that the UV/optical
variability of 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.
On longer time-scales 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 are also
recovered for the X-ray variability of local Seyferts based on either the
normalised excess variance at different time scales (Ponti et al., 2012;
Akylas et al., 2022) or the Power Spectral Density of the observed variations
(McHardy 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 $\nu_{b}$ scales with both black hole mass and
Eddington ratio (see Eq. 12, McHardy et al., 2004, 2006; Paolillo et al.,
2023) 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). We
note however, that such a PSD model is not supported by the recent work of
Paolillo et al. (2023). They combined archival observations of local Seyferts
and high redshift X-ray selected 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. They show that the
observations are consistent with a universal PSD law with a constant amplitude
that is independent of redshift or Eddington ratio and a break frequency that
scales with black hole mass and possibly Eddington ratio. Although further
work is needed to resolve this discrepancy, our analysis demonstrates the
power of ensemble X-ray flux variability studies of hight 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.

##  6 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 turnover in the X-ray structure function to the
longest time-scales of about 10 yr rest-frame probed by our sample. This has
implications on the size of accretion disk and may suggest that 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 increases
with decreasing black hole mass and toward lower Eddington ratios. This
dependence is best represented by PSD models with 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.

## Acknowledgements

The research leading to these results has received funding from the Hellenic
Foundation for Research and Innovation (HFRI) project "4MOVE-U" grant
agreement 2688, which is part of the programme "2nd Call for HFRI Research
Projects to support Faculty Members and Researchers". This work is based on
data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-
German science mission supported by the Russian Space Agency (Roskosmos), in
the interests of the Russian Academy of Sciences represented by its Space
Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt
(DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its
subcontractors, and is operated by NPOL with support from the Max Planck
Institute for Extraterrestrial Physics (MPE). The development and construction
of the eROSITA X-ray instrument was led by MPE, with contributions from the
Dr. Karl Remeis Observatory Bamberg and ECAP (FAU Erlangen-Nuernberg), the
University of Hamburg Observatory, the Leibniz Institute for Astrophysics
Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the
University of Tübingen, with the support of DLR and the Max Planck Society.
The Argelander Institute for Astronomy of the University of Bonn and the
Ludwig Maximilians Universität Munich also participated in the science
preparation for eROSITA. This research made use of Astropy,333www.astropy.org
a community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013, 2018). Based on observations obtained with XMM-Newton, an ESA
science mission with instruments and contributions directly funded by ESA
Member States and NASA.

## Data Availability

The code and data used in this paper are available at GitHub
https://github.com/ageorgakakis/StructureFunction and Zenodo
https://zenodo.org/record/10560969.

## References

  * Adami et al. (2018) Adami C., et al., 2018, A&A, 620, A5
  * Akylas et al. (2022) Akylas A., Papadakis I., Georgakakis A., 2022, A&A, 666, A127
  * Allevato et al. (2013) Allevato V., Paolillo M., Papadakis I., Pinto C., 2013, ApJ, 771, 9
  * Almeida et al. (2023) Almeida A., et al., 2023, ApJS, 267, 44
  * Arévalo & Uttley (2006) Arévalo P., Uttley P., 2006, MNRAS, 367, 801
  * Arévalo et al. (2023) Arévalo P., Lira P., Sánchez-Sáez P., Patel P., López-Navas E., Churazov E., Hernández-García L., 2023, MNRAS, 526, 6078
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Bauer et al. (2009) Bauer A., Baltay C., Coppi P., Ellman N., Jerke J., Rabinowitz D., Scalzo R., 2009, The Astrophysical Journal, 696, 1241
  * Becker et al. (1995) Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559
  * Breedt et al. (2009) Breedt E., et al., 2009, MNRAS, 394, 427
  * Brunner et al. (2022) Brunner H., et al., 2022, A&A, 661, A1
  * Cackett et al. (2021) Cackett E. M., Bentz M. C., Kara E., 2021, iScience, 24, 102557
  * Clerc et al. (2012) Clerc N., Sadibekova T., Pierre M., Pacaud F., Le Fèvre J. P., Adami C., Altieri B., Valtchanov I., 2012, MNRAS, 423, 3561
  * Czerny (2006) Czerny B., 2006, in Gaskell C. M., McHardy I. M., Peterson B. M., Sergeev S. G., eds, Astronomical Society of the Pacific Conference Series Vol. 360, AGN Variability from X-Rays to Radio Waves. p. 265
  * Dey et al. (2019) Dey A., et al., 2019, AJ, 157, 168
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos D., McHardy I. M., Papadakis I. E., 2013, MNRAS, 433, 907
  * Georgakakis et al. (2021) Georgakakis A., Papadakis I., Paolillo M., 2021, MNRAS, 508, 3463
  * González-Martín & Vaughan (2012) González-Martín O., Vaughan S., 2012, A&A, 544, A80
  * Hardcastle et al. (2023) Hardcastle M. J., et al., 2023, arXiv e-prints, p. arXiv:2309.00102
  * Hernández Santisteban et al. (2020) Hernández Santisteban J. V., et al., 2020, MNRAS, 498, 5399
  * Jiang & Blaes (2020) Jiang Y.-F., Blaes O., 2020, ApJ, 900, 25
  * Kalberla et al. (2005) Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M., Bajaja E., Morras R., Pöppel W. G. L., 2005, A&A, 440, 775
  * Kellermann et al. (1989) Kellermann K. I., Sramek R., Schmidt M., Shaffer D. B., Green R., 1989, AJ, 98, 1195
  * Kollmeier et al. (2017) Kollmeier J. A., et al., 2017, arXiv e-prints, p. arXiv:1711.03234
  * Kozłowski (2016) Kozłowski S., 2016, ApJ, 826, 118
  * Kozłowski (2017) Kozłowski S., 2017, ApJ, 835, 250
  * Liu et al. (2022) Liu A., et al., 2022, A&A, 661, A2
  * Luo et al. (2017) Luo B., et al., 2017, ApJS, 228, 2
  * Lyke et al. (2020) Lyke B. W., et al., 2020, ApJS, 250, 8
  * MacLeod et al. (2012) MacLeod C. L., et al., 2012, ApJ, 753, 106
  * Markowitz et al. (2003) Markowitz A., et al., 2003, ApJ, 593, 96
  * Massaro et al. (2015) Massaro E., Maselli A., Leto C., Marchegiani P., Perri M., Giommi P., Piranomonte S., 2015, Ap&SS, 357, 75
  * McHardy et al. (2004) McHardy I. M., Papadakis I. E., Uttley P., Page M. J., Mason K. O., 2004, MNRAS, 348, 783
  * McHardy et al. (2006) McHardy I. M., Koerding E., Knigge C., Uttley P., Fender R. P., 2006, Nature, 444, 730
  * Mehrtens et al. (2012) Mehrtens N., et al., 2012, MNRAS, 423, 1024
  * Menzel et al. (2016) Menzel M. L., et al., 2016, MNRAS, 457, 110
  * Merloni (2019) Merloni A., 2019, in Preparing for 4MOST. p. 23, doi:10.5281/zenodo.3244952
  * Middei et al. (2017) Middei R., Vagnetti F., Bianchi S., La Franca F., Paolillo M., Ursini F., 2017, A&A, 599, A82
  * Nandra et al. (1997) Nandra K., George I. M., Mushotzky R. F., Turner T. J., Yaqoob T., 1997, The Astrophysical Journal, 476, 70
  * Paolillo et al. (2004) Paolillo M., Schreier E. J., Giacconi R., Koekemoer A. M., Grogin N. A., 2004, ApJ, 611, 93
  * Paolillo et al. (2017) Paolillo M., et al., 2017, MNRAS, 471, 4398
  * Paolillo et al. (2023) Paolillo M., et al., 2023, A&A, 673, A68
  * Papadakis (2004) Papadakis I. E., 2004, MNRAS, 348, 207
  * Papadakis, I. E. et al. (2002) Papadakis, I. E. Brinkmann, W. Negoro, H. Gliozzi, M. 2002, A&A, 382, L1
  * Papadakis et al. (2001) Papadakis I. E., Nandra K., Kazanas D., 2001, ApJ, 554, L133
  * Papadakis et al. (2003) Papadakis I. E., Reig P., Nandra K., 2003, MNRAS, 344, 993
  * Papadakis et al. (2008) Papadakis I. E., Chatzopoulos E., Athanasiadis D., Markowitz A., Georgantopoulos I., 2008, A&A, 487, 475
  * Pâris et al. (2017) Pâris I., et al., 2017, A&A, 597, A79
  * Pâris et al. (2018) Pâris I., et al., 2018, A&A, 613, A51
  * Ponti et al. (2012) Ponti G., Papadakis I., Bianchi S., Guainazzi M., Matt G., Uttley P., Bonilla N. F., 2012, A&A, 542, A83
  * Predehl et al. (2021) Predehl P., et al., 2021, A&A, 647, A1
  * Prokhorenko et al. (2024) Prokhorenko S. A., et al., 2024, arXiv e-prints arXiv:2401.12860, p. arXiv:2401.12860
  * Richards et al. (2006) Richards G. T., et al., 2006, The Astronomical Journal, 131, 2766
  * Ruiz et al. (2022) Ruiz A., Georgakakis A., Gerakakis S., Saxton R., Kretschmar P., Akylas A., Georgantopoulos I., 2022, MNRAS, 511, 4265
  * Sánchez-Sáez et al. (2018) Sánchez-Sáez P., Lira P., Mejía-Restrepo J., Ho L. C., Arévalo P., Kim M., Cartier R., Coppi P., 2018, ApJ, 864, 87
  * Sartori et al. (2018) Sartori L. F., Schawinski K., Trakhtenbrot B., Caplar N., Treister E., Koss M. J., Urry C. M., Zhang C. E., 2018, MNRAS, 476, L34
  * Sartori et al. (2019) Sartori L. F., Trakhtenbrot B., Schawinski K., Caplar N., Treister E., Zhang C., 2019, ApJ, 883, 139
  * Secunda et al. (2023) Secunda A., Greene J. E., Jiang Y.-F., Yao P. Z., Zoghbi A., 2023, arXiv e-prints, p. arXiv:2306.05455
  * Shimwell et al. (2017) Shimwell T. W., et al., 2017, A&A, 598, A104
  * Shimwell et al. (2022) Shimwell T. W., et al., 2022, A&A, 659, A1
  * Tang et al. (2023) Tang J.-J., Wolf C., Tonry J., 2023, Nature Astronomy, 7, 473
  * Timlin et al. (2020) Timlin John D. I., Brandt W. N., Zhu S., Liu H., Luo B., Ni Q., 2020, MNRAS, 498, 4033
  * Tortosa et al. (2023) Tortosa A., et al., 2023, MNRAS, 526, 1687
  * Turner et al. (2001) Turner M. J. L., et al., 2001, A&A, 365, L27
  * Uttley et al. (2002a) Uttley P., McHardy I. M., Papadakis I. E., 2002a, MNRAS, 332, 231
  * Uttley et al. (2002b) Uttley P., McHardy I. M., Papadakis I. E., 2002b, Monthly Notices of the Royal Astronomical Society, 332, 231
  * Uttley et al. (2014) Uttley P., Cackett E. M., Fabian A. C., Kara E., Wilkins D. R., 2014, A&ARv, 22, 72
  * Vagnetti et al. (2016) Vagnetti F., Middei R., Antonucci M., Paolillo M., Serafinelli R., 2016, A&A, 593, A55
  * Webb et al. (2020) Webb N. A., et al., 2020, A&A, 641, A136
  * White et al. (1997) White R. L., Becker R. H., Helfand D. J., Gregg M. D., 1997, The Astrophysical Journal, 475, 479
  * Wu & Shen (2022) Wu Q., Shen Y., 2022, ApJS, 263, 42
  * Xu et al. (2022) Xu W., Ramos-Ceja M. E., Pacaud F., Reiprich T. H., Erben T., 2022, A&A, 658, A59
  * Yao et al. (2023) Yao P. Z., Secunda A., Jiang Y.-F., Greene J. E., Villar A., 2023, ApJ, 953, 43
  * Zhu et al. (2020) Zhu S. F., Brandt W. N., Luo B., Wu J., Xue Y. Q., Yang G., 2020, Monthly Notices of the Royal Astronomical Society, 496, 245
  * de Vries et al. (2005) de Vries W. H., Becker R. H., White R. L., Loomis C., 2005, AJ, 129, 615

##  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. (2018); Sartori et al. (2019) based on
the methods of Emmanoulopoulos et al. (2013). For a given rest-frame time-
scale, $\Delta T$, we can then randomly choose $N_{\rm sample}$ number of data
point pairs on the light curve each of which are separated by the time
interval $\Delta T$. These pairs represent the 2-epoch X-ray photometry of a
population of QSOs with size $N_{\rm sample}$ that can be used to determine
their structure function.

Before drawing $N_{\rm 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 11.
We use this historgam to randomly draw $N_{\rm sample}$ fluxes that are
representative of the SDSS DRQ16 QSO sample. These are used to normalise the
mean of the simulated light curve. The end product of this step are $N_{\rm
sample}$ light curves with the same power spectrum but different mean fluxes.
It these $N_{\rm sample}$ renormalised light curves that we use to draw (from
each) two random fluxes separated by the time interval $\Delta T$. 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 $\Delta T$. The inferred SF can
then be compared with the known structure function of the input light curve.

Figure 10 shows the results of this experiment by comparing the inferred SF
with the input one as a function of the sample size, $N_{\rm 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 $\Delta T$
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 10. This figure demonstrates
that when increasing $N_{\rm sample}$ the inferred SF converges to the input
one. Additionally, the inferred SF is always consistent with the input one
within the $1\sigma$ 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.

Figure 10: Inferred structure function from the simulations described in
Appendix A as a function of $N_{sample}$, i.e. the simulated number of 2-epoch
X-ray observation pairs. The squares are the inferred structure functions for
simulations where the photon counts are from XMM-Newton EPIC for one of two
epochs and eROSITA-DE DR1 for the second epoch. The circles show the inferred
structure functions in the case of simulations in which both epochs are
assumed to be observed by XMM-Newton EPIC. The uncertainties associated with
each data point correspond to $1\sigma$ and are determined from the posteriors
of the Bayesian approach described in Section 3. The dashed lines correspond
to 3 different input structure function levels that correspond to either
different time intervals $\Delta T$ or different PSD normalisations. These
input structure function levels are marked with different colors (orange,
blue, green) and should be compared with the inferred structure function data
points of the same colour.

##  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
$i$ of a given XMM-Newton observation for which the RapidXMM database lists
$N_{i}$ total number of photons, a background level $B_{i}$ and an exposure
time $t_{i}$. The probability that this source has flux $f_{X}$ is given by
the Poisson probability distribution function

$P(f_{X})=\frac{e^{-\lambda}\times\,\lambda^{N}}{N!},$ (14)

where $N=\sum_{i=1}^{3}N_{i}$ 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

$\lambda=f_{X}\times\sum_{i=1}^{3}\bigl{(}ECF_{i}\times EEF_{i}\times
t_{i}\bigr{)}+\sum_{i}B_{i},$ (15)

where $ECF_{i}$ is the count to flux conversion factor for the camera $i$ and
$EEF_{i}$ is the encircled energy fraction of the RapidXMM aperture for the
detector $i$ (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 $\Gamma=1.9$ absorbed by the appropriate Galactic hydrogen column density
(Kalberla et al., 2005) in the direction of the SDSS QSO under consideration.
Equation 14 can be solved numerically to determine the X-ray flux probability
distribution function. This then used to determine the mode for each SDSS QSO
and build the distribution plotted in Figure 11.

Figure 11: SDSS QSO 0.2-2 keV X-ray flux distribution inferred from the
RapidXMM aperture photometric products (see text for details).