Tracing the evolving X-ray reverberation lags within an individual AGN light curve

We present the Granger causality (GC) test for the X-ray reverberation analysis of Active Galactic Nuclei (AGN). If the light curves in the continuum-dominated band help predict (Granger cause) those dominated by reflection, the Granger lags that associate to the intrinsic reverberation lags can be inferred. We focus on six AGN observed by XMM-Newton, including the sources well-known to exhibit clear X-ray reverberation lags (IRAS 13224-3809 and 1H 0707-495) and those in which reverberation signatures are not well confirmed (MCG-6-30-15, IZW1, Mrk 704 and Mrk 1040). We employ the sliding-window algorithm and estimate the Granger (intrinsic) Fe-L lags along the light curve as the window moves through. This reveals the evolving lags towards the end of some individual observations, suggesting that the corona varies progressively. Occasionally, we observe two clearly separate lags that suggest an extended corona consisting of two zones while producing competing reverberation of two lags. While the GC test is purely hypothetical and might not explain true causality, our conclusion is that the lags are present and could be understood as reverberation lags. Assuming the lags changing solely with the corona, we find that the IRAS 13224-3809 corona varies between $\sim 10$-$25$ $r_{\rm g}$ and sometimes move to $\gtrsim 50$ $r_{\rm g}$. The corona of 1H 0707-495 and MCG-6-30-15 may be analogous to that of IRAS 13224-3809, while in IZw1, Mrk 704 and Mrk 1040 a more compact corona is expected.


INTRODUCTION
In active galactic nuclei (AGN), the central supermassive black hole accretes large amounts of matter, releasing vast amounts of energy across the electromagnetic spectrum.The X-rays are thought to primarily originate from a cloud of hot electrons close to the black hole, called the corona, through the inverse Compton scattering of low energy (optical/UV) photons from the accretion disc with the coronal electrons (Rybicki & Lightman 1986).This produces the X-ray emission which is usually called the primary X-ray continuum, whose spectrum is characterized by a cut-off power-law where its slope and cut-off energy provides information about the physical processes related to the coronal properties (Pozdnyakov et al. 1983).Some X-rays from the corona can interact with the accreting gas, resulting in the reflection spectrum where characteristic features such as photoelectric absorption and fluorescent line emission are imprinted (George ★ E-mail: pchainakun@g.sut.ac.th & Fabian 1991;Ross et al. 1999;Ross & Fabian 2005;García et al. 2014).Studying the reflection features can reveal information about the geometry and properties of the innermost region closest to the host black hole (e.g.Fabian et al. 2000;Fabian & Vaughan 2003;Reynolds & Nowak 2003).
AGN also exhibit large amplitude X-ray variability on various timescales.While the corona produces continuum X-rays, the subsequent reflection X-rays from different disc regions lead to the observed reverberation time delays due to their longer distances travelled to the observer.The delays of the reflection in response to continuum variations (i.e.reverberation lags) occur on short timescales for the inner disc reflection, and measuring these lags allows us to study the dynamics and structure of disc-corona system (Reynolds et al. 1999;Uttley et al. 2014;Cackett et al. 2021).

OBSERVATIONS AND DATA REDUCTION
The AGN sample analysed in this work was previously observed by XMM-Newton observatory and can be downloaded from XMM-Newton Science Archive (XSA). 1 Here, the same observational data of the AGN IRAS 13224-3809 and 1H 0707-495 studied in recent work (i.e., Mankatwit et al. 2023) were also alternatively analysed by a novel method proposed in this paper, along with the additional sample of the AGN MCG-6-30-15, 1ZW1, Mrk 704 and Mrk 1040; all sample data are summarised in Table 1.Indeed, the data reprocessing were basically the same with that performed in Mankatwit et al. (2023; see also Chainakun et al. 2022b).In brief, the data were reprocessed as well as removed periods which were highly affected by background flaring following the official XMM-Newton data analysis guide;2 only pn data were used to utilise their high effective area and temporal resolution.Then, the background-subtracted light curves of each observation were extracted from the circular region with encircled energy fraction of ∼ 90 per cent centred on the AGN location while the background region are the nearby, source free region; these were extracted in two different wavebands -i.e., the soft (0.3-1 keV) and hard (1.2-5 keV) energy bands -with the timing resolution of 1 second.The obtained light curves then were represented as the observational data for analysing by method purposed in this work.Note that for the observational ID 0743050301 and 0743050801 of the AGN IZW1, each observation was not operated continuously, but was split into five segments of observation; this is because the original purpose of these observations were optimized for high-resolution spectroscopic study of outflows from the AGN (PI E. Costantini; see also Wilkins et al. 2017).This results in the observational data of ∼ 20 ks for each segment, and, in this study, we analyzed each individual segmental light curve separately -i.e., treating each segment as individual data similar to that was obtained from different observations -to avoid the data gaps between the segments.
In addition, to allow us to compare the quality of real data with that of the simulated data, we also determined the quality of each dataset -i.e.how much the source signal is affected by the actual background level -by calculating the signal to noise ratio (S/N) defined as: where  is the count rate, being equivalent to count () divided by the exposure time ().Here, the S/N in the hard band was calculated, and used to represent the quality of each dataset, as shown in column 6 of Table 1.

GC test
The association between the time series can be studied using the the GC, by determining whether there is a relationship in which one time series precedes (or Granger-causes) the other (Granger 2004).We define the soft (reflection dominated, 0.3-1 keV) and hard (continuum dominated, 1.2-5 keV) band light curves to be   and ℎ  , respectively.Let us first assume that they are simply generated from the autoregressive process given by:  (3) The summation term describes how the current observation depends on the  preceding value, while   and   are the corresponding coefficients for the soft and hard bands, respectively.  and  ℎ are parameters that provide new information into both energy bands.The prediction errors (uncertainties or residuals) which occurred in the soft and hard bands are denoted by   , and   ,ℎ , respectively.
To test for the Granger causality under the X-ray reverberation framework, the   is modified so that it also contains the past values of ℎ  .This can be given by a bivariate regressive process: where  is a constant and   , is the error from this estimating model.Note that the parameter  is the upper limit of the lag value we want to test and some of the coefficients   and   can be zero.While the primary power-law component from the coronal emission is the main constituent in any of the X-ray bands,   is usually more dominated by the reverberation flux than ℎ  (e.g. the fraction of reflection found in   is relatively large compared to what found in ℎ  ).We then expect   to lag behind ℎ  .In other words, ℎ  should Granger-cause   (ℎ → ) since the observed   depends on past values of ℎ  .
Here, we perform ℎ →  analysis by testing if any lag values of ℎ  are statistically significant in the bivariate regressive model.All the lag values of 1, 2, 3, . . .,  are used to estimate the two regression equations which are eqs. 2 and 4. The corresponding errors obtained in both equations for each lag value are evaluated using the GC indicator to see if the error   , from eq. 4 is significantly smaller than the error   , from eq. 2 (i.e. if the past values of ℎ  help predict   ).
In other words, the GC is assessed by comparing the variances of the errors from both bivariate and univariate models.When the past value of ℎ  does not predict   ,   (  , ) ≈   (  , ) so  ℎ→ = 0, meaning that ℎ  does not Granger-cause   .Contrarily, if ℎ  Granger-causes   , so ℎ  helps improve the prediction for   , then   (  , ) <   (  , ) and  ℎ→ increases.
Note that the GC test must be performed only when the light curves are stationary (i.e.having a constant mean, variance, and coherence, and no seasonal component).This is because the GC test infers the lags based on building a forecasting model (i.e.statistical testing if there are the causal linkages where the past values of ℎ() help predict () with any particular lag values).The non-stationary time series variables should not be included in a regression model, otherwise a problem known as spurious regression may arise.Therefore, most time series forecasting models assume stationary in the first place.
Here, the X-ray light curves,   and ℎ  , are simulated from the red noise PSD (via stingray.simulator(Huppenkothen et al. 2019)) convolving with the disc response functions, which should link to the X-ray variability mechanism in AGN.The autoregressive model is used later as a forecasting model only to predict   and ℎ  after being transformed to stationary ones (normally, stationary series are easier to analyze by simpler forecasting models).In practice, the stationary light curve means that we minimize its components dependent on time, so that the forecasting time series model can be made.We follow the method outlined in Chainakun et al. ( 2023) where differencing together with the unit root test via Augmented Dickey-Fuller (ADF) method are used to ensure the stationarity.The GC tests are carried out using the grangercausalitytests module in statsmodels.tsa.stattools(Seabold & Perktold 2010).We employ the -value to infer the significance of the Granger lags.Statistical hypothesis tests are arranged such that  < 0.01 means that the obtained lags are significant.

Sliding window
A sliding window is used to perform the GC test on different segments of the light curve.The window moves through the light curve, from left to right, allowing the Granger lags to be estimated for each segment.The window size determines the length of the segment considered at each specific time, whereas the moving step-size determines the amount of shift applied to the sliding window as it progresses each step.The choice of step size (e.g.overlapping or non-overlapping) depends on several factors and may affect the computational resolution and accuracy of the specific problem.
To avoid redundant computations, we define the window size to be a fixed value equal to the half length of the light curve, as illustrated in Fig. 1.If the light curve contains  time bins, each of which has the size of Δ, the window size is then Δ/2.When the Granger lags for the current window have been analyzed, the window is moved by one position to the right, with the moving step equal to the bin size Δ.We do not consider the stationarity for an entire light curve or for all segments simultaneously, but instead look at the GC lag for each individual segment where the sliding-window moves through one by one segment (i.e. the mean, variance and autocovariance should be constant within each segment covered by the window, but different segments can have different statistical properties).Therefore, the estimated lags between two light curves can be varied across all segments, and stacking them together can reveal how time lags evolve continuously along the light curve within a single XMM-Newton orbit.

Error estimates
For the GC test in the previous section, we obtain the Granger lags of each observation with the condition of the hypothesis test (e.g.-values).Different significant lag-values can be found based on different choices of binning.The error of the lag may be estimated by, for example, dividing the light curve for each window further into several sub-segments.Then, the lag is given by the average of these sub-segments, and the standard error of the mean lag is quoted as an error.This method is similar to the one used for the crosscorrelation function which does not always work especially when the light-curve segments are short.Therefore, it might not be obvious how to estimate the error.For consistency, we use the time-binning size as the approximate error for the Granger lags.For example, if the time bin is 100 s (i.e. the Granger tests are performed with the time step of 100 s with  = 1, 2, 3, . . .), the obtained lags can possibly be 100 s, 200 s, 300 s, . . .with the estimated error of +/− 100 s.

Identifying the lags in the light curve
To illustrate the Granger-lag profile used in our analysis, we construct 50 ks light curves with a binning of 1 s from the primary variation   generated by the red noise processes in stingray.simulator(Huppenkothen et al. 2019).The reverberation response is modelled, for simplicity, using the top-hat response function   rev whose profile is controlled by the centroid and the width which are set to be  = 600 s and   = 50 s, respectively.On the other hand, we include the response due to the disc propagating fluctuations that operate on relatively long timescales in the hard band using the top-hat function   prop , with  = 1000 s and   = 200 s.The area under these responses are normalized to 1.The soft reflection-dominated and hard continuum-dominated band light curves are computed by where   and  ℎ represent the fraction of reverberation and discpropagating fluctuations contributed in the soft and hard bands.We introduce the effects of random noise to the light curves by adding uncorrelated variability with random mean and the standard deviation using numpy.random.normal.The signal-to-noise ratio (S/N) is estimated from the mean count rate over the average error in the way that enables us to compare it with the S/N of the real data (eq.1).The light curves are binned with the bin size of Δ = 50 s.Note that the bin size will affect the significance of the lags.If the bin size is chosen to be comparable to the width of the response, the characteristic response profile will become close to the delta-function and the Granger lags can be interpreted as the average reverberation time determined by the centroid of the top hat response.In reality, we do not know in advance the width of the reverberation response and binning the data too large may remove information about the width of the response.In our current approach, the variable bin-size is applied to investigate all possible significant lags.We perform the ℎ →  test and analyze the lags following the method outlined in Chainakun et al. (2023).Briefly speaking, we investigate all lag values  and the corresponding errors obtained from the two regression equations (eqs. 2 and 4).The -value is used to justify if the error from eq. 4 (where   also contains the past values of ℎ  ) is significantly smaller than that from eq. 2 (where   is explained by its own past values).We plot each lag value  in terms of the corresponding -values in order to infer their significance.
Fig. 2 shows examples of the Granger-lag profiles obtained from simulated 10 pairs of the light curves.It can be seen that the reverberation-lag amplitude ( = 600 s) could be estimated from the minimum Granger-lag at  ≤ 0.01, especially when the S/N ≳ 100.The inferred reverberation lag is also an intrinsic one since the dilution can only make the lags become less significant but does not alter the value of the lags.Note that the -value can stay smaller than 0.01 even at longer values of the targeted lag (i.e. at > 600 s here).This is because the Granger causality is a statistical test to determine whether "at least" one of the lags of ℎ  up to that particular value is significant to explain   .In other words, when at least one of the lags of < 600 s of ℎ  can help predict   , it is also true that at least one of the lags of < 800 s can help explain the data.This is why the -values for long Granger lags can be < 0.01, and why we should use the minimum Granger lag that is first significant as a proxy of the reverberation lags.However, the two lags can still be probed using variable bin sizes, such as in Fig. 6 of Chainakun et al. (2023).The random noise and effects of the competing process such as the disc propagating fluctuations that operate on much longer timescales affects more on the significance of the lag than on manifesting the lag values (see also Appendix A and B in Chainakun et al. 2023).

Identifying the lags using the sliding window
Rather than considering an entire light curve to have a single lag value, the sliding window is applied to estimate the lags along the light curve.The Granger-lags along the light curves are probed using several bin sizes, so that the results can be compared and discussed.The bin sizes used in this illustration start from the smallest size of Δ  = 100 s and increase with the step size given by  step = 10 s until it reaches the maximum size of Δ  = 500 s.All lags found when the window moves through the sequence are recorded and plotted against the the mid-time of the sliding window.We simulate the light curve as before and assume the   rev of  = 500 s in the soft band.The result is shown in Fig. 3.It can be seen that the almost constant Granger-lags of ∼ 500 s can be detected along the segments where the window moves through.We employ the Kernel Density Estimation (KDE) function available in seaborn (Waskom 2021) to estimate the density of the population of the obtained lags with a Gaussian kernel3 .The middle and bottom panels of Fig. 3 represents the contour lines depicting the density levels of the lag distribution and the corresponding smoothed version of the density, respectively.It is normalized such that all the probabilities for all the lag values found in this map add up to 1, meaning 100 per cent of the lag falls under the distribution plot.The units of the density are arbitrary where lighter colors show higher density values for the obtained lag distribution.The almost constant lags that cluster around ∼ 500 s is analogous to, e.g., the reverberation lags from the lamp-post corona staying at a fixed height on the symmetry axis of the black hole.Note that we also investigate the effects of the chosen bin size in Appendix A).Now, we investigate a dual lamp-post case intended to explain the vertical extent of the corona by adding a second top-hat function (  rev,2 ) mimicking the reverberation response due to a second (upper) lamp-post source (Chainakun & Young 2017;Lucchini et al. 2023).The soft reflection-dominated band light curve is then given by where  ,1 and  ,2 regulate the importance of the reflected flux between the lower and the upper source, respectively.For the lower source, there is a larger amount of photons incident the disc so  ,1 >  ,2 .We simulate the light curve using  1 = 600 s and  ,1 = 80 s for   rev,1 , and  2 = 1200 s and  ,2 = 50 s for   rev,2 .The obtained Granger-lag profiles are presented in Fig. 4. The results demonstrate that a dual lamp-post corona can produce the competing lags in two separate zones specifically to the lags due to each X-ray source.Note that the width of the response as well as the factor  ,1 and  ,2 affect the significance of the lags, but the lags can and when the top-hat response representing reverberation and propagation lags is included in the soft and hard bands, respectively.We also include uncorrelated variability to show how the estimated lags depend on the S/N of the light curves.For the reverberation response, we fix the centroid of the response time to be 600 s.The reflection fraction is fixed at  = 1.The significance level of  ≤ 0.01 is shown in the horizontal dashed line.The intrinsic lags of 600 s can be estimated using the minimum lag-value at the point where the curve crosses the significance level line.It can be seen that the error in estimating the lags increases with decreasing the S/N ratio, and the preferred S/N in this method is approximately more than 100.
be approximated as the intrinsic reverberation lags with amplitudes given by  1 and  2 .Smaller  ,2 does not decrease the lag amplitude but it can lead to the less significance of the detected time lags as in Fig. 4. the sliding window is used to record the lag value seen in that segment.When the mid-time of the sliding window changes from 49 ks to 70 ks, the lags are found to decrease from ∼ 1570 s to ∼ 700 s.By applying the sliding window and binning the data using various sizes of the time bins, we can identify the lags continuously along the light curve as well as the distribution of the significant lags dependent on the choice of the time bins.once the sliding window is applied.We use the KDE to highlight how the distribution of the lags is found and how the lags change across the light curve.Here, the trend of decreasing lags towards the end of the light curve is evidenced.Next section, we explore the general trends of these lags for all observations of our targeted AGN.

Global features of variable lags
General trends of the Granger-lags seen in the light curves can be summarized as an example in Fig. 7.The lags in some observations of these AGN remain almost constant as time evolves.These include, e.g., IRAS 13224-3809 rev. nos. 3043 and 3044, and 1H 0707-495 rev. no. 1972.Their corona then should be located at an almost constant height above the black hole during the observation.Meanwhile, some observations show increasing and decreasing trends of the lags towards the end of the observation, analogous to the lags from the dynamic corona changing its height.
The irregular trends of the variable lags are also seen in a number of particular observations.For example, in IRAS 13224-3809 rev.no.2127, we see two separate lags at ∼ 500 s and ∼ 1200 s.Perhaps this suggests the corona is substantially extended, clustering in two zones unfolded before us as the competing reverberation of two lags.Interestingly, in the case of Mrk 704 rev.no.1630, the increasing 3053 where the window of fixed size moves through (top panel).In this example, the lags (  < 0.01) from two specific sequences covering by red and blue windows are found to be ∼ 1570 s to ∼ 700 s, respectively (middle panels).The obtained lags are plotted against the light-curve time given by the mid-time of the corresponding sliding window (bottom panel).While the time-bin size is used as the approximate error for the lags, we note that the corresponding KDE plot (as shown in Fig. 6) should also provide the idea of the distribution of all possible lags estimated with this method.
lags are observed but we see the hint of separate, competing lags that may suggest the extended corona.

Evolving lags in 6 AGN
From the visual representations of the data, we attempt to determine the trend of the evolving lags by considering the Spearman's rank correlation coefficient (  ).We analyze whether the Granger-lag values generally increase (or decrease) as time progresses for each individual light curve.Fig. 8 shows the obtained monotonic correlations between the lags and the light curve time for all AGN investigated here.Only the results with  < 0.01 are presented.The correlations are significant in 1H 0707-495, IRAS 13224-3809, 1Zw1, MCG-6-30-15, and Mrk 704.When   > 0 (or   < 0), the lags significantly increase (or decrease) towards the end of the observation.We find that the lags significantly increase (decrease) with time in 8 (or 11) observations of these AGN.

DISCUSSION AND CONCLUSION
A traditional way to study the X-ray reverberation delays is to analyze them in the Fourier-frequency domain by calculating, e.g., the lagfrequency spectrum (Fabian et al. 2009;Zoghbi et al. 2010;Wilkins & Fabian 2013;Cackett et al. 2014;Chainakun & Young 2015;Caballero-García et al. 2018) where the use of the full light curve is usually necessary.In fact, a dynamic of the X-ray corona was robustly constrained in, e.g., IRAS 13224-3809 by considering the lag observed in each of the full light curves from 16 XMM-Newton observations (Alston et al. 2020;Caballero-García et al. 2020;Chainakun et al. 2022b).However, analysing the lags from the full light curve can relegate them to the average value, hence preventing us from studying how the lags evolve throughout the individual observation.This work builds on the GC analysis on the intrinsic X-ray reverberation reported by Chainakun et al. (2023) where the hint of the coronal evolution towards the end of some individual observations of IRAS 13224-3809 was suggested.
Here, the study is expanded to six AGN and the sliding window is applied to probe the variable lags along the light curve.We attempt to determine and visualize the trends of the evolving lags by considering the distribution of the obtained lags via the KDE (Waskom 2021).The lags seen along the individual light curve where the sliding window moves through are likely variable (Fig. 7).This suggests a presence of the dynamical corona in the individual observations of these AGN.The simplest way to induce such intrinsic-lag variation is that the corona changes its height along the rotational axis of the black hole, as explained by the light bending model (Miniutti & Fabian 2004).
From 14 XMM-Newton observations of IRAS 13224-3809, Chainakun et al. (2023) found the Granger-lags are significant in 12 observations while the evolving lags are seen only in 4 observations.By adopting variable bin sizes and applying the slidingwindow technique, we find the variable and significant lags in all 14 IRAS 13224-3809 observations.Furthermore, we find the tendency of the lags to significantly increase (or decrease) towards the end of the observation in 8 (or 11) light curves from six AGN (Fig. 8).When this trend is not observed, the lags tend to concentrate around a specific value suggesting a corona which is less variable in size and location but, perhaps, still moving up or down along the rotational axis.Nevertheless, for the majority of the cases, it is very likely that the scatter of the obtained lags are intrinsic to the extended corona.This is because the distributions of the lags revealed to us look very different from one observation to another and may cluster in more than one zone.This can be inferred as having in some sense a competing reverberation process due to different coronal parts.
While the exact shape of the corona is under debate, it may vary among different AGN or different observations.According to Fig. 7 (bottom left panel), we see in IRAS 13224-3809 rev.no.2127 two almost-persistent but separate lags at ∼ 500 s and ∼ 1200 s along the light curve.These features may be captured by a vertically extended corona simplified using two lamp-post sources as demonstrated in Fig. 4. Therefore, this may simply be explained by dual lamp-post model (Chainakun & Young 2017;Hancock et al. 2023;Lucchini et al. 2023) where there are two X-ray blobs inducing two separate reverberation lags.In terms of the GC test, the two lags from dual lamp-post source are not average as one, but instead show that both lag values can provide useful information of reverberation that statistically help explain the light curve.A more complex and extended corona (e.g.Wilkins et al. 2016;Chainakun et al. 2019) may be required to explain a more complex trend of the evolving lags.
The IRAS 13224-3809 corona constrained using different methods, such as the lag-frequency spectra (Alston et al. 2020), the combined spectral-timing model Caballero- García et al. (2020), the power spectral density analysis (Chainakun et al. 2022b) and the machine learning approach via the random forest regressor (Mankatwit et al. 2023), was quite consistent in that it was found to vary between ∼ 5 − 20  g .A study on IRAS 13224-3809 using a dual lamp-post model also suggested the corona could extend up to ∼ 20  g (Hancock et al. 2023).Here, the intrinsic lags of ∼ 200-500 s are found in all IRAS 13224-3809 observations, and sometimes the lags can be as large as > 1000 s in particular segments of each individual light curve.In some observations the lags may stay almost constant at < 500 s towards the end of the observation (e.g.rev.nos.3043 and 3044 in top panels of Fig. 7).Assuming a face-on disc and the central mass of 2 × 10 6  ⊙ (Alston et al. 2020), we can convert the intrinsic lags directly to the true light-travel distance.The inferred distance suggests that the corona spends most of its time at ∼ 10-25  g and may occasionally move up to ≳ 50  g , consistent with Chainakun et al. (2023).The disc density as well as the X-ray luminosity that determines the ionization of the disc may affect the dilution of the lags, but they should not significantly affect the intrinsic lags that are probed by the Granger test as long as the geometry is the same.However, the Granger lags interpreted here in the context of the height-changing corona is a rough approximation since the disc can have a geometric thickness that varies among different AGN (Taylor & Reynolds 2018).It is still possible that the changes in Granger lags observed here may be the results of the combined effects due to the changes in both corona and disc geometry.
For 1H 0707-495, the central mass can possibly be ∼ 2.3 − 10 × 10 6  ⊙ (e.g.Zhou & Wang 2005;Zoghbi et al. 2011;Pan et al. 2016) while the observed Granger-lags are in the same range of those seen in IRAS 13224-3809.With the larger mass and the comparable intrinsic lags, the coronal height in 1H 0707-495 should be smaller than that in IRAS 13224-3809.Wilkins & Fabian (2012) investigated the emissivity profiles using various coronal geometries and IRAS 13224-3809 (rev. no. 3043) IRAS 13224-3809 (rev. no. 3044) 1H 0707-495 (rev. no. 1972) 1H   suggested a presence of the 1H 0707-495 corona at ∼ 2  g extending radially outwards to ∼ 30  g .Recent studies on the relativistic reflection spectra of 1H 0707-495 suggested that its corona should be very compact and located at ℎ ∼ 3  g with an extended size of ≲ 1  g (Szanecki et al. 2020).A proposed corona at ∼ 3-4  g was also consistent with Dauser et al. (2012) andCaballero-García et al. (2018).However, the coronal extent must be large enough to receive the seed photons from the accretion disc in order to produce the continuum X-rays.Regarding this, Dovčiak & Done (2016) proposed that the 1H 0707-495 corona may locate at ∼ 2 − 5  g and its size may extend outwards to ∼ 20  g during its maximum flux.Recently, Hancock et al. (2023) showed that the frequency-dependent time lags of 1H 0707-495 could be explained using a dual lamp-post model.They fixed the lower source height at 2  g and found the upper coro-nal height varying at ∼ 3-20  g .Our finding also supports that the 1H 0707-495 corona is dynamic and extended as in IRAS 13224-3809, but with smaller physical distance due to their mass difference.
Previous studies based on the lag-frequency spectra suggested a compact corona for MCG-6-30-15 that is located at ∼ 3  g (Emmanoulopoulos et al. 2014;Epitropakis et al. 2016).The extended corona model has never been adopted to formally fit their lagfrequency spectra before.However, our results conclude that the corona in MCG-6-30-15 should also be varied since we observe the variable Granger-lags as in IRAS 13224-3809 and 1H 0707-495.Interestingly, the averaged lags of MCG-6-30-15 are mostly within ∼ 500-1500 s, similar to IRAS 13224-3809.The MCG-6-30-15 mass is also ∼ 2 × 10 6  ⊙ (Ponti et al. 2012) which is comparable to the IRAS 13224-3809 mass.This further suggests that its size and geometry may be analogous to that of IRAS 13224-3809.
In the case of IZw1, the Fe-K reverberation lags were probed using the lag-energy spectra, but no Fe-K reverberation feature was detected (Kara et al. 2016).Nevertheless, Wilkins et al. (2021) applied a different method to the standard Fourier technique and detected the Fe-K lag amplitude of ∼ 746 ± 157 s due to the X-ray flares occurring near the black hole.Here, the evidence of the X-ray reverberation lags in IZw1 is also revealed but in the Fe-L band, with the lags varying mostly within ≲ 1000 s.The IZw1 mass is approximately It is important to note that the GC method is a purely hypothetical test that is physical-model independent.The GC test by itself does not provide an explanation of true causality, so it does not provide information of the origin of these lags.Here, we perform theoretical simulations to show that these lags can possibly be produced by reverberation.The results reveal that the lags are variable in a complex way.The main parameter driving the variable lags might be intrinsic to the corona itself.The size and geometry of the corona may vary progressively even within one individual observation.All of these results support a dynamic corona that, sometimes, can clus-ter in more than one zone, resulting in the separate competing lags.Using more than two light curves would help reveal more additional hidden factors influencing the target light curve.
The GC technique is a statistical test so justifying the errors of the lags is not straightforward.In fact, the AGN light curves points can be correlated to past values, so the time-domain statistics themselves can be correlated.Therefore, one cannot easily distinguish if it is the light curves themselves being variable or it is the variations of the intrinsic lags that mainly produce the variation of the Granger lags (i.e. the lags from reverberation are indistinguishable from one that happens by chance).It could be useful in the future to take into account the GC test results in conjunction with some of the statistical requirements in order to develop a robust physical model behind the test.Here, the meaning and significance of the lags rely on a given interpretation specific to the reverberation framework.Therefore, the observed variable-lags obtained by the GC test should be thought in the way that they all have a possibility to be produced by reverberation.The results from this technique still needed to be considered together with prior knowledge and known information from Fourier-based time-lags to draw a complete conclusion.Also, the GC technique provides a single number of the lag, hence the information about the structure of the cross-correlation between two light curves is limited.Improving the GC technique itself might help evade this issue but this is beyond the scope of this work.Last but not least, we note that the specific implementation details of the GC test as well as the sliding window algorithm can be further advanced by, e.g., applying it to directly fit the data or using a variable windowsize specifically for any targeted segment of the light curve.New and improved statistical analysis may provide us a better way to robustly detect the significance of the lag.

Figure 1 .
Figure 1.Sketch of moving window with a step size equal to the bin size of the light curve, Δ.The size of the window is fixed at  Δ/2.

Figure 2 .
Figure 2. Granger-lag profiles (ℎ → ) from 10 simulated pairs of light curves plotted in terms of the -value when the light curves are identical (top-left panel)and when the top-hat response representing reverberation and propagation lags is included in the soft and hard bands, respectively.We also include uncorrelated variability to show how the estimated lags depend on the S/N of the light curves.For the reverberation response, we fix the centroid of the response time to be 600 s.The reflection fraction is fixed at  = 1.The significance level of  ≤ 0.01 is shown in the horizontal dashed line.The intrinsic lags of 600 s can be estimated using the minimum lag-value at the point where the curve crosses the significance level line.It can be seen that the error in estimating the lags increases with decreasing the S/N ratio, and the preferred S/N in this method is approximately more than 100.

Figure 4 .
Figure 4. Same as in the bottom panel of Fig. 3, but when the soft band contains two top-hat response functions of  = 600 s and 1200 s representing a dual lamp-post case.While the fraction of the reflected flux of the lower source is fixed at  ,1 = 0.6, that of the upper source is varied to be  ,2 = 0.3 (top panel) and 0.4 (bottom panel).A dual lamp-post X-ray corona can produce two separate lags appearing in the Granger-lag profiles.

Figure 5 .
Figure 5. Illustration of the sliding window used to process and analyze the Granger-lags from each sequence of the light curve of IRAS 13224-3809 rev.no.3053where the window of fixed size moves through (top panel).In this example, the lags (  < 0.01) from two specific sequences covering by red and blue windows are found to be ∼ 1570 s to ∼ 700 s, respectively (middle panels).The obtained lags are plotted against the light-curve time given by the mid-time of the corresponding sliding window (bottom panel).While the time-bin size is used as the approximate error for the lags, we note that the corresponding KDE plot (as shown in Fig.6) should also provide the idea of the distribution of all possible lags estimated with this method.

Figure 6 .
Figure 6.Top panel: Granger-lag values (data points) detected along the light curve of IRAS 13224-3809 rev.no.3053 and the contour lines showing the density distribution of the data.Bottom panel: Corresponding smoothed density distribution of the lags where lighter colors show higher density values.We can see the tendency of the lags that decrease towards the end of the observation.

Figure 7 .
Figure 7. General trends of the distribution of the Granger-lags found as seen in our light curves.The panels in the first, second, third and fourth rolls represent the (almost) constant, increasing, decreasing and irregular trends, respectively.Note that they are selected and catalogued by eyes.All sources presented here have the S/N ≳ 100 except IRAS 13224-3809 rev.no.3043 (S/N = 94), and 1H 0707-495 rev.no.159 (S/N = 52) and 1387 (S/N = 43).Higher density values for the lag distribution are plotted in lighter colors.See text for more details.

Figure 8 .
Figure 8. Distribution of the obtained Spearman's rank correlation coefficient (  ) for the correlation between the Granger-lags and evolving time.Only the samples that have significant correlation (  < 0.01) are included.The significant GC lags correlated with the light curve time can be found in all AGN except Mrk 1040.  > 0 means the lags increase towards the end of the observation, and vice versa.

Δ𝑡Figure A1 .
Figure A1.Density maps showing the distribution of the Granger-lags detected at a significance level of  < 0.01 along the light curve of IRAS 13224-3809 rev.no.2131.Note that the units on the density are arbitrary and higher probability-density values are associated with lighter colors.Top panels: Bin size of the light curve is varied between Δ  = 50 s and Δ  = 250 s, with the step size of  step = 10 s (left panel) and 20 s (right panel).Bottom panels: Bin size is varied between Δ  = 100 s and Δ  = 300 s, with the step size of  step = 10 s (left panel) and 20 s (right panel).Adjusting the bin sizes used in the investigation may affect the Granger-lag profiles, but the trend of how the evolving lags appear as well as the cluster of the lag confidence do not significantly change.

Table 1 .
XMM-Newton observations of the AGN sample.
The data S/N in the hard band defined by eq. 1. * The observation was split into five individual segments which each one has the GTI of ∼20 ks (see text for further detail).