The universal shape of the X-ray variability power spectrum of AGN up to z ∼ 3

Aims. We study the ensemble X-ray variability properties of Active Galactic Nuclei (AGN) over a large range of timescales (20 ks ≤ T ≤ 14 yrs), redshift (0 ≤ z (cid:46) 3), luminosities (10 40 erg s − 1 ≤ L X ≤ 10 46 erg s − 1 ) and black hole (BH) masses (10 6 ≤ M (cid:12) ≤ 10 9 ). Methods. We propose the use of the variance–frequency diagram, as a viable alternative to the study of the power spectral density (PSD), which is not yet accessible for distant, faint and / or sparsely sampled AGN. Results. We show that the data collected from archival observations and previous literature studies are fully consistent with a universal PSD form which does not show any evidence for systematic evolution of shape or amplitude with redshift or luminosity, even if there may be di ﬀ erences between individual AGN at a given redshift or luminosity. We ﬁnd new evidence that the PSD bend frequency depends on BH mass and, possibly, on accretion rate. We ﬁnally discuss the implications for current and future AGN population and cosmological studies.


Introduction
Flux variability is a defining characteristic of Active Galactic Nuclei (AGN).AGN vary on all timescales, and across the whole electromagnetic spectrum.The fastest and largest amplitude variations are observed at the highest energies (X-rays and γ-rays), strongly suggesting that such radiation is mainly generated in small-size regions, close to the central engine.
One of the most frequently used tools to study the observed variations is the power spectral density function (or power-spectrum for simplicity; PSD).Early studies of the Xray variability of AGN with EXOSAT showed that the PSD has a power-law shape with a slope of ∼ −1.5, and an amplitude that scales with the source luminosity (Green et al. 1993;Lawrence & Papadakis 1993).Long RXTEobserving campaigns over many years combined with shorter XMM-Newton observations (mainly) have allowed the detailed study of AGN X-ray PSDs over a large frequency range, revealing at least one, and in some cases two, breaks in the PSD of nearby AGN (e.g.Uttley et al. 2002;Papadakis et al. 2002;Markowitz et al. 2003;McHardy et al. 2004;González-Martín & Vaughan 2012).These should represent characteristic timescales linked to the physical process producing the observed emission.
Most of our knowledge about AGN power-spectra in the X-ray band is derived from extensive observations of nearby and mostly low-luminosity AGN.It is not possible (yet) to estimate the PSD of AGN at larger redshifts because the available lightcurves have few points and are sparsely sampled.For this reason, our knowledge of the variability properties of the overall AGN population is mainly based on studies of the excess vari-ance as a function of redshift and luminosity, using lightcurves from large samples of X-ray detected AGN in various surveys (e.g Paolillo et al. 2004;Papadakis et al. 2008;Young et al. 2012;Shemmer et al. 2014;Lanzuisi et al. 2014;Yang et al. 2016;Middei et al. 2017;Zheng et al. 2017;Ding et al. 2018;Thomas et al. 2021).
Recently, Paolillo et al. (2017, P17 hereafter) studied the Xray variability properties of distant AGN in the Chandra Deep Field-South region (CDF-S) over 17 years.They used the normalized excess variance σ 2 nxs (i.e. the average lightcurve variance corrected for the noise, see eq. 1 in P17) as a measure of the X-ray variability amplitude of the sources, and they studied the dependence of σ 2 nxs on X-ray luminosity in various redshift bins.They assumed power-spectrum models based on PSD analysis of nearby, bright X-ray Seyferts, and they found that the variability properties of high-z AGN are consistent with a PSD described by a bending power-law, where the bend frequency (and perhaps the PSD amplitude as well) depends on the accretion rate, expressed in terms of the Eddington ratio λ E = Ṁ/ ṀEdd , where Ṁ is the mass accretion rate and ṀEdd is the Eddington accretion rate.
In this work we expand this study, collecting several complementary AGN samples with available excess variance measurements, including the CDF-S (P17), COSMOS (Lanzuisi et al. 2014), CAIXA (Ponti et al. 2012), TARTARUS (O'Neill et al. 2005) and RXTE (Zhang 2011), in order to cover a wide range of timescales, redshifts, luminosities and black-hole (BH) masses.We also estimated the excess variance of numerous additional local AGN, using Swift/XRT and RXTE lightcurves, in order to cover timescales in-between the shortest and the longest ones probed by the σ 2 nxs data from the literature.Our objective is to study the PSD itself using excess variance measurements.In most previous works, σ 2 nxs has been used to investigate the dependence of AGN variability amplitude on Xray luminosity, BH mass and redshift.However the dependence of σ 2 nxs on total lightcurve duration T itself is also important, because the excess variance is (approximately) equal to the integral of the intrinsic power spectrum in the range of frequencies 1/T ≤ ν ≤ 1/(2∆t min ), where ∆t min is the minimum time difference between successive points in the lightcurve 1 (see §2).Due to the close relation between the PSD and σ 2 nxs , we can compute σ 2 nxs from lightcurves with different duration T , and then plot σ 2 nxs as a function of ν T (≡ 1/T ). 2 We refer to the σ 2 nxs vs. ν T plot as the "variance-frequency" plot (VFP).The VFP provides information closely related to the PSD, with the advantage that it can be directly derived for large samples of faint and/or distant AGN, and thus can be studied with the aim of constraining the intrinsic PSD properties.
Although σ 2 nxs is easy to compute, and hence we can create a VFP even from sparsely sampled lightcurves (which are not sufficient to measure the PSD), the use of the VFP is difficult on an individual AGN basis.In fact, given the statistical properties of the σ 2 nxs in the case of a single object we would need many lightcurves in order to estimate, reliably, the intrinsic variance on various timescales (see Allevato et al. 2013, and references therein).
1 Although the excess variance is a biased estimate of the PSD integral, in the case of red-noise PSDs and sparsely sampled lightcurves, it is possible to correct σ 2 nxs for this effect as discussed in section 5 of Allevato et al. (2013). 2 The additional dependence on 1/∆t min is less relevant due to the steep slope of the PSD at high frequencies, and it is anyway fully taken into account in both our simulations, modeling and fitting, as we discuss in detail in the following sections.
On the other hand, we could consider samples of AGN which have been monitored in the same way (i.e.where T and ∆t min are the same for all sources) to compute the mean excess variance, and use it to create the VFP.However, in the case where we use lightcurves of many AGN, one has to consider the dependence of σ 2 nxs on BH mass (M BH ) as well: for a given lightcurve duration, the variance decreases with increasing BH mass (e.g.Papadakis 2004;O'Neill et al. 2005;Ponti et al. 2012).Therefore, assuming we know the BH masses for all AGN in the sample, we must first model the excess variance dependence on M BH to create VFP for AGN at a fixed mass.
Following this approach, the questions we aim to investigate are the following: a) is the measured VFP consistent with the hypothesis that the X-ray PSD has the same form in all AGN (i.e., the hypothesis that the X-ray variability mechanism is the same in all of them), and b) if yes, what are the characteristics of this "universal" X-ray PSD of AGN?
The paper is organized as follows: in §2 we explain the relation between the VFP and the PSD, and how we can measure the VFP for sparsely sampled, low S/N AGN.In §3 we present the BH mass and variability measurements for high redshift sources in CDF-S and COSMOS samples, and the best-fit results to their σ 2 nxs -M BH relations.In §4 we present the best-fit results for low redshift sources from literature or archival data.In §5 and 6 we present the VFP of AGN up to redshift ∼ 3, the method we use to fit the observed VFP and the best-fit results.Finally, in §7, we summarize our work and we discuss the implications of our study.

The PSD vs Variance-Frequency plot
Let us assume that the AGN X-ray PSD follows a relation of the form: where A is the PSD normalization (A = PSD(ν b ) × 2ν b ), and ν b is the bend frequency3 .The PSD thus defined has a logarithmic slope of −1 at low frequencies (ν << ν b ), which steepens to −(1+ s) at higher frequencies (ν >> ν b ).This model PSD is based on the results from PSD studies of nearby, low-luminosity (but Xray bright) AGN (e.g.McHardy et al. 2004).Let us also assume that ν b depends on the BH mass as follows (e.g.McHardy et al. 2006;González-Martín & Vaughan 2012): where B is a constant.According to Eq. ( 4) in Allevato et al. 2013, the excess variance of a lightcurve of duration T will be a measure of the normalized "band" variance, defined as follows4 : σ 2 (ν T , ν max ) = where ν T =1/T max and ν max = 1/(2∆t min ).In the case of equally sampled light curves, with a few missing points (like the RXTE, Swift, ASCA and XMM-Newton light curves we use -see below) the shortest frequency sampled by the data is well defined, with ∆t min = 2∆t, where ∆t is the bin size of the light curves.However, in the case of the unevenly sampled light curves (like CDF-S and COSMOS), ∆t min is not so obvious (see for example the discussion in Appendix D in Scargle 1982).For these light curves, we decided to accept as ∆t min the shortest time difference between successive observations.For a fixed ∆t min , we expect a negative correlation between σ 2 and ν T , i.e. the variability amplitude should increase with increasing lightcurve duration in AGN (due to the red noise nature of the AGN PSD).Equation 3 shows that σ 2 (ν T , ν max ) depends on the shape (and normalization) of the PSD.Consequently, a plot of σ 2 as a function of ν T holds similar information to the PSD.We call the plot of σ 2 (ν T , ν max ) versus ν T , the Variance-Frequency Plot (VFP) of an AGN. 5igure 1 shows the PSD and the VFP plot (upper and lower panels, respectively) for an AGN with a BH mass of 10 8 M .
The PSD is computed using Eq.(1), while σ 2 (ν T , ν max ) is computed using Eq. ( 3), for various PSD parameters.The variance plotted in this figure should be equal to the (intrinsic) variance of lightcurve segments with duration of T = 1/ν and ∆t = 250 s, so that ν max = 4 • 10 −3 Hz.
The figure shows that all the major features of the PSD are apparent in the VFP plot as well.For example, both the PSD and VFP have a power-law like shape at high frequencies, with the VFP being flatter than the PSD (note the different scale of the y-axis of the two panels in Fig. 1).At frequencies lower than ν b , both the PSD and VFP flatten, to a slope of −1 in the case of the PSD, while σ 2 ∝ − ln(ν T ), below ν b .The various lines in the same figure also show that when varying the PSD parameters, the PSD and the VFP shapes also vary in similar ways.
It is always better to estimate the power-spectrum itself, even from a statistical point of view.The statistical properties of the periodogram (the estimator of the PSD) are far superior to the statistical properties of the excess variance as a measure of the intrinsic variance of a single source (see Allevato et al. 2013 for the statistical properties of the latter).However, as we already mentioned in §1, we cannot use the available lightcurves of the high-z AGN to estimate their PSD.On the other hand, we can measure their excess variance on different timescales (i.e., different ν T ), hence we can compute the VFP and, as Fig. 1 shows, infer the intrinsic PSD of the sources.

Measuring the VFP of active galaxies
Ideally, we would need long and short, well sampled lightcurves of an AGN to reliably measure variance on a large range of timescales, and construct the VFP.But that is not possible with the currently available data, specially for high−z sources.The available lightcurves are not good enough (either because they are too sparse and/or their signal-to-noise ratio is too low) to measure the VFP of a single AGN.For that reason, we will follow a different approach to construct the VFP of active galaxies, as we describe below.
We can consider samples of AGN, with known BH mass, and lightcurves with the same duration (T obs ) and bin size (∆t), for all AGN in each sample.We can use the lightcurves to measure the excess variance σ 2 nxs for each AGN in the sample.One way to study the VFP would be to choose objects with the same BH mass in each sample, and then plot σ 2 nxs versus 1/T obs for all of them.However, since the excess variance measurements are not Gaussian distributed and their error is unknown (see Allevato et al. 2013), we will not be able to fit the resulting VFP and compare it with model predictions.In order to overcome this serious problem, we have to average, somehow, the measured exsccess variances.
On the other hand, we cannot simply compute the mean excess variance of all the objects in the sample, if they host different mass BHs, since the excess variance will depend on the BH mass.However, we can take advantage of this property, and produce a σ 2 nxs -M BH plot for all sources in each sample.The excess variance is a measure of the band variance, σ 2 (ν T , ν max ), defined by Eq. ( 3).This equation, together with Eq. ( 2), shows how σ 2 nxs should depend on BH mass, depending on the sampled frequencies: if neither ν max nor ν T are much smaller than ν b , σ 2 nxs should decrease with increasing BH mass, roughly in a linear way, in log-log space; if instead ν max and ν T << ν b , σ 2 nxs will not depend on BH mass (it will depend on ν T and ν max , only; see Eq. ( 3)) and the σ 2 nxs versus M BH plot will be flat.This behavior is illustrated by the models plotted in Figs.3-7 as dotted lines, as discussed further in §3 and §4.
We can fit the σ 2 nxs versus M BH plot with a linear function of the form: where M is the mean BH mass of the sample, and α(T obs ) is the variance of an AGN with a mass of M, computed using a lightcurve of duration T obs .The key point here is that, since α(T obs ) is computed by fitting all the data in the σ 2 nxs versus M BH plot, its distribution will be much closer to the Gaussian distribution, and its error will be known (from the fitting procedure).
Thus, we can consider AGN (with known BH mass) that have been observed by various satellites, over different timescales, T obs , to compute σ 2 nxs , plot σ 2 nxs versus M BH , and fit the data with Eq. ( 4).If M is the same for all samples, then the plot of α(T obs ) versus 1/T obs will be representative of the VFP plot for the AGN with a mass of M. In this way, we can also take advantage of the study of the β(T obs ) versus 1/T obs plot as well.As we showed in §2, for a given M BH (and hence ν b ), the VFP shape, and thus β(T obs ) should vary with T obs .The way it varies depends on the relation between ν b and M BH .In other words, the study of the β(T obs ) versus 1/T obs plot will constrain the parameter B in Eq. ( 2).
We plan to follow the approach we outlined above in order to study the VFP of AGN, as we describe in detail in the following sections.We will use data from various X-ray surveys for high−z objects, as well as lightcurves from pointed observations of nearby objects, in order to construct the VFP of AGN, both at high and low redshift.In this way, we will be able to directly compare the low and high−z objects, and investigate whether their PSDs are the same or not.

Black-hole mass measurements for the CDF-S sources
The first AGN sample we considered is derived from the CDF-S X-ray catalog of Luo et al. (2017), and the CDF-S variability measurements used in this work are described in P17; we refer to those works for specific details.BH mass measurements for CDF-S sources are primarily based on the measurements published by Suh et al. (2015).These are derived from optical/nearinfrared spectroscopic measurements of Hα, Hβ and Mg-II line widths of X-ray sources in the E-CDF-S region.BH masses were obtained from scaling relations of M BH with FWHM and luminosity of the broad-line components in the spectra.We refer to the original paper for details on the method and the assessment of the reliability of the masses.The median uncertainty is claimed to be ∼ 0.1 dex with an additional 0.3 dex due to calibration uncertainties in the scaling relations.This data set was integrated with Hα BH mass measurements from Schulze et al. (2018), after re-calibrating the scaling relation to the same one adopted in Suh et al. (2015, their Eq.( 1)), and with Mg-II BH mass measurements from Schramm & Silverman (2013).In total, we collected masses for 40 sources: 35 from Suh et al. (2015), 3 from Schulze et al. ( 2018) and 2 from Schramm & Silverman (2013).In the case of multiple BH mass measurements for the same source, we used the average value; however our results do not change if we adopt individual BH mass measurements instead (giving priority to Hα measurements).For the subsequent analysis we defined a "robust sample" of 15 sources with available M BH measurements, average S /N > 0.8 per point and more than 90 points in their X-ray lightcurves in order to have reliable σ 2 nxs measurements.In fact, we verified that including more sources with lower quality lightcurves does not yield significant improvements in the final analysis and increases the uncertainties on the measured variability (see P17 for further details).
In Fig. 2 we plot the luminosity-redshift distribution of CDF-S sources with available M BH measurements and those in the robust samples.The sources in the robust sample are distributed between ∼ 0.5-1.7 in redshift, ∼ 5 × 10 42 -2 × 10 44 erg s −1 in (2-8 keV rest-frame) X-ray luminosity, 10 7 -10 9 M in BH mass, and have low absorption N H < 2 × 10 22 cm −2 as expected from their type 1 spectra.

The CDF-S σ 2
nxs -M BH relation We created σ 2 nxs -M BH plots for the CDF-S data using the excess variance measurements of P17 and the BH mass measurements for the sources in the "robust sample" discussed in the previous section.P17 calculated the excess variance of the CDF-S AGN using T obs = 45, 128, 654 and 6005 day-long (observer's frame) lightcurves (see the discussion in their §4.2).Here we used only their excess variance measurements from the 128 and 654 day-long lightcurve segments, i.e. the segments plotted in the two rightmost panels in Fig. 1 of P17.They are separated by almost 5 years, and the excess variance measurements based on them should be uncorrelated.On the contrary, the σ 2 nxs measurements from the longest and shortest timescales use lightcurve parts which overlap with each other (see Fig. 1 in P17), hence the resulting σ 2 nxs will be correlated.Note that since the AGN PSD is known to depend on the energy (e.g.McHardy et al. 2004), as well as to minimize the effect of absorption, the lightcurves are extracted in the 2-8 keV rest-frame band.
The panels in Fig. 3 show the σ 2 nxs -M BH plots for the two different timescales.A weak anti-correlation between σ 2 nxs and M BH can be observed, although the scatter is considerable.This is due to several reasons.First, the error on M BH and, more importantly, on the individual σ 2 nxs measurements (Allevato et al. 2013), will introduce considerable scatter around the intrinsic relation.In addition, although T obs is the same for all objects, the lightcurve variance depends on the duration in rest-frame, T rest , which is not the same for all objects, as their redshifts differ.We note that, as we explained in §2, if both the duration and bin size of the lightcurve are much longer than the bend timescale, then the intrinsic σ 2 nxs -M BH relation will be flat (unless the PSD amplitude depends on M BH ).In other words, a lack of correlation between excess variance and BH mass is representative of an intrinsically flat σ 2 nxs -M BH relation, thus holding important information regarding the shape of the VFP and the PSD at low frequencies.
We cannot use χ 2 to fit the data plotted in Fig. 3 (as well as in all figures which show σ 2 nxs -M BH relations) and then test whether the model fits the data well, or not, because the excess variance measurements are not Gaussian (and in any case their error is unknown).For that reason, we assumed that a straight line provides a good fit to the data, and we used the ordinary least-squares regression of Y on X (OLS(Y|X)) method of Isobe et al. (1990) to fit the data in Fig. 3 (in log-log space), with the linear model defined by Eq. ( 4), with M = 10 8 M , which is similar to the average M BH in the CDF-S sample.As such, the line normalization will be best determined at this mass.In this case, α(T obs ) corresponds to the intrinsic variance of an AGN with M BH = 10 8 M , when measured from a lightcurve of duration T obs .
Best-fit results are listed in Table 1 and black solid lines in Fig. 3 show the best-fit linear models.Not surprisingly, given the large scatter of the points around the best-fit lines, the error on the slope is large; however, the error on the normalization is small.This is due to the fact that, given the rather flat best-fit slope, α(T obs ) is representative of the mean excess variance of all the points in each plot, which is reasonably well determined from the 15 measurements at each timescale.

The COSMOS σ 2
nxs -M BH relation Lanzuisi et al. (2014) used XMM observations in the COSMOS field over a period of ∼ 3.5 years to study the long-term X-ray variability of a large sample of AGN.We used the data plotted in their Fig. 5 to fit the σ 2 nxs -M BH relation.To improve the accuracy of the measured variances, we selected only objects with at least three points in their lightcurves, and with a total rest-frame duration between 100 and 560 days (see their Fig. 1).Since the final sample spans a wide range in T obs , we further divided it in two bins based on the rest-frame lightcurve length: 100 days ≤ T rest < 330 days and 330 days ≤ T rest < 560 days, with a median duration Trest of 240 and 413 days respectively.There are 82 AGN in both groups, with a median redshift of 1.5 and 1.0, in the first and second group, respectively (see Fig. 2).Their X-ray luminosity ranges from 6 × 10 42 to 3 × 10 45 erg s −1 in the 2-10 keV band.We fitted both σ 2 nxs -M BH plots with the model defined by Eq. ( 4), and the same OLS(Y|X) routine as above (Fig. 4).The timescales and best-fit results for the COS-MOS data are listed in Table 1.

The variance-BH mass relation of low-redshift AGN
4.1.Compilation of σ 2 nxs -M BH relations from literature In order to sample shorter and longer timescales, which are not accessible for high-redshift AGN, we used both published and archival data.We first considered the σ 2 nxs -M BH data from the CAIXA sample (Ponti et al. 2012).These authors presented the results from a systematic study of the excess variance of a large AGN sample using XMM-Newton lightcurves.We used

Survey
T T obs and ∆t min,obs represent the maximum and minimum sampled timescales in the observer reference frame, while Trest is the median value of the maximum restframe timescale over all sources; in square brackets we quote the T rest range for high-z samples (differences in T rest for the low-z samples are not significant).α(T obs ) and β(T obs ) represent the best-fit intercept and slope of Eq. ( 4).their σ 2 nxs measurements (2-10 keV) from the lightcurves with T obs = 80 ks.They measured the excess variance on three shorter time timescales as well, but the use of the same lightcurves when measuring σ 2 nxs on different timescales would imply that their σ 2 nxs measurements would be heavily correlated (for the same source).We constructed the respective σ 2 nxs -M BH plot (see the bottom panel of Fig. 5) using the data from their "Rev" AGN sample. 6In doing so, we updated their BH mass estimates with the measurements listed in the AGN BH mass database (Bentz & Katz 2015).There are 11 radio-quiet AGN in this sample with excess variance measurements based on 80ks-long lightcurves.For consistency, we fitted the CAIXA σ 2 nxs -M BH plot using Eq. ( 4), with M = 10 8 M , and the same OLS(Y|X) routine that we used to fit the respective CDF-S plots.The best-fit results are listed in Table 1, and the black solid line in the lower panel of Fig. 5 shows the best-fit line.
In order to get information on shorter timescales, we considered the excess variance measurements from the TARTARUS sample of O'Neill et al. (2005).They used ASCA40 ks long lightcurves, and measured the excess variance of nearby, X-ray bright Seyferts (in the 2-10 keV band).We chose 16 (radio-quiet) AGN from their sample with BH mass measurements based on the reverberation mapping technique, as listed in the database of Bentz & Katz (2015).We added six sources from the CAIXA sample for which Ponti et al. (2012) provide 40 ks σ 2 nxs measurements (these are: NGC 4151, Mrk 110, Mrk 279, Mrk 590, NGC 4593 and PG 1211+143); they are not part of the 80 ks sample and have BH mass estimates based on reverberation mapping.The respective σ 2 nxs -M BH plot is shown in the top panel of Fig. 5.As above, we fitted the data using Eq. ( 4) and the same OLS(Y|X) routine that we used to fit the respective CDF-S plots.Best-fit results are listed in Table 1, and the solid line in the upper panel of Fig. 5 indicates the best-fit line.
We also considered the σ 2 nxs -M BH data from the RXTE lightcurves of Zhang (2011), to get information on longer timescales.They used data from the ASM on board RXTE to study the X-ray variability amplitude of 27 AGN over a period of T obs = 14 years.Using their σ 2 nxs and M BH measurements, we fitted the resulting variability-M BH relation with the same model and fitting routine as above.The σ 2 nxs -M BH relation, together with the best-fit line, are shown in Fig. 6 while the best-fit re-

Medium-term σ 2 nxs measurements from archival observations
There is a considerable gap between the duration of the shortest CDF-S lightcurves (128 days) and the ∼ 1 day-long CAIXA lightcurves.In order to measure the excess variance on intermediate timescales, we used RXTE/PCA and Swift/XRT data of nearby AGN, and we computed their excess variance using lightcurves which are 9.45 days long (rest-frame).This timescale   2. Symbols have the same meaning as in Fig. 3. is ten times longer than the CAIXA lightcurves and ∼ 10 times shorter than the CDF-S.The sources, together with information on the lightcurves we used and the resulting excess variance, are listed in Table 2.
Column (2) lists the start/end time of the lightcurves we used (in MJD) and, in parentheses, the (minimum) lightcurve bin size ∆t min , in days.We did not reduce the Swift data ourselves; instead we used the lightcurves from Edelson et al. (2019) We divided the lightcurves into segments with a (rest-frame) duration of 9.45 days (as determined by the shortest lightcurve in the sample).We computed the excess variance of each segment in the usual way (i.e., using Eq.(1) in P17).We then computed the mean of the individual σ 2 nxs measurements for each source, which is listed in the last column of Table 2.As for all the other samples, we fitted the σ 2 nxs -M BH relation using Eq. ( 4) and the OLS(Y|X) routine.The result is shown in Fig. 7, together with the best-fit model, while the best-fit results are listed in Table 1.

The observed VFP of AGN
As we already explained, the best-fit α(T ) values listed in Table 1 are representative of the variance of a 10 8 M AGN, on the timescales that are listed in the second column of the same table.We therefore used the α(T ) values listed in this table and we constructed the VFP for the 10 8 M AGN. Figure 8 shows the best-fit α(T ) and β(T ) values plotted as a function of ν T (top and bottom panels, respectively).The two timescales probed by the CAIXA+TARTARUS data (filled squares) constrain the high frequency part of the VFP, and the Swift+RXTE data (filled stars) allow us to sample intermediate timescales.The results from CDF-S (filled circles) provide information on long timescales, while the COSMOS and the RXTE points (filled diamonds and filled triangle, respectively) further improve the accuracy of the VFP at low frequencies.
The limits of the PSD integral in Eq. ( 3) depend on restframe T max and ∆t min .All AGN in the CAIXA, TARTARUS, RXTE and Swift samples are nearby and all lightcurves are of the same duration, hence T obs = T rest (and ∆t min,obs = ∆t min,rest ) in these samples.Things are more complicated in the CDF-S and COSMOS samples.T obs is the same for all AGN in the CDF-S; however, T rest is not, because they do not have the same redshift.In the case of the COSMOS AGN, even T obs differ for different sources due to the survey strategy (see Lanzuisi et al. 2014).Due to differences in T rest in these samples, in Fig. 8 we plot all quantities as a function of ν T = 1/ Trest , where Trest is the median lightcurve duration in the AGN rest-frame ( Trest values are listed in the fourth column of Table 1).Note that this choice is for visualization purposes only, since in the fitting procedure we properly take into account the actual rest-frame timescales sampled for each individual source (see §,6).
The overall VFP (upper panel in Fig. 8) appears to be very well described by a single function, from the lowest to the highest sampled frequencies.We can see the logarithmic The observed VFP of a 10 8 M AGN as a function of the rest-frame frequency ν T , using our measurements, from Table 1 (filled symbols).The red empty symbols indicate the best-fit model values (α(T max ) in Eq. ( 4)).For reference we also show the theoretical expected trends for a 10 8 M BH monitored with an average ∆t min = 250 s (dashed line) and ∆t min = 300 days (dotted line), but we stress that the comparison should be done between the filled and empty points.Lower Panel: Measured and best-fit model slope of the σ 2 nxs -M BH relations (β(T max ) in Eq. ( 4)); symbols have the same meaning as those in the upper panel.
rise of the variance with decreasing frequency (i.e., increasing timescale), with a slope which is roughly equal to −1, from the CAIXA+TARTARUS to the high-frequency CDF-S point.This implies a PSD slope of ∼ −2 at high frequencies.The variance-frequency slope flattens at lower frequencies, indicating the presence of a bend frequency, analogous to the PSD bend frequency, ν b , somewhere between (1-20 days) −1 .The low and high-frequency parts of the observed VFP are determined by the high−z and low−z AGN samples, respectively, but the important observational result is that the low-frequency VFP appears to be the continuation of the high-frequency VFP, without any hints of a normalization mismatch between the two parts.
At low frequencies, the RXTE variance may appear to underestimate the value expected from a simple extrapolation of the CDF-S and the COSMOS measurements in the VFP, although T rest of the RXTE lightcurves is longer than the rest-frame duration of the longest CDF-S lightcurves.In principle, this could suggest that the PSD normalization at low frequencies is not the same in the distant and local AGN.However, ∆t min in the RXTE lightcurves is also significantly larger than the (rest-frame) minimum timescale in all other lightcurves.In fact, it is almost certainly much longer than the PSD bend frequency of a 10 8 M AGN.This implies that we are missing a significant part of the intrinsic variance in the RXTE lightcurves, hence the smaller α(T ) value.
The bottom panel in Fig. 8 shows how the slope of the σ 2 nxs -M BH relation varies with frequency.As with the VFP plot, the best-fit slopes at high frequencies appear to connect smoothly, without any normalization discontinuities, with the best-fit β(T ) values at lower frequencies.The slope of the σ 2 nxs -M BH plots at low frequencies (long timescales) approaches zero, i.e. the bending frequency of the 10 8 M AGN is probably higher than both the mean bin size and duration of the available lightcurves in the RXTE sample and, to some extend, in the CDF-S and the COSMOS sample.

Model fitting procedure and best-fit results
We considered a grid of A, B, and s values, and for each (A, B, s) combination we used Eqs.( 2) and (3) to compute the variance σ 2 for the AGN in the CDF-S, COSMOS, CAIXA+TARTARUS, RXTE and Swift+RXTE samples.To this end we derived ν T , ν max and ν b using the M BH and z values of each source, and the timescales T obs and ∆t min,obs in Table 1 (we assumed z = 0 for the AGN in the CAIXA+TARTARUS, RXTE and the Swift+RXTE samples which only contain nearby AGN).
In order to properly fit the data we must take into account the biases due to a sparse and/or irregular sampling of the lightcurves.To this end we used Eq. ( 11) in Allevato et al. (2013) to compute the bias affecting the measured σ 2 nxs and we include the same factor in the model variance σ 2 .The bias correction depends on the PSD slope itself, so we used the "average" slope (from Eq.( 1)) over the range of sampled rest-frame timescales for each source.Furthermore we assumed a sparse sampling for the COSMOS sources, and a continuous sampling pattern for the remaining samples.
In this way, for each model parameter combination we ended up with a pair of model (σ 2 , M BH ) values, for each source and timescale in the CDF-S, COSMOS, the CAIXA+TARTARUS, and the Swift+RXTE samples.If the observed σ 2 nxs measurements were Gaussian variables with known errors, then we could use χ 2 statistics to fit all the observed σ 2 nxs versus M BH data plotted in Figs. 3 -7 with the model (σ 2 , M BH ) curves.However this is not the case, so we cannot fit directly the model to the observed excess variance -BH mass plots.Instead, we followed a different way to fit the model to the data.
For each model parameter combination we fit the resulting σ 2 -M BH points with Eq. ( 4), with M = 10 8 M , using the same OLS(Y|X) routine that we used to fit the data.In the end, each set of model parameters, (A, B, s), would result in two best-fit values, α mod (A, B, s, T obs ) and β mod (A, B, s, T obs ), for each σ 2 -M BH relation.
The overall model fit: Based on what we discussed above, each set of model parameters, (A, B, s), would result in eight α mod (A, s, B, T obs ) and β mod (A, s, B, T obs ) values, i.e. one for each of the observed σ 2 nxs -M BH relations plotted in Figs. 3 -7.To get the best-fit model, we minimized χ 2 , defined as follows: where δ[α(T obs,i )] and δ[β(T obs,i )] are the errors on the best-fit normalization and slope values, α(T obs ) and β(T obs ), listed in Table 1.Defined in this way, the best-fit model is the one that minimizes the differences between the model and the observed VFP (i.e. the points plotted in the upper panel of Fig. 8), as well as the model and the observed slope of the σ 2 nxs -M BH relations (bottom panel of Fig. 8).
Hz −1 , B = 3.4 +3.1 −1.4 × 10 −6 Hz, and s = 1.7 +0.9 −0.4 .Here the errors are the 90% uncertainties for two interesting parameters.The errors of the best-fit parameter values are both asymmetric and correlated as shown in Fig. 9. Open diamonds in Figs. 3 -7 show the best-fit, model (σ 2 ,M BH ) predictions for each source in these plots.The dashed red lines show the best linear fit to the model points.Clearly, in some cases (e.g.Figs. 5 and 7) the model predicts a curved σ 2 nxs -M BH relation (dotted curve), which actually may be closer to the observed one.But, even in these cases, a straight line can provide a reasonably good fit to the data (and the model predictions).Since we fit a straight line to both the observed and the model relations, we can compare the model predicted and the observed best-fit line parameters, and search for the parameter values which provide the best agreement between data and the model predictions.The open symbols in Fig. 8 show the normalization and slope of the best-fit lines to the model σ 2 nxs -M BH relations that are closest to the data.Strictly speaking, the VFP plotted in the upper panel of Fig. 8 is representative of the average VFP of an AGN with a mass of 10 8 M .The best-fit parameters also refer to such an object.For example, if we had normalized the σ 2 nxs vs. M BH relation at a different M BH , then the expected VFP would be different (see the black and red lines in Fig. 1) and the best-fit bending frequency would be different if ν b depended on BH mass.However, the fact that we consider the VFP of an AGN at 10 8 M is merely a "technicality" as this mass is close to the mean BH mass of the sources in the samples we considered.In reality, the best-fit α mod (A, B, s, T obs ) and β mod (A, B, s, T obs ) model parameters are computed by fitting a straight line to all the model points (i.e.open red points in Figs. 3 -7), Hence, the best-fit model VFP plotted in Fig. 8 holds information about the PSD of all AGN in each sample, and the same would be true if we had normalized the best-fit lines to another BH mass.
We find it impressive that a single PSD model can fit the observed VFP so well, given that we constructed σ 2 nxs -M BH plots using lightcurves of 160 AGN, nearby and distant, observed with different satellites, with different sampling and duration.This result strongly suggests that, on average, the X-ray PSD, over five orders of magnitude in frequency (i.e. from timescales of ∼ 40 ks up to ∼ 10 − 15 years) is described by the same form for all AGN at z ≤ 2 − 3.If there were significant differences between the amplitude, and/or bend frequency, or high-frequency slope of the high−z and low−z PSDs, then the low and high frequency parts of the VFP plotted in the top panel of Fig. 8, which are determined by the high−z and low−z AGN respectively, would differ significantly (see for example the various curves in Fig. 1).

Summary and discussion
We used excess variance measurements computed using short (CAIXA+TARTARUS, Ponti et al. 2012;O'Neill et al. 2005), intermediate (Swift+RXTE, this work) and long (RXTE, Zhang 2011) timescale lightcurves, as well as lightcurves from the COSMOS and the CDF-S surveys (Lanzuisi et al. 2014, P17), to construct σ 2 nxs -M BH plots on various timescales.We fitted them with a simple linear model (in the log-log space), and we studied the resulting variance-frequency plot (VFP), together with the slope of the variance-BH mass relations as a function of frequency.Our main result is that the hypothesis of a common Xray PSD form in all AGN, which remains the same irrespective of redshift and luminosity, is fully consistent with the observed VFP.
The variance vs. timescale relation of the CDF-S sources was used in P17 as well (see their Fig.7), to show that the excess variance measurements were indeed consistent with the assumption of a bending power-law PSD, and in Zheng et al. (2017) to study the low-frequency PSD slope.In this work, we consider a much larger data set to create a more detailed VFP, and we use it together with the slope of the variance -BH mass plots, to actually constrain the X-ray PSD.
We use excess variance measurements from over a hundred AGN, with a wide range of BH mass (∼ 10 6 − 10 9 M ) and X-ray luminosity (∼ 10 40 − 10 46 erg s −1 ) to determine the VFP over timescales from a few hours up to 5 -14 years.Based on the fact that the VFP and the PSD hold the same information, we were able to determine the average PSD of the AGN in our sample.Detailed PSDs have been determined (and fitted with models) only for nearby AGN.To the best of our knowledge, this is the first time that the average X-ray PSD of a representative sample of (X-ray selected) AGN is accurately determined, up to a redshift of ∼ 3.
The fact that the observed VFP is well fitted by a single function is remarkable.The average variance, α(T ), has been measured using lightcurves of different objects, obtained from different instruments, at different times, and over entirely different timescales.Our results thus strongly suggest that the shape of the average PSD (defined by eq. 1) is the same in all AGN, irrespective of luminosity and/or redshift, and is consistent with almost all nearby AGN.We note that the PSD shape of least one nearby AGN, namely Ark 564, is different.It has a power law shape and shows two breaks, at high and low frequencies (Papadakis et al. 2002;McHardy et al. 2007).This shape implies that Ark 564 may be in a state which is equivalent to the so-called "high/intermediate state" in black hole X-ray binaries.Our results suggest that such a state should be rare among AGN.
Regarding the low-frequency slope we assumed a fixed value.If we consider a more general extension of Eq. (1), i.e.

PSD(ν) = Aν
, this means fixing l = −1.On the other hand Zheng et al. (2017) suggested that a steeper low-frequency slope l −1.2 is more appropriate to fit the CDF-S data.In order to test this possibility we repeated all our fits using the more general expression above; we find that our data are better fit by the canonical slope of −1 and, while we cannot rule out a marginally steeper low-frequency slope, a value as steep as -1.2 is excluded at the 99% significance level.
Our results are not meant to imply that all AGN should have the same high frequency slope and PSD amplitude.Most likely, both the PSD amplitude and high frequency slope will be distributed over a range of values, and the best-fit values we report should be indicative of the means of these distributions.For example, González-Martín & Vaughan (2012) studied in detail the X-ray PSD in X-ray bright, local AGN and their results do show a rather broad range of high-frequency slopes (between 1.8 and 4.6, see their Table 4).Similarly, their best-fit PSD amplitude values (also listed in their Table 4) show a broad range of values, from 0.001 to 0.04.What is interesting is that the means of these distributions are the same as our results.
Indeed, the best-fit, high frequency PSD slope from the modeling of the mean VFP is −(1 + s) = −2.7 +0.9 −0.4 (90% confidence limits).Although the PSD slope may depend on the energy band (e.g.McHardy et al. 2004), we point out that in our case we used whenever possible the hard X-ray band (i.e., E 2 keV) and even in the COSMOS sample the large median redshift of the sources implies that we tend to sample energies 1 keV.The median of the high frequency PSD slopes reported by González-Martín & Vaughan (2012) in their Table 4 is -2.57, for sources where a bend frequency has been detected8 .This is consistent with our results.
In addition, our best-fit PSD normalization of A = 0.016 +0.002 −0.003 Hz −1 , is also consistent with the mean PSD normalization as determined by the PSD modeling of local AGN.According to Eq. ( 1), the PSD amplitude at the bend frequency, in terms of ν b × PSD(ν b ), is equal to A/2.Based on our best-fit results, this is 0.008 ± 0.001, which is fully consistent with the mean PSD amplitude of ∼ 0.009 reported by González-Martín & Vaughan (2012) 9  The broad-band VFP plotted in the upper panel of Fig. 8 shows a clear flattening below ∼ 10 −6 Hz, which implies (from our best-fit result) a PSD low frequency bend at b = 3.4 +3.1 −1.4 × 10 −6 Hz.This corresponds to a bend timescale of (1.8−5.8)days (90% confidence).According to González-Martín & Vaughan (2012), log(T b ) log(M BH ) −1.7, where T b is in days and M BH is the BH mass in units of 10 6 M .10For a 10 8 M AGN, this relation predicts T b = 2 days, which is consistent with our results.On the other hand, according to McHardy et al. (2006), the bend timescale should also depend on the source luminosity.These authors find that log(T b ) = 2 × log(M BH ) − log(L bol ) − 2.33, where T b is in days, M BH is the BH mass in units of 10 6 M , and L bol is the bolometric luminosity in units of 10 44 erg s −1 . 11We have assumed that the constants A and B in the McHardy et al. (2006) equation are equal to 2 and 1, respectively, which means that the bend timescale is proportional to the BH mass and inversely proportional to accretion rate (in units of the Eddington limit).For T b = 3.5 days, as derived in this work, and M BH = 10 8 M , then L bol = 1.3 × 10 45 erg s −1 , which is 10% of L Edd .The VFP analysis thus allows for a dependence of the PSD bend frequency on the AGN accretion rate as well.
We note that our results suggest that the average bend frequency and PSD amplitude are the same for the high−z and the low−z sources.This implies that either these two PSD parameters do not depend on accretion rate, or that the accretion rate is the same for both the nearby and the distant AGN in our sample.P17 found that the average accretion rate of the CDF-S sources is ∼ 0.05−0.1 of the Eddington limit.Regarding the low−z objects, if we use the data listed in Table 1 of (Zhang 2011) for the RXTE sample (which are representative of all our low−z AGN) we find an average accretion rate of ∼ 0.06 of the Eddington limit, which is comparable with the accretion rate of the CDF-S sources.We will need to study the VFP of AGN with significantly different accretion rates, in order to investigate the dependence of the PSD amplitude and bending frequency on the accretion rate.
Our results should be useful in future variability studies of large AGN samples, using X-ray lightcurves from, e.g, the eROSITA all-sky survey, and future surveys that may be conducted with the eXTP (in't Zand et al. 2019), Einstein probe Yuan et al. (2018Yuan et al. ( , 2022) ) and Star-X (Saha et al. 2017) proposed missions, as well as Athena (Nandra et al. 2013), Lynx (Gaskin et al. 2018) or AXIS (Mushotzky 2018), provided that deep surveys are properly planned to probe the time domain as well (see e.g.Paolillo et al. 2012).Previous studies, such as P17, assumed that the shape of the X-ray PSD was the same both in nearby and distant, luminous AGN.We show in this work that this is indeed the case.Thus our result can be useful to any study that involves the modeling of the ensemble X-ray variability of AGN in order to, e.g., use them as cosmological probes (La Franca et al. 2014;Lusso et al. 2019Lusso et al. , 2020;;Demianski et al. 2020), or to constrain the AGN demographics through their ensemble X-ray variability properties (Sartori et al. 2019;Georgakakis et al. 2021).From a more physical point of view, a common X-ray PSD shape implies that the same variability mechanism operates in all luminous AGN, and that the mechanism does not evolve with time until up to at least z ∼ 2 − 3.This result is also in agreement with the lack of evolution observed in AGN spectral features, such as the UV vs. L X ratio α ox (e.g.Lusso et al. 2010Lusso et al. , 2020) ) or the Xray spectral slope Γ (e.g.Young et al. 2009), and indicates that it is the underlying X-ray emission mechanism in general that does not evolve with time.Although we do not have a well developed physical model for the X-ray emission and variability in AGN, these implications put another observational constraint for any future attempts.

Fig. 1 .
Fig. 1.The PSD and the VFP (upper and lower panels, respectively) of an AGN with a BH mass of 10 8 M .The PSD model parameters are: A = 0.02 Hz −1 , s=1, and ν b = 3 • 10 −6 Hz (black solid lines).Black dashed lines indicate the PSD and VFP when when we increase the PSD amplitude by a factor of 2.Red lines show the PSD and VFP when we increase the bend frequency (by a factor of 10), while the blue lines show the changes when we increase the high frequency PSD slope from -2 to -3 (all the other parameters in this case are like those of the red lines).For the computation of the variance, we assumed ν max = (1/250 s).

Fig. 2 .
Fig. 2. Luminosity -redshift distribution of CDF-S and COSMOS sources.Solid symbols represent sources with available BH mass measurements while empty ones highlight the final samples used in this work (see text for details).

Fig. 3 .
Fig.3.Plots of σ 2 nxs vs. M BH for the CDF-S "robust" sample (filled circles) for the two different timescales discussed in the text.The solid lines show the linear fit to the data.Open diamonds show the best-fit model predicted variance, computed by taking into account the restframe T rest , and ∆t min,rest of each source.The dashed lines mark the best fit to the model predictions (see §6 for details).The dotted curves represent the model trends for the sources at z = 0.5 (top one) and z = 2.0 (bottom one), i.e the approximate redshift range of the CDF-S sample (see Fig.2).

Fig. 7 .
Fig. 7. Plots of σ 2 nxs vs. M BH for Swift+RXTE sources in Table2.Symbols have the same meaning as in Fig.3.
, for all sources, and Cackett et al. (2020) for Mrk 142.The Swift lightcurves are in the 1.5-10 keV band, except for Mrk 142, where the published lightcurve is in the 0.3-10 keV band.All the other lightcurves were taken from the RXTE AGN Timing & Spectral Database 7 (Rivers et al. 2013) and they are in the 2-10 keV band.BH mass estimates (from Bentz & Katz 2015) are listed in column (3).

Fig. 8 .
Fig. 8. Upper panel:The observed VFP of a 10 8 M AGN as a function of the rest-frame frequency ν T , using our measurements, from Table1(filled symbols).The red empty symbols indicate the best-fit model values (α(T max ) in Eq. (4)).For reference we also show the theoretical expected trends for a 10 8 M BH monitored with an average ∆t min = 250 s (dashed line) and ∆t min = 300 days (dotted line), but we stress that the comparison should be done between the filled and empty points.Lower Panel: Measured and best-fit model slope of the σ 2 nxs -M BH relations (β(T max ) in Eq. (4)); symbols have the same meaning as those in the upper panel.

Table 1 .
Results of the linear fits to the σ 2 nxs -M BH relation (in log-log space), for the different datasets used in this work.

Table 2 .
The RXTE and Swift lightcurves we used to measure the excess variance of AGN on timescales of ∼ 10 days.