Tracing the accretion history of supermassive black holes through X-ray variability: results from the ChandraDeep Field-South

Abstract

1 INTRODUCTION

Flux variability is a defining characteristic of active galactic nuclei (AGNs), reflecting the small spatial region in which the observed emission is produced and the production mechanism itself (Fabian 1979; Rees 1984). AGNs are observed to vary on all time-scales, and across the whole electromagnetic spectrum, although the maximum power and fastest variations are found at the highest energies (X-rays and γ-rays), due to the fact that such radiation is mainly generated close to the central engine and over small spatial regions (Ulrich, Maraschi & Urry 1997).

In the X-ray band in particular, extended and detailed observations of nearby sources have revealed that the variability is characterized by ‘red noise’ behaviour, with more power existing on the longest time-scales, in close resemblance to what is observed in binary accreting systems containing smaller, stellar-mass black holes (BHs). The origin of the X-ray variability itself is not well understood, and both internal (instabilities of the accretion flow, flaring corona, orbiting hotspots) and external (variable obscuration, micro lensing) phenomena have been proposed to explain the flux variations. Early investigations of AGN variability did not show any distinct features in their power spectral density (PSD), suggesting that the PSD has a pure power-law shape (Green, McHardy & Lehto 1993; Lawrence & Papadakis 1993) that unfortunately has little power to discriminate among variability models. However several factors (analogies with Galactic binaries, dependence of the emission on the physics of the accretion process, unphysical behaviour when extrapolated to long time-scales) indicated that some characteristic time-scale should be observable in the PSD.

More recently, the combination of long observing campaigns over several decades, and shorter high-quality XMM–Newton observations, has allowed the discovery of at least one, and in some cases two, breaks in the PSD of nearby AGNs (Papadakis et al. 2002; Uttley, McHardy & Papadakis 2002; Markowitz et al. 2003; McHardy et al. 2007). Such features seem linked to both the BH mass and the properties of the accretion flow, and have enabled using variability to test the properties of the accretion flows and to measure the main physical parameters (BH mass, accretion rate) of the AGN. Similarly, we have been able to see extreme cases of variability induced by varying column densities of obscuring material, in type 1 (Yang et al. 2016) and type 2 AGNs (e.g. Risaliti et al. 2009, 2011; Giustini et al. 2011; Nardini & Risaliti 2011) and in BAL quasars (e.g. Lundgren et al. 2007; Gibson et al. 2008, 2010).

Most of our knowledge about AGN variability in the X-ray band is derived from extensive observations of nearby and mostly low-luminosity AGNs, as these were the only ones initially accessible by low effective area and/or low spatial resolution instruments. Such facilities have been the ones allowing the long and regular monitoring campaigns required to avoid the problems introduced by low statistics and irregular sampling in the temporal analysis. The extension of such studies to a larger population of distant sources requires both a large effective area and a good angular resolution to avoid crowding effects. Progress was made adopting a less sophisticated approach, measuring the integrated power over long time-scales in an attempt to investigate the variability properties of AGNs over cosmological volumes. For instance Almaini, Lawrence & Sharks (2000) and Manners, Almaini & Lawrence (2002) studied samples of quasi-stellar objects (QSOs) selected from ROSAT surveys up to z ≲ 4; Paolillo et al. (2004) analysed the variability properties of AGNs in the ChandraDeep Field-South (CDF-S) using Chandra data, Papadakis et al. (2008) and Allevato et al. (2010) used XMM–Newton observations to study the variability of AGNs in the Lockman Hole and the CDF-S, respectively, Vagnetti, Turriziani & Trevese (2011) and Middei et al. (2017) investigated serendipitous XMM–Newton, Swift and ROSAT samples, while Shemmer et al. (2014, 2017) explored a group of luminous quasars combining ROSAT, Chandra and Swift observations.

These works have shown that variability is ubiquitous in AGNs and that it has similar properties to those of nearby and less luminous AGNs, but also suggested that its amplitude may increase with lookback time and is possibly a tracer of the higher average accretion rates present in the earlier Universe. The results so far are not conclusive and suffer from biases due to sparse sampling and low statistics, as well as randomness intrinsic to red noise processes. For instance, Gibson & Brandt (2012), studying a serendipitous sample of Sloan Digital Sky Survey (SDSS) spectroscopic quasars, confirmed several results of previous works but failed to detect any clear evidence of increased variability at large redshifts. Similar conclusions were reached by other authors: Mateos et al. (2007) and Lanzuisi et al. (2014) studying the XMM–Newton light curves of AGNs in the Lockman Hole and Cosmic Evolution Survey (COSMOS) fields, and by Vagnetti et al. (2016) using the MEXSAS serendipitous sample.

The CDF-S represents the deepest observation of the Universe in X-rays. As discussed above, the first 1 Ms data were used by Paolillo et al. (2004) to investigate the nature of variability in distant AGNs, but many of their results were only marginally significant due to the low number of sources, the limited availability of spectroscopic redshifts and the limited time-scale coverage. This data set has grown over time to span, with the 7 Ms data presented in Luo et al. (2017), a time interval of ∼17 yr, reaching a depth of 1.9 × 10⁻¹⁷, 6.4 × 10⁻¹⁸ and 2.7 × 10⁻¹⁷ erg cm⁻² s⁻¹ for the 0.5–7, 0.5–2 and 2–7 keV bands, respectively, and has accumulated a wealth of ancillary multiwavelength data. This work is thus intended to test and extend the previous results, and link them to our knowledge based on nearby samples. We have already exploited these data in part in Shemmer et al. (2014) where we used the 2 Ms CDF-S light curves to compare the bulk variability of radio-quiet AGNs to bright high-redshift quasars, in Young et al. (2012) where we used the 4 Ms data to detect the faint AGN population in normal galaxies by means of variability and in Yang et al. (2016) to investigate the long-term variability of AGNs. Here we present a more refined analysis of the 7 Ms light curves probing different temporal time-scales in order to understand the connection between variability and AGN physical properties at z > 0.5.

The paper is organized as follows: in Section 2 we discuss the data and the light curves extraction process; Section 3 explains how we detect and characterize the population of variable sources; Section 4 explains how we measure the average variability and study its dependence on the properties of the AGN population; Section 5 discusses how we test different variability models and use them to constrain the AGN accretion history with lookback time. Finally, in Section 6, we discuss our results and present our main conclusions.

Throughout the paper we adopt values of H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3 and |$\Omega _\Lambda = 0.7$| (Spergel et al. 2003).

2 THE DATA

The CDF-S data used here are those described in detail in Luo et al. (2017, also see Luo et al. 2008a; Xue et al. 2011). For completeness we shortly summarize here the main properties of the data set, referring the reader to the above papers for a thorough discussion of the data properties. The data set consists of 102 observations collected by Chandra between 1999 and 2016, adding up to a total exposure time of 6.727 Ms; the individual observations have exposure times ranging from ∼9 up to 141 ks, and have very similar aimpoints within ∼1 arcmin, although different roll angles (see table 1 in Luo et al. 2017). The data reduction procedure adopted in order to create event lists and exposure maps is described in detail in Luo et al. (2017) and we refer the reader to that paper for details.

In order to study the temporal behaviour of the AGN population, we extract AGN light curves following the same procedure adopted by Paolillo et al. (2004) for the 1 Ms data set. We start from the main source catalogue of Luo et al. (2017), consisting of 1008 X-ray sources, which represent our ‘main sample’; we ignore instead the supplementary near-infrared bright catalogue in Luo et al. (2017), as these sources are all too faint to be useful in our analysis. For each source we measured counts within a circular aperture with variable radius R_S depending on the angular distance θ in arcsec, from the average aimpoint: |$R_S=2.4\times \textit{FWHM} $| arcsec, where |$\textit{FWHM}=\sum_{i=0,2} a_i\theta^i $| is the estimated full width half maximum of the point spread function (PSF) and a_i = {0.678, −0.0405, 0.0535} (see Giacconi et al. 2002, for further details). The only difference with respect to Paolillo et al. (2004) is that the minimum radius was set to 3 arcsec, in order to exploit fully the sharp Chandra PSF in the field of view (FOV) centre, and minimize the cross-contamination between nearby sources.¹ Similarly the local background for each source was measured in a circular annulus of inner (outer) radius R_S + 2 (R_S + 12) arcsec. Neighbouring objects were always removed from the source or background region when they overlapped. This approach is less sophisticated than the one adopted by Luo et al. (2017) who used the acis extract software to model the Chandra PSF and extract fluxes within polygonal regions; however, it has the advantage of being simpler and avoiding the low signal-to-noise ratio (S/N) wings of the Chandra PSF (also see Vattakunnel et al. 2012). A comparison between our total fluxes and those derived in Luo et al. (2017) shows that on average our extraction procedure recovers 95 per cent of the total source flux. In any case we stress that we use our aperture photometry only for the variability analysis, while turning to the more accurate Luo et al. (2017) photometry to obtain total fluxes and luminosities.

We binned the data into individual observations: this allows derivation of light curves with 102 points over a 17 yr interval. Sources near the edge of the detector may be missing in several observations due to the different aimpoint and roll angles of each pointing. We follow Paolillo et al. (2004) in retaining only the epochs where >90 per cent of the source region and >50 per cent of the background region fall within the FOV. In any case 884 (88 per cent) of our sources have light curves with at least 50 bins and 758 (75 per cent) are sampled by all 102 observations. The light curves were extracted both in the full 0.5–8 keV band, and in the 2–8 keV rest-frame band for the 986 sources with available redshifts.² Examples of CDF-S light curves are shown in Fig. 1. A short movie showing the variability of sources in the entire CDFS field can be found at http://people.na.infn.it/paolillo/MyWebSite/CDFS.html.

Figure 1.

Example of CDF-S light curves in the 0.5–8 keV band, for sources with different total counts and temporal behaviour; given the sparse cadence of CDF-S observations, each panel groups nearby observations separated by large temporal gaps; specifically the first three panels represent the 1 Ms data presented by Giacconi et al. (2002) and used in Paolillo et al. (2004), the fourth, fifth and sixth panels cover the additional 1, 2 and 3 Ms presented, respectively, by Luo et al. (2008a), Xue et al. (2011) and Luo et al. (2017). The average count rate is marked by a dotted line; a solid line shows the zero count rate level. The source number from the Luo et al. (2017) catalogue and the total source counts are shown in the upper left-hand corner of each row. The squares, in source 367, mark observations where part of the source background falls outside the FOV. The dark, medium and light grey bars in the top row highlight the temporal segments used to compute variability on short, intermediate and long time-scales, as discussed in Section 5.1.

Open in new tabDownload slide

3 FINDING VARIABLE AGNs

Figure 2.

Cumulative fraction of variable sources in the CDF-S, as a function of the net source counts. The number of variable sources increases at higher counts (i.e. S/N). We show both the entire main sample, and the subsample of sources with L_X > 10⁴² erg s⁻¹ that minimizes contamination by normal galaxies. Error bars represent the 95 per cent binomial uncertainty (calculated using the Bayesian approach of Cameron 2011).

Open in new tabDownload slide

We note that Paolillo et al. (2004) implemented also an additional variability estimate (⁠|$\chi ^2_{\text{max}}$|⁠) based on the maximum deviation from the mean, observed among all the bins in the light curve. This approach was useful both because in the 1 Ms data set we had at most 11 points in each light curve and thus low variability levels could be hard to detect, and also to detect short transient events. The 7 Ms data set however has much better sampled light curves, and spans longer time-scales, thus probing lower frequencies where the variability power of AGNs is expected to be larger due to the red noise PSD. Moreover Luo et al. (2008b) searched the first 2 Ms of CDF-S data, finding no evidence that short-duration events (durations of few months, such as stellar tidal disruptions) may dominate the observed variability, and only a few fast transient has been observed over the full 7 Ms data (Bauer et al. 2017; Zheng et al., in preparation). In this work we thus concentrate on variability estimates that are averaged over many epochs.

The X-ray luminosity versus redshift distribution of sources from the 7 Ms Luo et al. (2017) sample is presented in Fig. 3. The rest-frame 0.5–8 keV luminosity was calculated by Luo et al. (2017) modelling the X-ray emission using a power law with both intrinsic and Galactic absorption; the column density was constrained finding the value that best reproduced the observed hard-to-soft band ratio, assuming an intrinsic power-law photon index of Γ_int = 1.8 for AGN spectra. Variable sources (solid circles) are detected up to z ∼ 5, but lie preferentially among the brightest sources at any redshift due to the large number of counts required to detect flux changes (see also Fig. 4). As was the case also in the 1 Ms data, there are several variable sources below the L_X = 10⁴² erg s⁻¹ limit, often adopted to separate AGNs from normal galaxies; this is not surprising since many galaxies are expected to host low-luminosity AGNs whose emission significantly contributes to the overall galaxy X-ray luminosity. In fact Young et al. (2012) already searched the 4 Ms data to identify low-luminosity AGNs, using longer integration time-scales (4 × 1 Ms bins) in order to increase the likelihood of detecting variability in faint sources.

Figure 3.

X-ray luminosity versus redshift. Solid dots mark variable sources. The horizontal dotted line marks the L_X = 10⁴² erg s⁻¹ limit generally used to discriminate AGNs from galaxies, while the solid line shows the 7 Ms flux limit.

Open in new tabDownload slide

Figure 4.

Excess variance versus total source counts (top panel) and S/N per bin (bottom panel) in the 0.5–8 keV band. Variable and non-variable sources are plotted as solid and open circles, respectively. Note that the some (faint) sources lie outside the plotted range. The vertical dotted lines show the counts and S/N limits adopted in this work to ensure that the observed variability is not dominated by statistical uncertainties.

Open in new tabDownload slide

4 MEASURING THE AGN VARIABILITY AMPLITUDE

The formal error on |$\sigma ^2_{\text{NXS}}$| is given by Turner et al. (1999), assuming stationarity and uncorrelated Gaussian processes. Subsequently Vaughan et al. (2003) provided an alternative approach, more suited to compare the temporal behaviour in different energy bands. These estimates however only account for measurement errors and not for the random scatter intrinsic to any red noise process. In addition, as shown in Allevato et al. (2013), the irregular sampling pattern in the case of sparsely sampled light curves should also introduce additional scatter that will depend on the sampling scheme and the intrinsic (and a priori unknown) PSD shape.

Fig. 4 shows the excess variance of the main sample as a function of total source counts and average S/N per bin⁴ in the full Chandra energy band. At low count rates or low S/N there is a very large scatter in the excess variance and a significant fraction of sources have negative values. As the S/N increases (average S/N per bin ≳0.8, corresponding roughly to total counts ≳350, see Fig. 4) the distribution skews significantly towards positive excess variances, reflecting the improved ability to measure the intrinsic source variance. We measure a median variance |$\sigma ^2_{\text{NXS}}=0.14_{-0.08}^{+0.16}$| corresponding to count rate fluctuations of ∼40 per cent (σ_NXS = 0.37) for variable sources with >350 counts, where the uncertainties are the lower and upper quartiles. This is ∼30 per cent larger than observed by Paolillo et al. (2004, where σ_NXS = 0.28), as expected if AGNs have a red noise PSD whose power increases on the longer time-scales sampled here, and possibly also due to the fact that we are probing fainter sources⁵ that tend to be intrinsically more variable (see Section 4.1).

Motivated by the discussion in the previous paragraph, we decided to create two new samples of sources with (1) S/N per bin >0.8, in the observed (0.5–8 keV) and the rest-frame (2–8 keV) bands, and (2) more than 90 points in their light curves (to exclude sources at the edge of the FOV sampled only by part of the observations). We name them as the ‘bright-O’ and ‘bright-R’ samples, respectively. The S/N lower limit value of 0.8 is reinforced by the results of Allevato et al. (2013) who showed that such a threshold is necessary to measure accurately the excess variance in sparsely sample data, as long as we average 10–20 individual measurements. For the bright-R sample all quantities, including luminosities and excess variance, are computed in the rest-frame 2–8 keV band.

We also used the Extended Chandra Deep Field-South (E-CDF-S) radio catalogue by Bonzini et al. (2013), which classifies radio sources based on their infrared 24 μm to radio 1.4 GHz flux density ratio, to identify radio-loud AGNs in our samples (see their sec-tion 3.2). There are six and five radio-loud AGN in the bright-O and bright-R samples, respectively. Although we do not detect a significant difference in the average variability of such sources from the rest of the AGN population, we decided to remove them anyway from the subsequent analysis as the physical origin of the variability may be different (e.g. originating from the jet). The final number of sources in the bright-O and bright-R samples is 110 and 94, respectively. The different sizes of the two samples are due to the presence of sources with missing redshift and to the lower average S/N in the rest-frame band.

4.1 Dependence of AGN variability on luminosity, redshift and local density

In Fig. 5 we plot the measured |$\sigma _{\text{NXS}}^2$| against the (absorption-corrected) X-ray luminosity for sources in the bright-O and bright-R sample (grey points). The black crosses show the mean |$\sigma ^2_{\text{NXS}}$| (and luminosity) over 15 sources. This mean should be representative of the intrinsic excess variance (assuming that all sources in a given luminosity bin have a similar intrinsic |$\sigma ^2_{\text{NXS}}$|⁠). We point out that it is a common mistake to remove non-variable or negative |$\sigma ^2_{\text{NXS}}$| sources. As shown in Allevato et al. (2013), the |$\sigma ^2_{\text{NXS}}$| distribution extends to negative values especially for low variability/low S/N sources; removing these values would bias the variability ensemble estimates. We therefore used the excess variance measurements of all sources in each bin, irrespective of whether they are positive or negative, to estimate the mean |$\sigma ^2_{\text{NXS}}$|⁠. The error of the individual points in Fig. 5 takes into account the measurement error of the points in the light curves, only, and have been estimated using the equations in Turner et al. (1999), as discussed above. Instead, the error of the mean excess variance in each bin is estimated following Allevato et al. (2013), and should be representative of the true, overall uncertainty of the mean excess variance.

Figure 5.

Excess variance over the full 7 Ms (17 yr span) versus luminosity for the bright-O and bright-R samples (left- and right-hand panel, respectively). The individual sources are plotted as light grey circles of increasing size as a function of redshift. Black crosses show the average variance and its error in bins with ∼15 sources. Variable/non-variable sources are shown as filled/empty symbols, respectively; sources with excess variance ≤5 × 10⁻³ are shown at the 5 × 10⁻³ level. The S/N threshold reported in the Figure titles refers to the bins in the light curves, as explained in the text.

Open in new tabDownload slide

Fig. 5 shows that AGN variability is anticorrelated with (unabsorbed) X-ray luminosity, thus confirming the results obtained in previous investigations of the CDF-S on different time-scales by Paolillo et al. (2004), Young et al. (2012), Shemmer et al. (2014) and Yang et al. (2016). The two panels in the figure show that the dependence of the variability amplitude on X-ray luminosity is very similar between the bright-O and bright-R samples. However, in order to (1) avoid the effects due to the absorbing column and its variations in time, which mainly influence energies below 2 keV; (2) to eliminate complications in the interpretation of our results due to the differences in the PSD between the soft and hard band seen in some local AGNs (e.g. McHardy et al. 2004b) and (3) to allow a better comparison of our results with those from with variability studies of local AGNs on long time-scales, which are mainly based on RXTE data and focus on the 2–10 keV energy range, from now on, we will only use the rest-frame 2–8 keV measurements of the bright-R sample.

Fig. 6 shows the dependence of the average excess variance (including both variable and non-variable sources) on redshift. The dashed line in this figure shows the running average of the excess variance measurements of the sources in the bright-R sample and its error (estimated as explained above). We observe the variability amplitude to decrease with increasing redshift, which is consistent with the fact that at high redshift we probe higher luminosity sources (see Fig. 3) that are intrinsically less variable (Fig. 5).

Figure 6.

Running average excess variance (in bins of 20 sources) as a function of redshift; the error on the mean (68 per cent uncertainty) is shown by the grey shaded area. The solid line shows the volume density of sources, with an arbitrary normalization. The vertical dashed line indicates the position of the spatial density peak at z ≃ 0.6–0.7.

Open in new tabDownload slide

It is interesting to compare the dashed and the solid lines in Fig. 6. The solid line indicates the volume density of the CDF-S sources as a function of redshift (in arbitrary units, so that it can be easily compared with the dashed line). The CDF-S region is characterized by several overdensities in redshift space, due to the presence of large-scale structures (e.g. Gilli et al. 2003). The most prominent one, at z ≃ 0.7, also contains a large number of variable sources down to L_X ≃ 10⁴¹ erg s⁻¹. The comparison of the mean excess variance with the local volume density does not any correlation above the statistical uncertainty. The lack of any excess of the average excess variance, coincident with the z ≃ 0.7 density peak, suggests that the variability amplitude in these AGN is not affected by environmental effects that could trigger and enhance variability through, e.g. enhanced accretion processes or dynamical instabilities.

4.2 Variability dependence on time-scale

In addition to luminosity and redshift, if the intrinsic variability process has a ‘red noise’ character, the excess variance also depends on the rest-frame duration of the light curves, which usually span a fixed time interval in the observer's frame. To investigate this issue, we computed the excess variance of each object in the bright-R sample on four different time-scales (in the observer's frame): 6005, 654, 128 and 45 d (the ‘7 Ms’, ‘long’, ‘intermediate’ and ‘short’ time-scales, respectively). The first time interval is simply the total duration of the full 7 Ms data set. The ‘long’ time-scale measurement was obtained using data only from the last 3 Ms that correspond to the points covered by the light grey bar in the top row of Fig. 1, between 5350 and 6000 d. The ‘intermediate’ excess variance was measured over the interval 2–4 Ms, i.e. between 3800 and 3940 d (intermediate grey bar in Fig. 1). The excess variance of the shortest time-scales is the hardest to estimate since the variations on these time-scales are usually dominated by statistical noise. To increase the reliability of our measurement we averaged the variance measured from six different short time intervals (420–450, 2900–2950, 3800–3840, 3860–3900, 3910–3940 and 5455–5500 d, dark grey bars in Fig. 1) where the Chandra observations had a dense cadence and the sampling is more uniform. As shown by Vaughan et al. (2003) and Allevato et al. (2013) averaging over multiple observations allows to reduce the intrinsic scatter on |$\sigma ^2_{\text{NXS}}$|⁠. To reduce further the uncertainty of the excess variance estimates, we limited the analysis to objects whose lightcurves have an average S/N per bin >1.5.

Since the sampled time-scales correspond to different rest-frame time-scales at different redshifts, we grouped our sources (up to z ∼ 2) in the two redshift intervals listed in the legend of Fig. 7. They were defined in such a way that the interval width was kept as small as possible to reduce the internal difference in rest-frame time-scales, and at the same time there were at least 15 sources in each bin. The average |$\sigma ^2_{\text{NXS}}$| of all the sources in each redshift bin is shown in Fig. 7, as a function of the maximum rest-frame time-scale (estimated at the mean z of each bin).

Figure 7.

Excess variance versus maximum sampled time-scale for sources in the bright-R sample with S/N per bin >1.5, in two redshift ranges. The solid black line indicates the model excess variance in the case when PSD(ν) ∝ ν⁻¹. The discontinuous lines represent instead the prediction from our best-fitting bending power-law model (Model 4, Section 5.1).

Open in new tabDownload slide

At each time-scale, the higher redshift measurements are systematically smaller than the average variability amplitude in the lower redshift bins. This is due to the different luminosity ranges sampled at each redshift. However, the important result is that, for both redshift bins, the variability amplitude clearly decreases towards shorter rest-frame time-scales. Although the data plotted in Fig. 7 are not directly PSD measurements (since the excess variance estimates the integral of the PSD between the minimum and maximum sampled time-scale), the decrease of the excess variance with decreasing time-scale is direct observational evidence for the red noise nature of the variability process of the high-redshift AGNs.

To demonstrate that this is indeed the case, in Fig. 7 we overplot a model prediction based on the assumption that the average intrinsic PSD has a power-law shape. The black solid line shows |$\sigma _{\text{mod}}^2$| when the PSD is a single power law of the form PSD(ν) = Aν⁻¹ (this model is appropriate for local AGNs on long time-scales; e.g. Uttley et al. 2002; Markowitz et al. 2003; McHardy et al. 2004a, 2006). In this case |$\sigma _{\text{mod}}^2 = A\ln (\frac{\nu _{\text{max}}}{\nu _{\text{min}}})$|⁠, where ν_min and ν_max are the lowest and highest sampled rest-frame frequencies. We fixed |$\nu_{\max}=(1+z)/(86400\, \Delta t_{\min}^{\rm obs})\, {\rm s}^{-1}$|⁠, where |$\Delta t_{\min}^{\rm obs}=0.25\,\rm d $|⁠,⁶ using the average redshift of the low-z sample. We also chose A so that the model excess variance matches the excess variance of the longest time-scale for the low-redshift bin, in order to display the model behaviour.

Qualitatively, the model predictions (decrease of excess variance with decreasing t) are similar to what we observe. However, the model has a shallower slope than the observed one. This suggests that such a flat PSD, typical of local AGNs on long time-scales, is inadequate to describe the measured excess variance on short time-scales. A bending power-law model with a high-frequency cut-off, described in detail in Section 5.1, does a much better job in reproducing the observed trend. We believe that Fig. 7 not only demonstrates that the high-redshift AGNs have PSDs that are well represented, on average, by a power law, but also that their PSD ‘breaks’ above some characteristic frequency, as observed in several nearby AGNs.

5 TESTING VARIABILITY MODELS AND TRACING THE AGN ACCRETION HISTORY

The discussion above demonstrates that it is difficult to draw firm conclusions about the dependence of the variability amplitude on the underlying AGN physical parameters and its evolution, based on the excess variance versus luminosity/redshift/time-scale plots. However, with the CDF-S data, we can now study more accurately the |$\sigma _{\text{NXS}}^2\hbox{--}L_{\rm X}$| relation at different redshifts, by treating properly the differences in sampled luminosities and time-scales (due to differences in z). To this end we divided the bright-R sample sources in four redshift intervals and we computed the average |$\sigma ^2_{\text{NXS}}$| of sources in luminosity bins containing at least 15 sources. We considered the two redshift bins that we also considered in Section 4.2, and to increase the redshift range we also considered the [1.8–2.75] and [2.75–4] redshift bins.⁷

The four columns of Fig. 8 show the |$\sigma ^2_{\text{NXS}}$| measurements in the four different redshift intervals plotted as a function of X-ray luminosity, for the four different time-scales discussed in Section 4.2. The decrease in variability with increasing X-ray luminosity is confirmed on most time-scales and redshift bins, at least up to z ∼ 2 where we probe a large enough range of luminosities. Arguably, the uncertainty on the individual points is large, but we do not observe any significant increase of variability amplitude with redshift, in any of the time-scales we considered. The amplitude of the |$\sigma ^2_{\text{NXS}}\hbox{--}L_{\text{X}}$| relations increases with increasing time-scale, which is caused by the red noise character of the observed variations. On the shortest time-scales, the anticorrelation between |$\sigma ^2_{\text{NXS}}$| and L_X is steep. Then it flattens as we sample increasingly longer time intervals; this behaviour is in agreement with the scenario where the intrinsic PSD is represented by a bending power law with a high-frequency cut-off.

Figure 8.

The |$\sigma ^2_{\text{NXS}}$| versus L_X relation for the bright-R sources in different redshift intervals. The four columns correspond to the different time-scales discussed in Section 4.2. The CDF-S data are binned in groups of >15 sources (the exact number is printed above each point). The arrows in the panels of the first column represent the 3σ upper limit. The lines in each panel show the predictions of Models 1, 2, 3 and 4 with parameters set to the best-fitting solutions reported in Table 1 (which allow a variable Eddington ratio as a function of redshift). The open triangles and solid cyan line in the rightmost panels represent the local AGN data from Zhang (2011) (each point includes 14 sources) and the corresponding model prediction.

Open in new tabDownload slide

In order to understand the observed complex dependence of the variability–luminosity relation in various redshift bins and over different time-scales, we fitted the data shown in Fig. 8 with predictions of PSD models, which are frequently used to parametrize the observed power spectra of nearby, X-ray bright AGN, as we explain below. To constrain better the models at the lowest redshifts, which are not sampled by the CDF-S population, we also considered the data from the sample of local AGNs studied by Zhang (2011). The Zhang (2011) light curves are based on 14 yr long RXTE monitoring campaigns, closely matching the full 7 Ms observed time-scales. Note that the RXTE monitoring cadence only allows us to probe the longest time-scales, with no equivalent to the additional long, intermediate and short time-scales probed for the CDF-S sources. Therefore, the Zhang (2011) data are thus only shown in the rightmost panels of Fig. 8.

5.1 Modelling the AGN variability

To link the variability of the AGN to its physical properties we explore four different variations of the PSD model defined by equation (2).

• Model 1. The PSD amplitude at ν_b is constant, in terms of ν_b × PSD(ν_b) = 0.02, for all AGNs as suggested by Papadakis (2004, see also González-Martín et al. 2011), while the break frequency depends on the BH mass as ν_b = 580/(M_BH/M_⊙) s⁻¹ as suggested by Gonzalez-Martin & Vaughan (2012).

• Model 2. The PSD amplitude is constant as in Model 1, while the break frequency depends both on the BH mass and the accretion rate so that

  \begin{equation*} \nu _{\rm b}=(200/86\,400) (L_{44, {\rm bol}}) (M_{6, {\rm BH}})^{-2}\, \mbox{s}^{-1} \end{equation*}

  as in the prescription of McHardy et al. (2006), where L_{44, bol} is the bolometric luminosity in units of 10⁴⁴ erg s⁻¹ calculated applying the recipe of Lusso et al. (2012),⁸ and M_{6, BH} is the BH mass in units of 10⁶ M_⊙.
• Model 3. Similar to Model 1, where ν_b = 580/ (M_BH/M_⊙) s⁻¹, but the PSD normalization is itself dependent on the accretion rate as

  \begin{equation*} \nu _{\rm b}\times \mbox{PSD}(\nu _{\rm b})=3\times 10^{-3}\lambda _{{\rm Edd}}^{-0.8} \end{equation*}

  (with |$\lambda _{{\rm Edd}}=\frac{\dot{m}}{\dot{m}_{{\rm Edd}}}$|⁠) as proposed by Ponti et al. (2012).

• Model 4. Similar to Model 2, where ν_b = (200/86 400) (L_{44, bol}) (M_{6, BH})⁻² s⁻¹, but the PSD normalization depends on the accretion rate as in Model 3.

5.2 Model fit procedure

We fit each one of the four models presented in Section 5.1 to the data points plotted in Fig. 8, over all luminosities, time-scales and redshifts simultaneously. Equation (3) was integrated over the range of rest-frame frequencies sampled in each redshift interval, as described in Section 5.1, in order to derive the excess variance |$\sigma ^2_{{\rm mod}}$|⁠. The observed |$\sigma ^2_{{\rm NXS}}$| is an estimator of the intrinsic light-curve variance |$\sigma _{{\rm mod}}^2$| provided that we take into account the biases introduced by the sampling pattern (red noise leakage and uneven sampling). We adopted the Allevato et al. (2013) recipe to correct for such biases, deriving the predicted excess variance as |$\sigma _{{\rm pred}}^2= \sigma ^2_{{\rm mod}} / (C \times 0.48^{\beta -1})$|⁠, where β is the PSD slope below the minimum sampled frequency ν_min, and C is a corrective factor dependent on the sampling pattern. In our models β was estimated as the average slope from ν = ν_min/5 and ν_min, i.e. the range from where most of the leaking power originates (cf. Allevato et al. 2013). The factor C ranges from 1 to 1.3 and allows us to account for the missing power due to the gaps in the light curve; we adopted the upper value of C = 1.3 adequate for the sparse sampling of CDF-S light curves and when β = 1.0–1.5, but we verified that changing the bias factor yields consistent results within Δλ_Edd ∼ 0.03, where higher accretion rates correspond to lower bias values.

The only free parameter in our models is the average accretion rate λ_Edd. For a given accretion rate we compute the bolometric luminosity L_bol, using X-ray luminosity value of each point, according to the Lusso et al. (2012) prescription. Using λ_Edd and L_bol, we then compute the average BH mass for all AGNs that contributed to |$\sigma _{{\rm NXS}}^2$| in each bin, and ν_b and the normalization A according to each model prescription. Knowing ν_b and A is then possible to compute |$\sigma _{{\rm mod}}^2$| using equation (3) for each point in the |$\sigma _{{\rm NXS}}^2\hbox{--}L_{\rm X}$| relations. The best-fitting λ_Edd is then found by χ² minimization of the differences between |$\sigma _{{\rm NXS}}^2$| and |$\sigma _{{\rm pred}}^2$| (we note again that we used |$\sigma _{{\rm pred}}^2$| and not |$\sigma _{{\rm mod}}^2$|⁠, to take into account the bias due to sampling and red noise leak).

We note that the CDF-S data used in the fit are not entirely independent. In fact, as discussed above, |$\sigma ^2_{{\rm NXS}}$| measures the integral of the PSD that, for the long time-scale bins, includes some contribution coming from the same high frequencies sampled on short time-scales. However this correlation is not expected to be strong, since the PSD has a steep negative slope so that, on each time-scale bin, most of the power comes from frequencies close to the minimum sampled value. In any case the correlation would reduce the degrees of freedom (d.o.f.) of our data thus yielding lower probabilities and strengthening our conclusions about which models should be rejected.

5.3 Model fit results

Initially we fitted the data allowing a different λ_Edd for each redshift bin. The best-fitting models are plotted in Fig. 8 and the best-fitting results are listed in Table 1, together with the minimum χ² and the likelihood of the model. For completeness, we report the best-fitting results in the case when we consider the CDF-S data only, and in the case when we add the Zhang (2011) data as well. The best-fitting λ_Edd values for the CDF-S sources are identical (as each redshift bin is independent from the others), however the addition of Zhang's low-redshift data constrains the models better (in terms of the null hypothesis probability). For that reason we discuss the best-fitting results when we fit both the CDF-S and the low-redshift data.

From Fig. 8 we see that the models become steeper on the shortest time-scales, due to the presence of the PSD break, and flatten on the longest one where most of the power comes from the ν⁻¹ part of the PSD. The Zhang (2011) light curves, on the other hand, are not affected by time dilation due to their redshift z ≃ 0, so that the variability is not expected to anticorrelate with luminosity, since the PSD break falls outside the range of probed frequencies.

The best-fitting results show that Model 3 is formally rejected at >99 per cent confidence level. Model 1 provides a statistically acceptable fit, but with an extremely low Eddington ratio, at all redshift bins except the [0.4–1.03] bin where it is unconstrained. In fact, the resulting best-fitting λ_Edd values are so low that even the low-luminosity CDF-S AGN should have BH masses larger than ∼10⁸ M_⊙ at all redshifts to explain their observed luminosity. When we force the model to have any value of λ_Edd > 0.03 results in the rejection of the model at the >99 per cent level. We therefore conclude that models where the break frequency depends on BH mass only (irrespective of whether the PSD amplitude depends on Eddington ratio or not) are not consistent with the data. Model 2 is also formally rejected and at the >99 per cent confidence level. We note that Models 1, 2 and 3 all show some tension with the CDF-S data at the lowest luminosities on long time-scales (see Fig. 8), as their normalizations are too low. On the other hand, Model 4 reproduces rather well the overall trends and dependence of variability on luminosity and time-scale, for the CDF-S data and its variable PSD normalization allows a better agreement with the observational measurements.

The behaviour of λ_Edd as a function of redshift is presented in Fig. 9 (we do not plot the Model 1 best-fitting results, as they suggest an accretion rate that is either unconstrained or extremely low). Model 3 predicts an increase of the accretion rate from ∼0.07 in the local Universe up to almost 0.3 at redshifts higher than 3. However the error of the best-fitting parameters are so large that we cannot claim a significant indication of an increasing accretion rate with increasing z. In fact, even in the case of Model 4 (which provides the best fit to the data) the best-fitting λ_Edd errors are so large that we cannot argue for a significant variation of the accretion rate with redshift.

Figure 9.

Eddington ratio estimates as a function of redshift from model fitting to the CDF-S and local AGNs variability (Table 1). The solid horizontal lines and shaded areas represent the results of the fit with a constant λ_Edd reported in Table 2.

Open in new tabDownload slide

For that reason, we repeated the fits (to both CDF-S and Zhang 2011 data), keeping the accretion rate fixed to a common value at all redshift bins. The best-fitting results are listed in Table 2. The results are consistent with those presented above. The horizontal solid line in Fig. 9 indicates the best-fitting λ_Edd, which is very similar to the mean of the accretion rate values listed in Table 1. Model 1 fits the data well, but with a very small accretion rate (as before), while Model 3 best fit is still rejected with a high confidence. This time, Model 2 is marginally accepted (at the 1 per cent level), but it is still Model 4 that provides again the best fit. In fact, using the F-test, we can see that the improvement of the Model 4 best fit when we let λ_Edd free is not significant when compared with the best fit in the case when λ_Edd is kept constant.

Summarizing, our analysis supports the view that the variability amplitude of the high-redshift AGN can be explained if we assume a power spectrum that is identical to the PSD of local AGNs (i.e. the variability mechanism is the same in local and high-z AGN). In particular, the variability amplitude of the CDF-S AGN, and its dependence on luminosity, z and time-scale, can be explained if the PSD break frequency ν_b depends on both BH mass and λ_Edd, in agreement with McHardy et al. (2006). Most probably, the PSD amplitude also depends on λ_Edd as proposed by Ponti et al. (2012). The Eddington ratio of the AGN population seems consistent with a constant value, at all redshifts. The quality of our data cannot allow us to detect a dependence of λ_Edd on z.

6 DISCUSSION AND CONCLUSIONS

6.1 Variability properties of high-redshift AGNs

In this work we have analysed the light curves of AGNs in the CDF-S region, using a data set spanning ∼17 yr. The observing strategy of the CDF-S survey allows us to derive light curves with similar sampling for all sources, thus minimizing the scatter introduced in timing analysis when different sampling patterns are used for different sources (Allevato et al. 2013). In order to assess the level of variability of our sources we used two different approaches, the first one based on Monte Carlo simulations suited to assess whether a source is variable within a certain confidence level, and the second one based on excess variance analysis in order to measure the intrinsic average variability of the AGN population and link it to the physical properties of the AGNs themselves.

We confirm results based on previous studies that virtually all AGNs are variable, and only the data quality prevents us from detecting variability in faint sources: 90 per cent of the sources with >1000 net counts (e.g. ∼20 counts epoch⁻¹) are detected as significantly variable at the 95 per cent confidence level. This result is due in large part to the long time-scales probed in the CDF-S data set, as the likelihood of detecting a source as variable increases with the sampled rest-frame time-scale, as expected for sources with a red noise PSD. In some local AGNs a low-frequency break has been observed in the PSD, below which the PSD shape becomes flat and the variability power becomes approximately constant. The fact that, on average, the variability of our sources is still increasing on the longest time-scales probed by our observations of ∼17 yr, constrains the position of this low-frequency break to be, on average over the sampled AGN population, at even longer time-scales, in agreement with the results of Zhang (2011) and Middei et al. (2017).

One of the most evident clues that AGN variability is related to the physical properties of the central BH is the discovery that variability anticorrelates with intrinsic AGN X-ray luminosity (and possibly also ultraviolet/optical luminosity; e.g. Collier & Peterson 2001; Kelly, Bechtold & Siemiginowska 2009; MacLeod et al. 2010; Simm et al. 2016). In X-rays this effect has been observed in samples of nearby AGNs and has been interpreted as the consequence of a larger BH mass in more luminous objects, which would also increase the size of the last stable orbit and thus smear the overall variability produced in the innermost parts of the accretion disc (Papadakis 2004). A clue that such a correlation holds also in higher redshift AGNs cames from ROSAT observations (Almaini et al. 2000; Manners et al. 2002), and has been confirmed by observations of the CDF-S (Paolillo et al. 2004; Young et al. 2012; Yang et al. 2016), of the Lockman Hole (Papadakis et al. 2008), of the COSMOS field (Lanzuisi et al. 2014) and by serendipitous XMM–Newton/Swift samples (Vagnetti et al. 2011, 2016).

In order to study the dependence of variability on the physical properties of the AGNs, we adopted an ‘ensemble’ analysis approach. Based on the results of the simulation performed by Allevato et al. (2013) we focused on a subsample of relatively bright (S/N per bin >0.8 or ≳350 counts) sources, in order to avoid introducing statistical biases in our analysis. We confirm the presence of an anticorrelation between variability and luminosity, but we fail to detect a (statistically significant) high-luminosity upturn suggested in some previous investigations of both the CDF-S sources (Paolillo et al. 2004) and of other samples (Papadakis et al. 2008; Vagnetti et al. 2011). Our analysis supports the view that high-redshift AGNs have a similar PSD to local sources, where the variability amplitude increases towards longer time-scales. Furthermore, while globally the PSD could be represented by a simple power-law model, it seems to be better reproduced by a bending power-law model, similar to many local AGNs where a high-frequency break ν_b is detected in the PSD. This is possibly the first direct indication that we can extend to high redshift the PSD behaviour observed for local sources.

Given the complex interplay among variability, luminosity, redshift and time-scale, we showed that a proper study of the behaviour of X-ray variability of high-redshift AGNs, its possible evolution and its dependence on the AGN physical parameters (i.e. BH mass and accretion rate), must account for all these dependencies simultaneously. To this end we compared the average excess variance, over three orders of magnitude in X-ray luminosity and over four different time-scales, further dividing the sample in four redshift bins, with the predictions from various PSD models. We accounted for both statistical uncertainties and systematic effects due to e.g. the different time-scales and energy ranges probed at each redshift, as well as from the sparse sampling pattern and red noise leakage, following the recipe of Allevato et al. (2013). To better constrain the variability at low redshift, where the CDF-S survey covers a volume too small to detect a significant number of AGNs (i.e. z ≲ 0.02), we further included in the analysis the sample of local AGN studied by Zhang (2011), whose light curves are comparable to those of the CDF-S sources on the longest time-scales. All models assume a bending power law where the break frequency ν_b depends on either the BH mass or both the BH mass and the Eddington ratio λ_Edd of the sources (McHardy et al. 2006). Additionally we tried models with either a fixed PSD normalization (Papadakis 2004; González-Martín et al. 2011) or a variable one depending on the Eddington ratio λ_Edd (Ponti et al. 2012).

Our results indicate that the variability of high-redshift AGNs is consistent, within the current uncertainties, with PSD represented by a bending power law where the break frequency ν_b depends on both the BH mass and the Eddington ratio λ_Edd, in agreement with the results found by McHardy et al. (2006) for local sources. Our best-fitting model suggests that the PSD amplitude should depend on the accretion rate as well, as proposed by Ponti et al. (2012); the tension between models with a fixed PSD amplitude and the CDF-S data on the longest time-scales, observed also by Young et al. (2012) and by Yang et al. (2016), seems to support this conclusion. The data rule out the possibility that ν_b depends only on the BH mass (Gonzalez-Martin & Vaughan 2012), irrespective of whether the PSD normalization is constant or not, unless we adopt implausibly low, average accretion rates. We note that understanding the dependence of the PSD on the physical parameters characterizing the AGN population is crucial if we intend to use AGN variability as a cosmological probe, as proposed for instance by La Franca et al. (2014), since if the measured variability has a dependence on the accretion rate this will introduce an additional scatter in the |$M_{{\rm BH}}\hbox{--}\sigma ^2_{{\rm NXS}}$| relation used to measure cosmological parameters, as well as a systematic redshift-dependent bias if the accretion rate changes with lookback time.

6.2 Constrains on the BH accretion history

We used the model fitting results to measure the average accretion history of AGNs. To probe a possible change of the Eddington ratio with lookback time we first allowed λ_Edd to vary, finding that models with fixed PSD normalization show no λ_Edd dependence on redshift, while models with a variable normalization suggest a possible increase up to z ∼ 2–3. Repeating the fits assuming a constant value for λ_Edd, we obtain that a λ_Edd ≃ 0.05–0.10 is consistent with the data for all the models except one.

Given the large uncertainties, an increase of the Eddington ratio with lookback time is not ruled out by our data, but in any case the increase cannot be as strong as suggested by previous studies of AGN variability in X-ray-selected samples, where it seemed that the high-redshift population at z > 2 could be dominated by near-Eddington accretors (e.g. Almaini et al. 2000; Manners et al. 2002; Paolillo et al. 2004; Papadakis et al. 2008). This difference is due in part to the improved statistical approach based on the simulation of Allevato et al. (2013), but also from a new and more physical modelling of AGN variability on different time-scales, luminosities and redshift simultaneously, allowed by the E-CDF-S data. This suggests again that extreme care should be used in interpreting variability results based on small, sparse samples without taking into account all sources of bias. Arguably, the large uncertainties due to the limited sample size do not yet enable drawing definitive conclusions on the evolution of accretion with redshift but we point out that we are quantitatively constraining λ_Edd(z) through X-ray variability measurements for the first time.

The average Eddington ratios derived for the CDF-S sample are in good agreement with estimates in the literature. For instance Lusso et al. (2012), probing AGNs in the XMM–COSMOS survey, find λ_Edd in the range 0.015–0.26, dependent on both bolometric luminosity and AGN type. Their results are in even better agreement with ours considering that our median log (L_bol) ≃ 45.5 and for such luminosity the COSMOS sample yields a λ_Edd = 0.07–0.12 inclusive of redshift and AGN-type dependence. On the other hand, Lusso et al. (2012) did not find any evidence of λ_Edd evolution with redshift, although they probe AGNs with z ≲ 2.3 and are thus less sensitive to the higher redshift population than us. Brightman et al. (2013), exploring the COSMOS and E-CDF-S, find Eddington ratios spanning a similar range of Lusso et al., but skewed towards a slightly higher median λ_Edd ∼ 0.15; this is likely a consequence of the somewhat higher median X-ray luminosity probed by their sample compared to this work, since their hard 2–10 keV log (L_X) ∼ 44.2 erg s⁻¹, while our median value is log (L_X) = 43.6 erg s⁻¹. Similarly Suh et al. (2015) extend this type of study to a sample including the Lockmann Hole, finding a broad Eddington ratio distribution described by a lognormal distribution peaking at λ_Edd ∼ 0.25. Again however this sample extends to higher log (L_bol) ∼ 47 than sampled by our work. Studies of bright quasars, such as the SDSS-Data Release 7 (DR7) sample of Shen et al. (2011), have also usually reported higher average Eddington ratios λ_Edd > 0.1 with a sizeable fraction of sources accreting close to the Eddington rate (see e.g. fig. 4 in Wu et al. 2015); however, quasar samples are characterized by bolometric luminosities that are between one and two orders of magnitude greater than all the studies mentioned above, including our own. In the X-ray band luminous quasars have been monitored by Shemmer et al. (2014), finding that their variability is generally larger than expected from their luminosities, based on an extrapolation from lower luminosity sources in the 2 Ms observations of the CDF-S. While this may result from higher average Eddington ratios, it must be noted that unfortunately monitoring observations of such sources in the X-ray band are sparse and heterogeneous, performed through different observatories and in non-dedicated campaigns. This results in large uncertainties that prevented them from drawing definitive conclusions. We further point out that if PSD amplitude depends on λ_Edd, as suggested by our data, the conclusions of Shemmer et al. (2014) (which used fixed PSD amplitude models) may need to be revised in favour of a lower Eddington ratio, although in at least two of their cases the Eddington ratio has been estimated, independently, from H β and the values are quite high. In fact, follow up observations of high-redshift QSO with Chandra are yielding results consistent with our CDF-S observations (Shemmer et al. 2017).

Our variability-based estimates allow us to derive average λ_Edd; however, it is clear that AGNs possess an intrinsic Eddington ratio distribution. For instance the work by Lusso et al. (2012) derives a lognormal λ_Edd distribution after correcting for selection effects, while other authors (Aird et al. 2012, 2013; Bongiorno et al. 2012) claim that the intrinsic λ_Edd distribution in galaxies can be represented by a power-law function, which is independent of the host galaxy stellar mass, possibly with a soft cut-off around the Eddington limit. In any case the observed λ_Edd is influenced by selection effects, and tends to display a lognormal distribution peaking at λ_Edd ∼ 0.01–0.1. A recent revision of the work by Aird and collaborators (Aird, Coil & Georgakakis 2017) suggests that the λ_Edd distribution is more complex than a simple power law, peaking somewhere in the range λ_Edd ∼ 0.01–0.1, although still ∼Eddington limited, and dependent on host galaxy type (e.g. star forming versus quiescent). The same study also suggests that the average AGN activity shifts towards higher λ_Edd with redshift, supporting the tentative trend observed in some of our models. All these results are based on X-ray-selected samples of AGNs, but Vito et al. (2016) recently showed that BH accretion in individually X-ray undetected galaxies is negligible compared to the accretion density measured in X-ray sources. Attempting to constrain the λ_Edd distribution through variability is difficult due to the large scatter of variability measurements, which is intrinsic to stochastic processes; furthermore, at present our analysis is strongly limited by the available statistics, as we have to divide our data into luminosity and redshift bins. This also limits our ability to explore the possible dependence of λ_Edd on other parameters such as the AGN type and host galaxy. Finally, BH spin is the one physical parameter missing from any of the analyses discussed here, so that we are assuming an underlying spin distribution independent of luminosity, redshift, host galaxy type, etc. Although this is a simplistic approach, unfortunately the study of BH spins is beyond the reach of current facilities, except for a few local AGNs.

We conclude pointing out that the limits of this study are mainly due to the sample size, which prevents us from reducing the statistical uncertainties. Future wide-field/large-effective-area facilities will enable making this method competitive with other tracers, allowing one to probe more effectively the luminosity, redshift and time-scale dependence of the intrinsic AGN variability, and also to assess whether the average accretion (and possibly mass) may differ between e.g. different galaxy populations at each redshift.

Acknowledgements

We thank Fausto Vagnetti and the anonymous referee for the helpful discussions and comments that allowed to improve the manuscript. This work was supported in part by PRIN-INAF 2014 ‘Fornax Cluster Imaging and Spectroscopic Deep Survey’. BL acknowledges support from the National Natural Science Foundation of China grant 11673010 and the National Key Program for Science and Technology Research and Development grant 2016YFA0400700. YQX and XCZ acknowledge the support from the National Thousand Young Talents program, the 973 Program (2015CB857004), NSFC-11473026, NSFC-11421303, the CAS Strategic Priority Research Program (XDB09000000), the Fundamental Research Funds for the Central Universities and the CAS Frontier Science Key Research Program (QYZDJ-SSW-SLH006).

REFERENCES