© The Authors 2024

1. Introduction

In the last decade, there have been many monitoring, multi-wavelength campaigns with the goal of studying the X-ray/UV/optical variations in active galactic nuclei (AGN). These campaigns are split, broadly speaking, into two categories. The first category includes the intensive monitoring of nearby AGN using the Neil Gehrels Swift Observatory, AstroSat (to a lesser extent), and in many cases ground based telescopes as well (e.g. McHardy et al. 2014; Edelson et al. 2015, 2017; Fausnaugh et al. 2016; Cackett et al. 2018, 2020; Chelouche et al. 2019; Edelson et al. 2019; Hernández Santisteban et al. 2020; Vincentelli et al. 2021; Kara et al. 2021, 2023; Donnan et al. 2023; Miller et al. 2023; Kumari et al. 2023, 2024). Only a few, bright Seyferts belong to this category of intensively sampled observations over a prolonged period of time, albeit their number is increasing. The second category includes long-term monitoring campaigns of quasars, that is AGN which are at higher redshifts and are brighter than their nearby counterparts (e.g. Jiang et al. 2017; Mudd et al. 2018; Homayouni et al. 2019; Yu et al. 2020; Guo et al. 2022; Jha et al. 2022; Homayouni et al. 2022). There are usually many AGN observed during each campaign, but their light curves are not as densely sampled as the ones from the monitoring of the nearby Seyfert galaxies.

In all cases, the UV/optical variations are well correlated, with the variations at longer wavelengths being delayed with respect to the variations at shorter wavelengths. These results are consistent with the hypothesis of X-ray reverberation in AGN. According to this hypothesis, some of the photons emitted by the X-ray corona are directed straight to the observer, while the rest are directed towards the accretion disc. A fraction of the latter is ‘reflected’ by the disc in the X-ray range. The remaining fraction is absorbed, increasing the temperature of the disc and the emitted UV/optical flux. Since X-rays are highly variable, the additional disc emission in the UV/optical bands is expected to be also variable, revealing a time lag (with respect to the X-ray), which should increase with wavelength (as observed). This phenomenon is referred to as the disc X-ray thermal reverberation (e.g. Kazanas & Nayakshin 2001; Cackett et al. 2007).

In the last few years, we have studied in detail X-ray reverberation in AGN, for instance Kammoun et al. (2021a, hereafterK21a) and Dovčiak et al. (2022). We have been able to successfully fit the time-lags spectra of nearby AGN (Kammoun et al. 2019, 2021b, hereafter K21b; Kammoun et al. 2023), the mean and the variable, optical/UV/X-ray spectral energy distribution in NGC 5548 (Dovčiak et al. 2022; Kammoun et al. 2024), and the radial flux profiles in micro-lensed quasars (Papadakis et al. 2022). Additionally, Panagiotou et al. (2020, 2022) utilised the same model to compute model power spectral density function (PSD) in the UV/optical range, and successfully fitted the optical/UV PSDs from the extensive monitoring campaign of NGC 5548.

So far, we have studied results from the intensive monitoring campaigns of nearby Seyferts. In this work, we use archival time-lag measurements of high(er) redshift and more luminous AGN, compute their mean time-lags spectrum, and we fit them with the analytical models of K21a and b, considering the modifications recently proposed by Kammoun et al. (2023). Our objective is to investigate whether X-ray reverberation can explain the observed time lags in these objects as well.

Our paper is organised as follows: In Sect. 2, we present the archival time-lag measurements that we use in our sample. In Sect. 3 we present the theoretical models we used, while in Sect. 4 we present the results of our analysis. Finally, we summarise and discuss the results and their implications in Sect. 5.

2. Archival time-lag measurements

2.1. The data

We considered the time-lag measurements of Jiang et al. (2017), Homayouni et al. (2019), Guo et al. (2022) and Jha et al. (2022). Below we provide a short description of their data.

Homayouni et al. (2019) studied quasar light curves from the Sloan Digital Sky Reverberation Mapping Project (SDSS-RM). They combined data from the SDSS-RM project and two ground based telescopes to construct light curves which were ∼200 days long in the g and i bands. They computed time lags between the two band light curves (with the g-band as the reference) using both the interpolated cross-correlation function (ICCF; Gaskell & Sparke 1986; Gaskell & Peterson 1987; White & Peterson 1994; Peterson 2004) and JAVELIN¹ (Zu et al. 2011, 2016). We used their ICCF time-lag measurements for their ‘well defined quasar sample’ (95 objects), as listed in their Table 1. Redshift, monochromatic luminosity at 3000 Å , L_3000 Å, and black hole (BH) mass, M_BH, data were also taken from the same table. In order to compute the bolometric luminosity, L_bol, we multiplied L_3000 Å by a factor of 5.15 (Richards et al. 2006). The time-lag measurements were not redshift corrected, so we multiplied them by $\frac{1}{1+z}$ (where z is the redshift of each source), in order to bring them to the rest frame of each source.

Guo et al. (2022) studied quasar light curves from the Zwicky Transient Facility survey (ZTF) in the g, r and i bands, which were ∼1200 days long. They defined a core sample of 38 quasars, and we considered the ICCF time-lag measurements listed in their Table 3 for those targets. The time lags are given in the rest frame of each source. Redshift, monochromatic luminosity at 5100 Å, L_5100 Å, and M_BH are also listed in the same table. In order to compute L_bol, we multiplied L_5100 Å by a factor of 9.26 (Richards et al. 2006).

Jha et al. (2022) studied quasar light curves obtained from the ZTF survey in the g, r and i bands. They used both the JAVELIN and ICCF methods to compute the time lags between the g-band and the r and i-band light curves. We considered the g-r and g-i ICCF time lags listed in their Table 3. There are 19 AGN with time-lag measurements in these bands, given in the rest frame of each source. Redshift, monochromatic luminosity at 5100 Å, L_5100 Å, and BH mass for these objects were taken from their Table 1. As before, we multiplied L_5100 Å by a factor of 9.26 in order to compute L_bol.

Finally, we also considered Jiang et al. (2017), who studied quasar light curves in the Medium Deep Fields (MDFs) of the Pan-STARRS1 survey (PS1), which were ∼1200 days long. They used JAVELIN to compute the time lags between the g-band and the r, i and z-band light curves. In this work we considered their g-r, g-i and g-z time lags, as listed in their Table 3, which were given in the rest frame of each source. There are 39 quasars with time-lag measurements in all bands. Redshift, bolometric luminosity, L_bol, and M_BH, for these quasars were taken from Table 1 in Jiang et al. (2017). We note that, mixing of time-lag measurements from ICCF and JAVELIN could introduce a systematic error into the analysis, as these are fundamentally different methodologies. As we already mentioned, we chose to use time lags computed with the use of the ICCF method whenever possible, but Jiang et al. (2017) do not use the ICCF method to compute time lags. Homayouni et al. (2019), Guo et al. (2022), and Jha et al. (2022) found that time lags computed by the two methods are broadly similar (see their Figs. 7, 3 and 4, respectively). This comparison implies that use of the Jiang et al. (2017) time lags in our study will not bias our results significantly.

In total, our initial sample consists of 326 time-lag measurements for 191 AGNe. All time lags are measured between the r-, i- and z-band light curves and the g-band light curve, which is used as the reference band for all quasars in the sample.

2.2. The observed time lags

The top-left panel of Fig. 1 shows the rest frame wavelength of the g-band light curves, λ_ref, plotted versus the rest frame wavelength of the secondary bands, λ_S. The sources have various redshifts, so (λ_ref, λ_S) are not just a single point in this plot. Instead, three main diagonals appear in this figure, which correspond to the g − r, g − i and g − z rest frame wavelengths. There seems to be some scatter in the g − r and g − i diagonals, due to the fact that the effective wavelength of the respective filters in the PS1, SDSS-RM and the ZTF surveys does not appear to be the same in all of them. For the Jiang et al. (2017) data we used the g, r, i and z-band effective wavelengths listed in Tonry et al. (2012). For the Homayouni et al. (2019) data we used the effective wavelengths in Fukugita et al. (1996). Guo et al. (2022) and Jha et al. (2022) used light curves from the ZTF survey and we assumed the effective wavelengths provided by the Spanish Virtual Observatory (SVO) Filter Profile Service.

Fig. 1. Upper left panel: The rest frame wavelength of the reference band (in our work always the g-band), λref, versus the rest frame wavelength of the secondary light curves, λS, used to measure the time lags. Upper right panel: The MBH as a function of λS for the objects in our final sample. Bottom panels: Same as the top right panel but for λEdd and L2 − 10 keV (left and right panels, respectively). The boxes in all panels indicate the measurements we used to compute the mean time-lags spectrum.

The boxes in the top-left panel of Fig. 1 indicate the data we use in our study. The data in the bottom left boxes correspond to sources with similar redshift. As a result, the spread of the rest frame λ_ref and λ_S is small for these objects. The top right boxes contain roughly the same number of time-lag measurements. Our objective is to calculate the average time lag in each box, and then compare the resulting mean time-lags spectrum with various time-lags models.

The original sample of 191 AGN includes sources with a large range in M_BH and L_bol. As we explain in the next section, we use the analytical time-lags models of K21b, which were developed for AGN with M_BH ≤ 10⁸ M_⊙ and accretion rate smaller than 0.5. For that reason, we kept in the sample sources with M_BH <  10⁹ M_⊙ and λ_Edd ²  ≤  1. If we accept lower limits on M_BH and λ_Edd then the number of time-lag measurements will decrease significantly together with our ability to constrain models. In any case, we repeated the same analysis to smaller samples, with (M_BH <  5 × 10⁸ M_⊙, λ_Edd <  SPSDOUBLEDOLLAR >  < 0.5), as well as with (M_BH <  2 × 10⁸ M_⊙, λ_Edd <  0.5). In both cases, the results are almost identical, although with larger errors.

Due to the reasons we discussed above, we finally kept 128 objects in the final sample. Figure 2 shows the distribution of their z, M_BH, L_bol and λ_Edd. Redshifts range from ∼0.05 to ∼1.1. Most of the quasars in the sample have redshifts smaller than ∼0.45, but there is also a large number of objects with redshift larger than ∼0.8. This bi-modal distribution is due to selection effects. We note that our sample is not a complete sample of quasars. It simply consists of sources with measured continuum time lags in the UV/optical bands. Nevertheless, their properties are representative of quasars in larger samples. For example, the upper right plot in Fig. 2 shows the distribution for the BH mass, and the bottom left plot shows the distribution for L_bol. The BH mass is larger than 10⁸ M_⊙, and L_bol is larger than 10⁴⁵ ergs/s for most sources. The range of the M_BH and L_bol values for the sources in our sample is within the range of the respective distributions in other quasar samples like, for example the SDSS quasars in Stripe 82 (see e.g. Petrecca et al. 2024). We are missing quasars with BH mass larger than ∼10⁹ M_⊙ and L_bol >  10⁴⁶ erg/s, as the optical/UV continuum time lags have not been measured yet in such objects. On the other hand, the sources we study are brighter, and host BHs which are more massive than the BHs in nearby AGN (see e.g. Koss et al. 2022). The bottom right plot in the same figure shows the distribution for the λ_Edd, which ranges from ∼1.5 × 10⁻³ up to 1. The accretion rate is smaller than ∼0.4 for the majority of the sources in the sample, which is comparable to the accretion rate of AGN in the nearby Universe (Koss et al. 2022).

Fig. 2. The redshift, MBH, Lbol and λEdd distributions for the sources in the final sample.

2.3. The observed, mean time-lags spectrum

The grey open circles in Fig. 3 indicate all time lags for the AGN in the final sample (231 measurements in total). The observed time-lags spectrum is very scattered, mainly because each point in this figure corresponds to the time lag between various wavelengths, for quasars with different M_BH, λ_Edd, and X-ray luminosity. This is shown in Fig. 1. The upper right and bottom panels in this figure show the BH mass, accretion rate and X-ray luminosity of the sources in the sample, plotted as functions of λ_S. Data appear in roughly six groups in these plots, which correspond to the six boxes we defined in the top-left panel in the same figure.

The max-to-min ratio of M_BH and λ_Edd is larger than ∼20 in some cases. Both of these parameters affect the magnitude of the time lag, hence the large spread shown in the relevant panels in Fig. 1 can account, to a large extent, for the spread of the observed time lags. The bottom right panel shows the spread of the 2–10 keV luminosity as well (computed as explained in Sect. 4). The X-ray luminosity also affects the magnitude of the time lag in the optical/UV bands (see e.g. Fig. 23 in K21a). Therefore, part of the scatter of the time lags plotted in Fig. 3 is also due to the spread of the L_{2 − 10keV} values that is shown in this panel. In addition, in some cases, the rest frame wavelength of the reference band, λ_ref, is not the same for a given λ_S: time lags at the same λ_S indicate delays between the observed variations in this wavelength and the variations in different λ_ref’s, depending on the redshift of each source (this is the case for sources which belong in the second and third lower boxes and in the first upper box in the top-left panel of Fig. 1).

Given all the above, modelling of the observed time lags plotted in Fig. 3 is not straightforward. In order to fit the data, we must use models which take into account the dependence of time lags on the BH mass, accretion rate, X-ray luminosity and the difference between the reference and the secondary wavelength. This is the case of the X-ray reverberation models of K21a and K21b. To demonstrate this, Fig. 4 shows the model X-ray reverberation time lags for various BH masses, accretion rates and reference wavelengths (top, middle and bottom panels, respectively). The time lags have been computed using the equations in Section 4 of K21b, for a non-rotating BH and a corona height of 20 R_g. The reference wavelengths in the top and middle panels is λ_ref = 2500 Å.

Fig. 3. The observed time lags as a function of λS. Red filled circles show the mean time lags (computed as explained in Sect. 2.3).

The top panel in Fig. 4 shows a plot of the model time lags for BH masses between ∼3.2 × 10⁷ M_⊙ and ∼8 × 10⁸ M_⊙ (this is roughly the BH mass range for most of the sources in our sample; see top right panel in Fig. 2). We have assumed an accretion rate of 0.18 in all cases (which is the mean accretion rate of the quasars in our sample). Clearly, we expect a large scatter in the observed time lags, if we consider a sample of sources with a large range of BH masses. Even if the BH mass is fixed, we still expect a (smaller) scatter in the observed time lags, if we study objects with a range of accretion rates (for reasons explained in Section 3.2 in K21a; see also their Fig. 17). This is shown by the middle panel in Fig. 4 (time lags in this figure are computed for an AGN with a BH mass of 2.6 × 10⁸ M_⊙, which is similar to the mean BH mass of the objects in the sample). The range of accretion rates we have considered is representative of the range of the λ_Edd estimates in the final sample (see the bottom right panel in Fig. 2). Clearly, a range of accretion rates will also contribute to the scatter of the observed time lags (albeit by a smaller factor, when compared to the scatter induced by the difference in BH mass).

Fig. 4. The dependence of the time lags on BH mass, accretion rate and reference wavelength (top, middle and lower panels, respectively), according to the K21a,b time-lags models (see text for details).

The bottom panel in Fig. 4 shows the expected scatter in the observed time lags due to the fact that, in some cases, the time lags at a given wavelength, λ_S, indicate delays between λ_S and various λ_ref. The model time lags in this panel are computed for M_BH = 2.6 × 10⁸ M_⊙ and λ_Edd = 0.18, while the range of reference wavelengths that we have considered is representative of the range of the reference wavelengths in the sample (see top left panel in Fig. 1).

The discussion above shows that X-ray reverberation can explain, qualitatively, the large scatter of the observed time lags in Fig. 3. However, it is necessary to also fit the model to the data to investigate whether X-ray reverberation can actually fit the observed time lags. It is not straightforward how this can be done, given the fact that the data plotted in Fig. 3 depend on various model parameters. One possibility would be to fit X-ray reverberation time-lags models to each individual source plotted in Fig. 3. However, individual time lags have large errors, which in many cases appear to be highly asymmetric (as reported in the literature).

Another possibility would be to compute the mean time lag, ${\overline{\tau }}_{\mathrm{lag},\mathrm{obs}}$, and the mean secondary wavelength, ${\overline{\lambda }}_S$, in each of the six boxes plotted in the top-left panel in Fig. 1. Since we can predict the theoretical model time lag for each source (assuming their BH mass, accretion rate, and X-ray luminosity), we could also compute the mean model time lags in the same way we compute the mean observed time lags, and then fit the mean model to the mean observed time lags. The advantage of this approach is that, by averaging many measurements, the error of the mean time lag will be significantly smaller than the error of each individual time-lag measurement and, at the same time, the mean will be approximately a Gaussian variable. As a result, we will be able to use traditional χ² statistics to fit the model to the data.

The red-filled squares in Fig. 3 show ${\overline{\tau }}_{\mathrm{lag},\mathrm{obs}}$ plotted as a function of ${\overline{\lambda }}_S$. The y-axis error indicates the standard error of the mean, while the horizontal ‘error-bar’ indicates the scatter of the secondary wavelengths around ${\overline{\lambda }}_S$. The ${\overline{\tau }}_{\mathrm{lag},\mathrm{obs}}$(${\overline{\lambda }}_S$) points in this figure define the observed, mean time-lags spectrum of the quasars in the sample, and we explain below how we fit it with the model X-ray reverberation time lags. We stress that the term ‘mean time-lags spectrum’ does not imply that the red points in Fig. 3 indicate the intrinsic, true mean time-lags spectrum of AGN as a whole. Instead, this is the ‘mean’ time-lags spectrum of this particular sample of objects (which, as we commented above, is not a complete sample of AGN), which was computed using time-lag measurements at various filters, as reported by the papers we considered.

3. Theoretical models

The X-ray reverberation time-lags models of K21a and K21b were based on the disc response functions of K21a. These functions were computed assuming an X-ray point like source illuminating a standard Novikov-Thorne (NT) accretion disc (Novikov & Thorne 1973). The X-ray source is located at a specified height, h, along the rotational axis of the BH and emits, isotropically in its rest frame, a power law spectrum with a photon index of Γ. Part of the X-ray flux that is received by the disc is reflected and re-emitted in X-rays, and the rest is absorbed. The absorbed X-rays act as an additional source of heating for the disc, leading to an increase of the local disc temperature and emission. Since the X-rays are variable, the extra disc emission should also be variable. X-ray reverberation predicts delays between the optical/UV variations, which should increase with increasing wavelength, since the X-rays initially illuminate the inner regions of the disc before propagating outwards, to the outer disc region. K21a computed model time lags, τ_mod(λ), according to the relation,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{mod}}(\lambda )=\frac{\int t\mathrm{\Psi }(t,\lambda )\mathrm{d}t}{\int \mathrm{\Psi }(t,\lambda )\mathrm{d}t},\end{array}$$(1)

where Ψ(t, λ) is the disc response function, which is defined as follows:

$$\begin{array}{c}\hfill \mathrm{\Psi }(t,\lambda )\propto \frac{F_{\mathrm{tot}}(\lambda ,t_{\mathrm{obs}})-F_{\mathit{NT}}(\lambda )}{L_X}.\end{array}$$(2)

where F_tot(λ, t_obs) in the equation above is the disc flux at wavelength λ that a distant observer will observe at time t_obs, while F_NT(λ) is the NT disc flux at the same wavelength. In effect, the disc response function is the extra, time variable disc flux due to the X-rays that are absorbed by the disc. The response functions were calculated taking into account all the relativistic effects in the propagation of light from the corona to the disc, and from the disc to the distant observer. K21a also considered the disc ionisation state in order to calculate as accurately as possible the X-ray flux that is reflected by the disc.

K21a computed model time lags using Eq. (1), and derived an analytic prescription for the time lags between optical/UV and X-rays. K21b slightly modified the K21a results, and provided an analytic model for the time lags between two wavelengths, λ_ref and λ_S, within the optical/UV bands. We use their models to fit the mean time-lags spectrum of the quasars in our sample. The main physical parameters of the model are the BH mass (in 10⁷ M_⊙), the accretion rate (in 0.05 of the Eddington limit), the height of the corona (in 10R_g), and the 2–10 keV luminosity (in units of 0.01 of the Eddington luminosity). Given these parameters, the model can predict the time lag between the disc emission at wavelengths λ_ref and λ_S within the UV/optical band, for a maximally or non spinning BH (i.e. a^* = 0³ and a^* = 1, respectively).

The K21a time lags were computed assuming a colour correction factor of f_col = 2.4. Colour correction factors are frequently used in order to take into account the effect of electron scattering of the disc photons, which may play a significant role and could lead to deviations from blackbody emission. A colour correction factor of f_col = 2.4 is appropriate in the case of AGN with a large accretion rate (Ross et al. 1992). Recently, Kammoun et al. (2023), provided an analytic correction for the K21a and K21b models, for any f_col value (see their Appendix B.1).

4. Fitting the observed time lags

We fitted the observed mean time-lags spectrum assuming six different X-ray reverberation time-lags models with a^* = 0 and f_col = 1, 1.7, 2.4 (models M01, M01.7, and M02.4, respectively), as well as a^* = 1 and f_col = 1, 1.7, 2.4 (models M11, M11.7, and M12.4). To compute the time lags when a^* = 0 we use Eqs. (2)–(5) in K21b, while for a^* = 1 we use their Eqs. (6)–(9). In order to compute the time lags for the f_col = 1 and f_col = 1.7, we added the term for the colour correction prescription, which is listed in Eqs. (B.1) and (B.2) in Appendix B.1 in Kammoun et al. (2023).

First we computed the model time lags for each of the sources in the boxes plotted in the top-left panel of Fig. 1. We used the M_BH and L_bol values from the literature (taken as described in Sect. 2.1), and we computed λ_Edd. We also computed the 2–10 keV luminosity assuming that log(L_bol/L_{(2 − 10) keV})=mlog(λ_Edd)+q, where the m = 0.75 ± 0.04 and q = 2.13 ± 0.04 (Lusso et al. 2012). Knowing M_BH, λ_Edd, a^*, and f_col, the only free parameter of the time-lags models is the X-ray corona height, h. For each source, we computed model time lags for corona heights in the range between 5 and 80 R_g, with a step of 1 R_g ⁴.

Once we computed time lags for the individual sources then, for each height, we computed the mean model time-lags spectrum in each box, ${\overline{\tau }}_{\mathrm{lag},\mathrm{mod}}({\overline{\lambda }}_S,h)$, and the χ² statistic, that is ${\chi }^2(h)={\displaystyle \sum _i}{[{\overline{\tau }}_{\mathrm{lag},\mathrm{obs}}({\overline{\lambda }}_{S,i})-{\overline{\tau }}_{\mathrm{lag},\mathrm{mod}}({\overline{\lambda }}_{S,i},h)]}^2/{\sigma }^2({\overline{\tau }}_{\mathrm{lag},\mathrm{obs},i}$), where i = 1, …, 6 (i.e. the number of boxes where we computed the observed mean time lags). The height which gives the minimum χ² (χ_min²) is the best-fit model for each of the M01, M01.7, M02.4, M11, M1.7 and M2.4 models. We note that the majority of the BH mass estimates for the sources in our sample are measured using the continuum luminosity and the width of specific emission lines in single epoch spectra. Those measurements are thought to have systematic errors of the order of ∼0.3 dex. Consequently, λ_Edd is also uncertain, even more so because L_bol is also uncertain as it is based on the use of bolometric correction factors, which have errors themselves. There is some additional uncertainty on the model time lags due to the error of the 2–10 keV luminosity as well (as parameters m and q in the Lusso et al. (2012) equations also have errors). Assuming that, on average, M_BH, λ_Edd, and 2–10 keV luminosity are uncertain by a factor of ∼2, then we estimate that the relative error of the model time lags can be up to 10–20 per cent (using the model time-lag plots in the Appendix of K21a). For that reason, we increased the error of ${\overline{t}}_{\mathrm{lag},\mathrm{obs}}$ by a factor of 1.1 when fitting the model to the observed time lags.

Fig. 5. The observed, mean time-lags spectrum (black circles, connected with a solid line) and the best-fit, M01.7 and M12.4 models (red and blue points, respectively, connected with dashed lines).

Our results are listed in Table 1. The first two columns list the various models we considered. The third column lists the best-fit corona height for each model. The respective χ_min² values are listed in the last column. Assuming that the corona height is the only free model parameter (for a given a^* and f_col), then the number of degrees of freedom (d.o.f.) is 5, and only models M01.7 and M12.4 can fit the data well (p_null = 0.04 for these models). The null hypothesis probability, p_null, is smaller than 10⁻⁵ for all other models.

Figure 5 shows the best fit results. Black filled circles connected with the solid line show ${\overline{\tau }}_{\mathrm{lag},\mathrm{obs}}$ plotted as a function of ${\overline{\lambda }}_S$ (their error is equal to the error of the mean multiplied by 1.1, as explained above). The red circles and the blue triangles, connected by the red and blue dashed lines, show the M01.7 and the M12.4 best-fit models, respectively. Both models fit well the data.

Fig. 6. Upper panel: The observed and best-fit model time lags (for model M01.7) for the individual AGN in our sample (grey and red circles, respectively). Bottom panel: The respective best-fit residuals (i.e. observed minus the model predicted time lags).

Fig. 7. The distribution of the best-fit residuals plotted in the lower panel of Fig. 6. The red-dashed line indicates a Gaussian distribution with zero mean and standard deviation equal to the mean error of the time-lag measurements we used in this work.

As a further indication that the models fit well the data, the upper panel in Fig. 6 shows the observed time lags of the individual AGN plotted as a function of λ_S (grey filled circles), and the respective model time lags for M01.7, assuming the best-fit value for the corona height of 42 R_g (red filled circles). The bottom panel in the same figure shows the best-fit residuals (i.e. observed – model time lag). The data points in both panels in Fig. 6 show that the model time lags match well, and can fully explain the range of the observed time lags.

Figure 7 shows the distribution of the best-fit residuals of the individual AGN, which are plotted in Fig. 6. The mean of the best-fit residuals is 0.36 ± 0.22 days, which shows that it is consistent with zero, as expected if the model fits the data well. Furthermore, the standard deviation of the residuals is 3.3 days. This is very similar to the mean time-lag error, which is 2.7 days (this was computed using the errors reported in the papers we used to collect data, for the ICCF time-lag measurements). We fitted the distribution of the residuals shown in Fig. 7 with a Gaussian of zero mean, and σ = 2.7 days, using the Kolmogorov-Smirnov test. The red dashed line in the same figure indicates this Gaussian distribution. The K-S statistic was 0.074, which implies p_null = 0.15, that is the distribution of the residuals is entirely consistent with a Gaussian distribution with mean zero, and standard deviation equal to the mean of the time-lag errors, as expected if the model fits the data well. This result further supports the hypothesis that the X-ray reverberation is fully consistent with the observed time lags.

Figure 6 shows that the best-fit residuals at wavelengths longer than ∼6200 Å appear to be predominately positive, indicating that the model time lags may underestimate the observations at these wavelengths. This is also shown in Fig. 5, where the observed mean time lag at the longest wavelength is ∼2σ larger than the model prediction. Although this difference is not statistically significant we note that this is not surprising. According to K21a, their model time lags are valid for time-lag measurements at wavelengths shorter than 5000 Å. For certain combinations of the model parameters, the disc may be hot enough to contribute to the disc response even at radii larger than 10⁴R_g. Since this is the largest outer disc radius that K21a considered, the model time lags could be smaller than the observed ones at long wavelengths, if the disc radius is larger than 10⁴R_g. This could explain the slight discrepancy between the observed and model time lags that we observe at λ ≥ 6200Å.

5. Summary and discussion

We have used time-lag measurements for 128 AGN, which are more luminous, and host BHs which are more massive than in nearby AGN, and we measured their rest-frame UV/optical, average time-lags spectrum. We assumed the BH mass and bolometric luminosity estimates that are listed in the papers we used to collect data, we assumed an X-ray luminosity which is typical of Type I AGN (as determined by Lusso et al. 2012), and we fitted the mean time lags with the analytical, X-ray reverberation models of K21b, taking also into account the colour correction added later on by Kammoun et al. (2023).

Our main result is that the observed time lags in quasars are fully consistent with the X-ray reverberation hypothesis. In particular, we find that the hypothesis of both spin zero and maximum spin for the BH in AGN can fit the data, as long as the colour correction factor is 1.7 and 2.4, respectively. To the best of our knowledge, this is the first time that the observed, continuum, UV/optical time lags in quasars are shown to be explained by X-ray reverberation of standard, Novikov-Thorne discs.

Our results are also consistent with the Papadakis et al. (2022) results, who studied the half-light radius versus BH mass relation in AGN, when the disc is illuminated by the X-ray corona. Contrary to past statements that the accretion disc size in quasars appears to be systematically larger than the standard thin disc model predictions by a factor of ∼2–4, Papadakis et al. (2022) showed that, when all the physical processes are taken into account, the half-light radii of X-ray illuminated discs are fully consistent with observations of gravitationally lensed quasars. To the best of our knowledge, X-ray reverberation is the only model that has been shown to fit well both the observed luminosity and accretion disc size in gravitationally lensed quasars and the optical/UV, mean time lags in quasars.

5.1. Comparison with previous work

Jiang et al. (2017), Guo et al. (2022) and Jha et al. (2022) found that the time-lag predictions provided from the X-ray illumination of the standard thin disc model were smaller by a factor of around 2 − 3 than the observed ones. Homayouni et al. (2019) also found similar results, if they would consider only the high S/N objects in their sample or only the ones with positive time lags. The results from these studies are similar to results that have been reported from the study of the continuum UV/optical time lags from the recent, multi-wavelength, intense monitoring campaigns of nearby Seyferts, when the time lags were fitted by less accurate, X-ray reverberation model time-lag relations (see e.g. Fausnaugh et al. 2016; Cackett et al. 2018; Edelson et al. 2019; Hernández Santisteban et al. 2020)

Actually, it was mainly because of the apparent disagreement between the X-ray reverberation models and the observed time lags that various other models were developed in order to explain the observations. For example, it was pointed out that reprocessed continuum emission originating in the broad-line region (BLR) gas could contribute to the wavelength-dependent continuum delays measured in AGN disc reverberation mapping experiments (e.g. Korista & Goad 2001; Lawther et al. 2018; Korista & Goad 2019; Netzer 2022). Such contribution should be expected and it is clearly evident in some (e.g. Cackett et al. 2018), but not in all sources. In our case, the mean time-lags spectrum does not show any indication of additional contribution from the BLR gas. The diffuse continuum contribution to the time-lags spectrum should be significant at wavelengths below ∼3500 Å, but we find no significant discrepancy between the model X-ray reverberation and the observed time lags at wavelengths smaller than ∼4000 Å (see Fig. 5).

Furthermore, alternative models were also developed in order to provide solutions to the (supposedly) ‘large-lags’ accretion disc problem. For example, Hall et al. (2018) studied the case when the atmospheric density in the accretion disc is sufficiently low, in which case the disc flux in the UV/optical wavelengths is lower and the discs appear to have larger size when compared to discs that emit as blackbodies (for fixed bolometric luminosity). Their model can match the AGN STORM observations of NGC 5548, but produces disc spectral energy distributions that peak at shorter wavelengths than observed in luminous AGN in general. Starkey et al. (2023) assumed a steep accretion disc height rim or rippled structures irradiated by the central X-ray source, and fitted successfully the NGC 5548 time-lags spectrum, when assuming the steep disc height rim, or multiple crests, are located at a distance ∼5 light days away from the X-ray source. We cannot comment on whether such modifications to the standard picture of the central engine in AGN (i.e. a compact X-ray source illuminating a standard, NT disc) are necessary because X-ray reverberation of canonical, NT accretion discs can fit the time lags we study. As we repeatedly shown in the past, X-ray reverberation can fit well the observed time lags in many nearby AGN as well (e.g. Kammoun et al. 2019; K21b; Kammoun et al. 2023). In fact, the X-ray reverberation model can fit the observed time lags for a large range of X-ray height and accretion rate values.

Most of the previous studies that we have referenced in the Introduction, have used analytical X-ray reverberation time-lags models which do not depend on the source height, the spin of the BH, they do not take into account relativistic effects, as well as the X-ray flux that is absorbed depending on the disc ionisation state. All these effects are taken into account by the X-ray reverberation models that we use. Perhaps the most important difference, which affects significantly the model predictions for the time lags, is that previous models assumed that the ratio of the disc heating due to X-ray absorption over the disc internal heating due to accretion is constant with radius. We suspect that the inability of the previous models to account for the observed time lags is due to the fact that this assumption is not correct. The K21b analytic models that we use in this work were developed using the K21a disc response functions, which were calculated by computing the external to internal heating at each disc radius, for any combination of the other physical parameters (i.e. BH mass, accretion rate, height, X-ray luminosity etc).

5.2. The corona height and the BH spin in quasars

Our results are based on the assumption that the X-ray source height is the same in all of AGN. We suspect this is not the case, and that there should be a distribution of heights in AGN. Our results then imply that, on average, the corona height is the same in all quasars, and it does not depend on parameters like the BH mass and/or the accretion rate. Interestingly, Papadakis et al. (2022) found that the height of the X-ray source should be larger than ∼40 R_g in order to explain the observed half light radii in micro-lensed quasars. This is exactly what we also find in our work. The corona height should be equal or larger than 40 R_g in order to explain the observed mean time lags in quasars. The fact that our model can explain these two observables, with the same model parameters, supports the hypothesis of X-ray thermal reverberation in AGN.

We note that the corona height depends on the model assumption that the disc is flat. In reality, even if H/R is roughly constant in accretion discs (where H is the disc height at radius R), we expect H to increase with increasing R. This will alter the incident angle of X-rays to the disc. It is possible then that X-ray corona with smaller heights will be able to fit the data equally well in this case. We plan to explore this possibility in the near future.

The X-ray reverberation time lags also depend on the BH spin. This was already discussed by K21a (see for example the difference in time lags for non-rotating and maximally rotating BH in their Figs. 16–23). It is for this reason that the analytical time-lags models of K21a and K21b are different for these two cases. Our results are based on the assumption that BH spin is the same in all AGN. Similarly to what we discussed above about the corona height, we expect there should be a distribution of spin parameters in AGN. Our results suggest that spin probably does not depend on other physical parameters of the system (like BH mass, or even redshift, perhaps, at least at redshifts smaller than 1). We find that both the hypothesis of zero spin and that of a maximally spinning BH can fit well the mean time lags.

We have also assumed that f_col is the same in all quasars. The best-fit model for a^* = 1 is the one where f_col = 2.4. This value is appropriate for highly accreting sources (Ross et al. 1992) and, since λ_Edd is lower than ∼0.4 for most of the sources in our sample, our results suggest that spin should not be equal to one for most quasars. On the other hand, the best-fit model for spin zero is the one with f_col = 1.7. This is very similar to f_col ∼ 1.71, which is predicted by Davis & El-Abd (2019) for the annulus with the highest effective temperature in the disc for M_BH = 2.6 × 10⁸ M_⊙ and λ_Edd = 0.18 (see their Eq. (10), for a viscosity parameter of α = 0.1). Therefore, the spin zero solution is certainly accepted, but this does not necessary imply that BHs in quasars are indeed non rotating. Our results do not exclude a non-zero spin solution, they simply appear to exclude the possibility of a maximally rotating BH in most sources in the sample.

The analytical models of K21b were developed just for two spin values, and the colour-correction factor applies to the full disc. We plan to use the codes of Kammoun et al. (2023) in a future work, in order to investigate in more detail the constrains that the existing time-lag data could impose on BH spin in quasars, as well on the colour-correction factors. It is not certain that f_col for the full disc will be well-approximated by the annulus with the highest effective temperature, and the use of Kammoun et al. (2023) could help us investigate this issue as it incorporates the colour-correction prescription of Done et al. (2012), which does not imply the same f_col for the full disc.

Acknowledgments

We thank the referee for the helpful comments. This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00. This work utilises Matplotlib (Hunter 2007), NumPy (Harris et al. 2020) and SciPy (Virtanen et al. 2020).

References

All Tables

All Figures

Fig. 1. Upper left panel: The rest frame wavelength of the reference band (in our work always the g-band), λref, versus the rest frame wavelength of the secondary light curves, λS, used to measure the time lags. Upper right panel: The MBH as a function of λS for the objects in our final sample. Bottom panels: Same as the top right panel but for λEdd and L2 − 10 keV (left and right panels, respectively). The boxes in all panels indicate the measurements we used to compute the mean time-lags spectrum.

In the text

	Fig. 2. The redshift, MBH, Lbol and λEdd distributions for the sources in the final sample.
In the text

	Fig. 3. The observed time lags as a function of λS. Red filled circles show the mean time lags (computed as explained in Sect. 2.3).
In the text

	Fig. 4. The dependence of the time lags on BH mass, accretion rate and reference wavelength (top, middle and lower panels, respectively), according to the K21a,b time-lags models (see text for details).
In the text

	Fig. 5. The observed, mean time-lags spectrum (black circles, connected with a solid line) and the best-fit, M01.7 and M12.4 models (red and blue points, respectively, connected with dashed lines).
In the text

	Fig. 6. Upper panel: The observed and best-fit model time lags (for model M01.7) for the individual AGN in our sample (grey and red circles, respectively). Bottom panel: The respective best-fit residuals (i.e. observed minus the model predicted time lags).
In the text

	Fig. 7. The distribution of the best-fit residuals plotted in the lower panel of Fig. 6. The red-dashed line indicates a Gaussian distribution with zero mean and standard deviation equal to the mean error of the time-lag measurements we used in this work.
In the text