© ESO, 2016

1. Introduction

The study of time-lags as a function of temporal frequency between X-ray light curves in different energy bands has frequently been used in the last two decades to probe the emission mechanism and geometry of the innermost region in active galactic nuclei (AGN; e.g. Papadakis et al. 2001; McHardy et al. 2004; Arévalo et al. 2006, 2008; Sriram et al. 2009) and X-ray binaries (XRBs; e.g. Miyamoto & Kitamoto 1989; Nowak & Vaughan 1996; Nowak et al. 1999). In the last few years, time-lag studies have revealed that “soft” X-ray variations in AGN are delayed with respect to “hard” X-ray variations at frequencies higher than ≈ 10^-4 Hz (e.g. Fabian et al. 2009; Zoghbi et al. 2010, 2011; Emmanoulopoulos et al. 2011; De Marco et al. 2013). These time-lags are commonly referred to as soft lags, and have been observed in ≈ 20 sources so far. They have attracted considerable attention, since they are thought to be direct evidence of X-ray reverberation. Further credence to such a scenario came with recent discoveries of time-lags between the Fe Kα emission line (≈ 5−7 keV) and the X-ray continuum (e.g. Zoghbi et al. 2012; Kara et al. 2013b,a,c; Zoghbi et al. 2013b; Marinucci et al. 2014), as well as between the so-called Compton hump (≈ 10−30 keV) and the X-ray continuum (e.g. Zoghbi et al. 2014; Kara et al. 2015) in a few AGN.

The study of X-ray reverberation time-lags can provide valuable geometrical and physical information on the X-ray source and reflector, since they should depend on, for example, their typical size, relative distance, proximity to the central black hole (BH), the mass and spin of the BH, and the inclination of the system. To obtain this information, the statistical properties of time-lags estimated from observed light curves must be known as accurately as possible. For example, one must know their bias (i.e. how accurately their mean approximates the intrinsic time-lag spectrum), error, and probability distribution. The later is necessary if one wishes to fit the observed time-lag spectra with theoretical models. To the best of our knowledge, such a detailed investigation has not yet been performed.

Results from preliminary studies along these lines have been presented by Alston et al. (2013) and Uttley et al. (2014). In this paper we report the results of a more detailed study regarding the statistical properties of Fourier-based time-lag estimates, based on both analytical work and extensive numerical simulations. Our primary goals are: a) to investigate whether the frequently used time-lag estimates are indeed reliable estimates of the intrinsic time-lag spectrum; b) to study the effects of light curve sampling patterns and duration, as well as measurement errors, on the statistical properties of these estimates; and c) to provide practical guidelines which ensure estimates that are minimally biased, have know errors, and follow a Gaussian distribution. The latter property would be desirable, for example, in the case of model fitting using traditional χ² minimisation techniques. Our work should have direct impact on current time-lag studies in the area of AGN X-ray timing analyses. We believe that our results can apply equally well to all scientific fields where similar techniques are employed to search for delays between two observed time-varying signals.

2. Definitions

Consider a continuous bivariate random time series { X(t),Y(t) }. We assume that the mean values (μ_X and μ_Y), as well as the auto-covariance functions (R_X(τ) and R_Y(τ); τ is the so-called lag) of the individual time series are finite and time-independent (i.e. they are stationary random processes). A random function that is frequently used to quantify the correlation between two time series in the time domain is the cross-covariance function (CCF), (1)where E is the expectation operator. We assume that the CCF depends on τ only, i.e. the CCF does not vary with time. Defined as above, R_XY(τ) is the CCF with Y(t) leading X(t). The Fourier transform of the CCF, (2)defines the cross-spectrum (CS) of the bivariate processes. The CCF is not necessarily symmetric about τ = 0, and hence the CS is generally a complex number that can be written as (3)The real functions c_XY(ν) and {−q_XY(ν) } represent the real, ℜ [ h_XY(ν) ], and imaginary, ℑ [ h_XY(ν) ], parts of the CS, respectively. The function c_XY(ν) is an even function of ν, while q_XY(ν) is an odd function of ν (Priestley 1981, P81 hereafter). The quantity | h_XY(ν) | is the CS amplitude, and φ_XY(ν) is the phase-lag spectrum, which is defined as (4)The phase-lag spectrum represents the average phase shift between sinusoidal components of the two time series with frequency ν. Since it is defined modulo 2π, it is customary to define its principal value in the interval (−π,π ], which is the convention adopted in the present work as well. The time-lag spectrum is defined as (5)and represents the average temporal delay between sinusoidal components of the two time series with frequency ν. Given the definition of the CCF by Eq. (1), a positive time-lag value at a frequency ν indicates that, on average, the sinusoidal component of X(t) lags behind the respective component of Y(t).

Another statistic that is often used to study correlations between two random processes in Fourier space is the so-called coherence function, (6)where h_X(ν) and h_Y(ν) denote the power-spectral density functions (PSDs) of the time series X(t) and Y(t), respectively (the PSD is defined as the Fourier transform of the auto-covariance function of a random process). The coherence function is interpreted as the correlation coefficient between two processes at frequency ν, as it measures the degree of linear correlation between them at each frequency. As with ordinary correlation coefficients, , where indicates a perfect (linear) correlation, and implies the absence of any correlation at frequency ν.

2.1. The effects of binning and discrete sampling

In practice, the data correspond to values of a single realisation of a random process that are recorded over a finite time interval, T (the duration). The recording is typically performed at regular time intervals, Δt_sam (the sampling period).

Consider a pair of observed time series { x(t_r),y(t_r) }, where t_r = rΔt_sam, r = 1,2,...,N, and N is the total number of points. Let us denote by X(t) and Y(t) the intrinsic, continuous time series, which we assume are stationary random processes (henceforth, the term intrinsic will always refer to { X(t),Y(t) }). The observed time series { x(t_r),y(t_r) } correspond to a particular, finite realisation of a discrete version of the intrinsic process, which we denote by { X(t_s),Y(t_s) }, where t_s = sΔt_sam and s = 0, ± 1, ± 2,.... In addition, light curves¹ in astronomy are the result of binning the intrinsic signal over time bins of size Δt_bin². In this case, the relation between X(t_s) and X(t) is given by (7)with an identical relation holding between Y(t_s) and Y(t).

It is rather straight-forward to show that { X(t_s),Y(t_s) } and { X(t),Y(t) } have the same mean values, but different CCFs and CS. As we show in Appendix A, the CS of the discrete process, h_XY(ν), which is defined only at frequencies | ν | ≤ 1/2Δt_sam, is given by (8)In other words, at each frequency ν within the aforementioned interval, h_XY(ν) is the superposition of the intrinsic CS values, h_XY(ν), at frequencies ν,ν ± (1/Δt_sam),ν ± (2/Δt_sam),.... This is entirely analogous to the so-called aliasing effect in the case of the PSD. However, although PSDs are always positive, hence aliasing in this case always implies transfer of “power” from higher to the sampled frequencies, the effects of aliasing on h_XY(ν) cannot be predicted without a priori knowledge of h_XY(ν). Aliasing affects both c_XY(ν) and q_XY(ν), which are not necessarily positive at all frequencies. The situation is further complicated by the fact that while c_XY(ν) is an even function of ν, q_XY(ν) is an odd function of ν. Therefore, aliasing should affect these two functions in different ways.

Equation (8) shows that aliasing is reduced if the light curves are binned, since in this case the values of h_XY(ν) that are aliased in the sampled frequency range are suppressed by the sinc² function, whose argument depends on the bin size, Δt_bin. Since the sinc² function asymptotically approaches unity in the limit whereby its argument goes to zero, if the observed time series are not binned (i.e. when Δt_bin → 0) the sinc² term in the right-hand side of Eq. (8) vanishes, and hence the effects of aliasing on h_XY(ν) are maximised.

3. Estimation of cross-spectra; the cross-periodogram

The discrete Fourier transform (DFT) of a time series x(t_r) is defined as (9)where is the sample mean of x(t_r). It is customary to estimate DFTs at the following set of frequencies: ν_p = p/NΔt_sam, where p = 1,2,...,N/ 2, and ν_Nyq = 1/2Δt_sam is the so-called Nyquist frequency, which corresponds to the highest frequency that can be probed for a given sampling period, Δt_sam (the symbol ν_p will henceforth always stand for this particular set of frequencies). The cross-periodogram of two time series is defined as (10)where the asterisk ^∗ denotes complex conjugation. The cross-periodogram is used in practice as an estimator of the CS. Given Eqs. (4) and (5), it seems reasonable to use the real and imaginary parts of I_xy(ν_p) and accept as estimators of the phase- and time-lag spectrum, respectively.

3.1. The bias of cross-spectral estimators

As we show in Appendix B, (13)where h_XY(ν) is the CS of { X(t_s),Y(t_s) } (as defined by Eq. (8)). The function F_N(ν′ − ν_p), which is called the Fejér kernel, has a large peak at ν′ = ν_p of magnitude NΔt_sam, and decays to zero as | ν′ | → ∞. Subsequent peaks appear at ν′ ≈ ν_p ± 3/2Δt_sam,ν_p ± 5/2Δt_sam,..., while the zeros of the function occur at ν′ = ν_p ± 1 /NΔt_sam,ν_p ± 2 /NΔt_sam,.... As N → ∞, F_N(ν′ − ν_p) → δ(ν′ − ν_p) (the Dirac δ-function). Therefore, in the limit N → ∞ the right-hand side of Eq. (13) converges to h_XY(ν_p). The cross-periodogram is therefore an asymptotically (when N → ∞) unbiased estimate of the modified (owing to the effects of binning and discrete sampling) intrinsic CS.

Since h_XY(ν) is not equal to h_XY(ν), it follows that φ_XY(ν) (τ_XY(ν)) will generally not be equal to φ_XY(ν) (τ_XY(ν)) either. Even if aliasing and binning effects are minimal (so that we can assume that φ_XY(ν) ≈ φ_XY(ν)), in practice it is difficult to predict whether the duration T = NΔt_sam of a given pair of observed time series will be sufficiently long for the effects of the convolution of h_XY(ν) with the Fejér kernel (Eq. (13)) to be significant of not, unless h_XY(ν) is known a priori. We define the “bias” of the cross-periodogram as b_I(ν_p) ≡ E [ I_xy(ν_p) ] − h_XY(ν_p) (henceforth, the term bias for a statistical estimator will always refer to the difference between its mean and intrinsic value). Even if b_I(ν_p) was known, and hence could be used to “correct” I_xy(ν_p), it would still not be straightforward to determine the bias of or , since E { arg [ I_xy(ν_p) ] } is not necessarily equal to arg { E [ I_xy(ν_p) ] }.

To quantify the bias of time-lag estimates based on the cross-periodogram, we performed an extensive number of simulations, which we describe below. The characteristics of the simulated time series (i.e. time bin size, sampling rate, duration, intrinsic CS, and PSDs) are representative of those observed in AGN X-ray light curves. However, most of our results should apply to any time series observed in practice (see discussion in Sect. 10).

4. Simulating correlated random processes

The simulations we performed are based on the procedure outlined by Timmer & Koenig (1995)³ to generate artificial realisations of a discrete stationary random process with a specified model PSD, h_X(ν), number of points, N, sampling period, Δt_sam, and mean count-rate, μ_X. We assumed a model PSD of the form (14)This function describes a power law with a low-frequency slope of −1 which smoothly “bends” to a slope of −2 at frequencies above the bend-frequency, ν_b, as in the case of inferred X-ray PSDs of most AGN. We chose an amplitude value of A = 0.01, and assumed that ν_b = 2 × 10^-4 Hz. Furthermore, we set N = 10.24 × 10⁶ and Δt_sam = 1 s. By construction, the intrinsic PSD of the generated light curves is discrete, and has non-zero values only at frequencies ν_j = j/NΔt_sam, where j = ± 1, ± 2,..., ± N/ 2.

For each simulated light curve we also created a corresponding “partner” with a mean count-rate μ_Y, assuming a specified model phase-lag spectrum, φ_XY(ν). This was achieved by multiplying the DFT of the first light curve by the factor (μ_Y/μ_X)e^{− iφ_XY(ν)}. The light curves thus have the same PSD shape, and can be considered as realisations of a discrete process whose intrinsic CS is equal to h_XY(ν) = (μ_Y/μ_X)h_X(ν)e^iφ_XY(ν). Since the CS, in effect, represents the average product of the Fourier component amplitudes of each light curve, and this amplitude is proportional to the PSD, the CS is, by construction, non-zero only at the same frequencies where the PSD is non-zero as well. The light curve pairs generated following this procedure will henceforth be referred to as original realisations.

4.1. The model time-lag spectra

We considered three different model time-lag spectra:

• 1.

  Constant delays: in this case we assumed that φ_XY(ν) = 2πνd, where d = 10 s, 150 s and 550 s is aconstant delay between the light curve pairs (henceforth,experiments CD1, CD2, and CD3, respectively).

• 2.

  Power law delays: in this case we assumed that φ_XY(ν) = (2πν)Bν^{− β}, where B is the normalisation and β the power-law index of the corresponding model time-lag spectrum. We considered the values { B,β } = { 0.001,1 }, { 0.01,1 }, { 0.1,1 }, { 0.01,0.5 }, { 0.01,1.5 } (henceforth, experiments PLD1, PLD2, PLD3, PLD4, and PLD5, respectively).

• 3.

  Top-hat response functions: in this case we assumed that φ_XY(ν) = arg [ 1 + fe^i2πνt₀sinc(πνΔ) ]. This phase-lag spectrum is expected if the intrinsic time series are related by the following equation: , where Ψ(t) (the so-called response function) is a simple so-called top-hat function, i.e. it has a constant value of f/ Δ in the interval | t − t₀ | ≤ Δ/2 and zero otherwise (Δ is the width of the top-hat and t₀ its centroid). We considered the model parameter values { f,t₀,Δ } = { 0.2,200 s,200 s } and { 0.2,2000 s,2000 s } (henceforth, experiments THRF1 and THRF2, respectively).

4.2. The sampling patterns of the simulated light curves

For each model time-lag spectrum listed above, we created 30 light curve pairs as per the specifications given in Sect. 4. To simulate light curves encountered in practice, the original realisations need to be properly “chopped”, binned, and sampled.

Most of the data in X-ray astronomy are provided by satellites in low-Earth orbit with a typical orbital period of ≈ 96 min and bin size of 16 s, such as ASCA, RXTE, Suzaku and NuSTAR. Light curves obtained from observations with these satellites are “affected” by periodic Earth occultations of a target during every orbit. As a result, they contain gaps that are typically ≈ 1−3 ks long. They are hence not appropriate for Fourier analysis using the “standard” techniques we considered in this work. One can bin the data at one orbital period to acquire evenly sampled light curves with no missing points. Such light curves can be used to probe variability on long time-scales. On the other hand, the data usually contain an appreciable number of continuously sampled segments, with a duration typically ≲3 ks. These segments can, in principle, be used to probe variability on short timescales.

Contrary to low-Earth orbit satellites, XMM-Newton observations result in continuously sampled light curves up to ≈ 120 ks long owing to its highly elliptical orbit. We can use such a pair of light curves to compute the cross-periodogram and time-lag estimates directly, or “chop” them into segments of shorter duration, calculate the cross-periodogram for each segment and then bin the resulting estimates at certain frequencies to estimate the time-lags (see Sect. 6 for a more detailed discussion on this issue).

The purpose of the simulations we performed is to study the sampling properties of time-lags when estimated using light curves with durations and sampling patterns similar to those described above. To this end, we adopted the following strategy:

• 1.

  Chop each of the original realisations into100 parts and bin them at 100 s. This process generates3000 light curve pairs with T = 102.4 ks and Δt_sam = Δt_bin = 100 s (LS102.4 lightcurves, hereafter).

• 2.

  Chop each of the LS102.4 light curves into 2/5/10/ 20 segments. This process generates 6000/15 000/30 000/ 60 000 light curve pairs with T = 40.8/20.4/10.2/5.1 ks and Δt_sam = Δt_bin = 100 s (LS40.8/20.4/10.2/5.1 light curves, hereafter).

• 3.

  Chop each of the original realisations into 33 × 96 = 3168 parts and bin them at 16 s. This process generates 95 040 light curves with T = 3.2 ks and Δt_sam = Δt_bin = 16 s (LS3.2 light curves, hereafter).

• 4.

  Chop each of the original realisations into 33 parts and bin them at 5760 s (≈ 96 min). This process generates 990 light curves with T = 305.3 ks and Δt_sam = Δt_bin = 5760 s (OB light curves, hereafter).

The reason for the original realisations being longer and with a finer sampling rate than the “final” light curves is to simulate the effects of binning and finite light curve duration on the bias of the time-lag estimates, (henceforth, the time-lag bias).

5. The bias of the time-lag estimates in practice

To quantify , we calculated using Eq. (12) for the 10 numerical experiments and each light curve type described in Sects. 4.1 and 4.2. We then computed the sample mean at each frequency, (quantities in angle brackets will hereafter denote their sample mean). We chose to study the properties of mainly, since the majority of works concerned with AGN X-ray timing studies use this estimator. In addition, we computed the quantity to quantify the time-lag bias in terms of its intrinsic value (henceforth, the relative bias). This allows us to directly compare the time-lag bias obtained from different light curve types in the various numerical experiments we considered.

Figures D.1–D.4 show for the LS102.4/LS3.2/OB (top rows; continuous black, red, and brown curves, respectively) and LS40.8/20.4/10.2/5.1 (bottom rows; continuous black, red, brown, and green curves, respectively) light curves. Each column in these figures corresponds to a different numerical experiment. The main results of the simulations we performed are summarised below.

5.1. Light curve binning and aliasing

A common feature in all mean sample time-lag spectra is that they decrease to zero at high frequencies, irrespective of the light curve sampling pattern, duration, and model time-lag spectrum. This decrease is due to the effects of aliasing and light-curve binning.

To demonstrate this issue, we considered experiment PLD2 (the results are similar for the other numerical experiments as well). We generated two additional ensembles of light curve pairs that have the same length as the LS102.4 light curves. The first ensemble was constructed using the original realisations of experiment PLD2, which were chopped into 3000 parts of length 102.4 ks each, and sampled (not binned) every 100 s (henceforth, the LS102.4-2 light curves). These light curves are still affected by aliasing, since the original realisations have a finer sampling rate of 1 s. For the second ensemble, we constructed 30 new original realisations with the same length as the rest of the numerical experiments (i.e. 10.24 Ms) and a sampling period of 100 s. We subsequently chopped them into 3000 parts of length 102.4 ks each (henceforth, the LS102.4-3 light curves). Since the original realisations in this case have the same sampling rate as the LS102.4-3 light curves, they should not be affected by either binning or aliasing.

We then calculated and for these two additional light curve types. Figure 1 shows ⟨ ℜ [ I_xy(ν_p) ] ⟩, ⟨ ℑ [ I_xy(ν_p) ] ⟩, , and (top left and right, bottom left and right panels, respectively) for the LS102.4, LS102.4-2, and LS102.4-3 light curves (black, red, and brown curves, respectively)⁴. The sample mean time-lag spectra for the LS102.4-2 and LS102.4 light curves decrease to zero at frequencies higher than ≈ 10^-3 Hz (= ν_Nyq/ 5) and ≈ 2.5 × 10^-3 Hz (= ν_Nyq/ 2), respectively. These frequencies are indicated by the vertical dashed lines in the same figure. This is not the case for the LS102.4-3 light curves, as shown by the brown curves in the bottom panels of Fig. 1. The time-lag bias of the sampled and binned light curves can be understood by the plots in the top panels of Fig. 1. On average, the gradients of the imaginary and real parts of the cross-periodogram enter steeper rates of decrease and increase, respectively, compared to their intrinsic values, around the frequencies indicated by the vertical dashed lines in Fig. 1. As a result, their ratio (which determines the phase-lag estimate) is decreased on average, hence the time-lag bias increases. This increase is most severe at high frequencies, where the effects of aliasing and light curve binning on the cross-periodogram are maximised. The time-lag bias owing to these effects is more pronounced for the LS102.4-2 light curves, since light curve binning suppresses the aliasing effect (see Sect. 3).

Figures D.1–D.4 show that the onset of the high-frequency increase in the relative time-lag bias occurs at ≈ ν_Nyq/ 2 for all light curve types, model time-lag spectra, and light curve durations. To further illustrate this effect, we define the function as the mean relative time-lag bias at all frequencies below ν_max. Figure 2 shows , evaluated at ν_max = ν_Nyq/ 3, ν_Nyq/ 2.5, ν_Nyq/ 2, and ν_Nyq/ 1.5, for the OB (diamonds), LS40.8 (squares), and LS102.4 (circles) light curves in experiments CD1 (black), PLD2 (red), and THRF1 (blue). The vertical red and black dashed lines indicate ν_Nyq/ 2 for the OB and LS102.4/40.8 light curves, respectively. At frequencies below ν_Nyq/ 2, the relative time-lag bias remains approximately constant, but increases markedly at higher frequencies (the same trend is observed in all numerical experiments we considered).

We conclude that time-lags should be estimated at frequencies ≲ν_Nyq/ 2 to minimise the effects of light-curve binning and aliasing on the time-lag bias.

5.2. Finite light curve duration

The time-lag estimates are biased even at frequencies ≲ν_Nyq/ 2. Figures D.1–D.4 show that the relative time-lag bias is generally ≲15% at frequencies ≲ν_Nyq/ 2 for the LS102.4/40.8/20.4 and OB light curves, and larger for the rest. It is usually positive, in the sense that the estimates are, on average, smaller in absolute magnitude than their corresponding intrinsic values. The bias at intermediate/low frequencies is not due to light curve binning and/or aliasing (we note, for example, that at frequencies ≲ν_Nyq/ 5 in Fig. 1 is identical for all three light curve types, which have the same length but different time bin sizes). This bias is due to the finite duration, T, of the light curves, or, technically speaking, due to the convolution of h_XY(ν_p) with the Fejér kernel (see Eq. (13)).

To quantify the dependence of this bias on T, Fig. 3 shows as a function of T for the LS5.1/10.2/20.4/40.8/102.4 and OB light curves (filled and open points, respectively) in all experiments that do not exhibit so-called phase-flipping (see Sect. 5.3 for details). The dashed red line indicates the best-fit relation to the LS-type data: . This relation fits the LS-type data quite well. We therefore conclude that the mean relative time-lag bias decreases with increasing light curve duration as . However, this relation is not consistent with the points corresponding to the OB light curves. This is unexpected at first, since their duration (≈ 300 ks) is larger than the duration of all LS-type light curves. This discrepancy arises because the frequency range probed by the OB light curves is significantly lower than by the LS-type light curves (owing to their longer duration and larger time bin size), as well as the fact that the relative time-lag bias increases at lower frequencies.

To illustrate this effect, we generated additional ensembles of light curve pairs that have the same time bin size (100 s) and a duration two, three, four, and five times longer than LS102.4 light curves in the case of experiment CD1 (the results are identical for the other numerical experiments as well). The solid lines in the top panel of Fig. 4 show the relative time-lag bias at all frequencies below ν_Nyq/ 2 for all these light curves (their duration is listed next to each curve). As expected, the relative bias at each frequency decreases with increasing T. For a fixed T, the relative bias increases with decreasing frequency. The filled boxes in the same plot indicate the mean relative bias over the whole frequency range, plotted at the mean logarithmic frequency. The mean relative bias, over the full frequency range, decreases with increasing T. In the bottom panel of the same figure we plot the same results in the case of OB-type light curves with a duration of ≈ 50−500 ks. The relative bias at each frequency is almost identical to the previous case, despite the larger time bin size of the OB-type light curves. Although the mean relative bias decreases with increasing T as before, owing to the larger time bin size in this case, we sample a lower frequency range, and hence the mean relative bias is larger than before (for the same T). We therefore conclude that the relative time-lag bias depends on , although the normalisation of this relation depends on the frequency range on which the time-lags are estimated; the lower the frequencies, the more biased the time-lag estimates will be. For the model CS we considered, our results suggest that, to obtain time-lag estimates with a relative bias ≲15%, the light curves must have T ≳ 20 ks.

5.3. Phase-flipping

Since the phase-lag spectrum is defined on the interval (−π,π ], there can be frequencies whereby this function will exceed the value of π and “flip back” to the value of −π (and vice-versa). This effect is known as phase-flipping (or phase-wrapping).

These kinds of events take place in experiments CD2, CD3, and PLD5 (see Figs. D.1 and D.3). The top left panel of Fig. 5 shows for the LS102.4 light curves in experiment CD3 (continuous black line), along with the corresponding model time-lag spectrum (black dashed line). The simulated light curves in this experiment are separated by a constant delay. As a result, the model phase-lag spectrum is a linearly increasing function of frequency, i.e. φ_XY(ν) = 2πνd. At ν = 1/2d ≈ 9.1 × 10^-4 Hz we get φ_XY(ν) = π and then, by definition, φ_XY(ν) “jumps” to the value −π. Since τ_XY(ν) = φ_XY(ν)/2πν, the model time-lag spectrum will likewise jump from 550 s to −550 s. Subsequent phase-flips occur at frequencies ν = j/ 2d (j = 2,3,...), where the model time-lag spectrum undergoes discontinuous jumps from positive to negative values with decreasing amplitude.

The mean sample time-lag spectrum also fluctuates from positive to negative values at these frequencies, although the transition is much smoother. To illustrate the reason for this behaviour, in the top right and bottom panels of Fig. 5, we plot the probability distribution of the phase-lag estimate at the following frequencies: 8.1 × 10^-4 Hz (top right panel), 9.1 × 10^-4 Hz and 1.0 × 10^-3 Hz (bottom left and right panel, respectively). These frequencies are indicated by the vertical dashed lines in the top left panel of the same figure. The middle frequency is very close to the frequency where the first phase-flip occurs, while the other two are lower and higher than that frequency. The solid and dashed vertical lines in the top-right and bottom panels of Fig. 5 show the model and sample mean phase-lag value at these frequencies, respectively. Since the phase-lag estimate is defined on the interval (−π,π ], the so-called wings of its distribution that exceed these boundaries (indicated by the “empty” histograms in the plots) are “folded back” into the allowed range. This causes the mean of the distribution to shift towards a value lower than that of the model phase-lag spectrum. This effect is most severe in the vicinity of frequencies where phase-flipping occurs, since, in this case, the mean of the phase-lag estimate is close to zero.

Phase-flipping does not necessarily occur at increasingly higher frequencies. If, as in the case of experiments PLD4 and PLD5 for example, the intrinsic phase-lag spectrum increases (in absolute magnitude) with decreasing frequency, phase-flipping takes place at increasingly lower frequencies. This is seen in the bottom left panel of Fig. 6, where we plot for the LS102.4 light curves in experiment PLD5 (open black circles). The model time-lag spectrum for this experiment is τ_XY(ν) = 0.01ν^-1.5 (indicated by the continuous brown line). The corresponding phase-lag spectrum is φ_XY(ν) = (2πν)0.01ν^-1.5, and phase-flipping occurs at (4 × 10^-4/j²) Hz (j = 1,2,...). At frequencies close to 4 × 10^-4 Hz, where the first phase-flip occurs, the mean sample time-lag spectrum exhibits a jump which is smoother than the abrupt jump of the model time-lag spectrum, for exactly the same reasons we discussed above. This is not very clear in the plot shown in the bottom left panel of Fig. 6, but is evident in the plot of (bottom-right panel in the same figure).

The dashed blue line in the top panels of Fig. 6 indicate the model real and imaginary parts of the CS. Interestingly, the sample mean of the real and imaginary parts of the cross-periodogram (indicated by the open circles in the top panels) are not biased at frequencies around ≈ 4 × 10^-4 Hz. In fact, the argument of their ratio divided by the angular frequency, i.e. arg { E [ ĥ_XY(ν_p) ] } /2πν_p, is very similar to the model time-lag spectrum, and exhibits a sharp jump at this frequency. This is a case where, for the reasons explained above, arg { E [ ĥ_XY(ν) ] } is not equal to E { arg [ ĥ_XY(ν) ] }.

The time-lag estimates at the two lowest frequencies are heavily biased owing to the bias of the real and imaginary parts of the cross-periodogram. This is reflected in the mean sample time-lag spectrum, which has a large relative bias at those frequencies (bottom right panel in Fig. 6). This bias originates from the convolution of the intrinsic CS with the Fejér kernel, which causes the mean real and imaginary parts of the cross-periodogram to diverge from the model CS because of its rapid oscillatory behaviour at low frequencies.

6. Smoothed/averaged cross-spectral estimators

The variance of the real and imaginary parts of the cross-periodogram is and , respectively (P81). They are unknown, and do not depend on the number of points in the observed time series (i.e. they will not decrease if we increase their duration). In practice, we usually average a certain number, say m, of consecutive cross-periodogram estimates, i.e. (15)and accept ĥ_xy(ν_k) as the CS estimator at the frequencies ν_k = (1 /m) ∑ _pν_p (k = 1,2,...,N/m). The process of binning over consecutive frequencies is called smoothing⁵. The real and imaginary parts of the smoothed CS estimator are correspondingly given by An alternative procedure is to partition the available time series into m shorter segments of duration T/m, compute the cross-periodogram for each segment, (l = 1,2,...,m, and ν_k = k/ (T/m), where k = 1,2,...,N/m), and then average the different cross-periodogram values at each ν_k: (18)with the real and imaginary parts estimated as in the case of the smoothed estimates (to differentiate between estimators that are smoothed or averaged over individual segments, we henceforth refer to the former as smoothed and the latter as averaged estimates).

One can show that the variance of the real and imaginary parts of the smoothed/averaged cross-periodogram (as well as their covariance) is inversely proportional to m. As the duration of the observed light curves increases (i.e. N increases for a fixed Δt_sam), m can be increased proportionally without degrading the frequency resolution. In this sense, the variance of the smoothed/averaged estimates decreases with increasing N (in the limit N → ∞, it approaches zero).

We can use the real and imaginary parts of the smoothed/averaged CS estimates to construct the following estimators of the phase- and time-lag spectrum: Their variance is given by the following asymptotic formulae (P81; Nowak et al. 1999; Bendat & Piersol 2011) Equations (21) and (22) highlight the importance of in constructing a reliable time-lag estimator. The lower the coherence, the higher the variance of the phase- and time-lag estimates will be (this is expected, as it must be more “difficult” to detect delays when the two light curves are highly incoherent). The “natural” choice for the coherence estimator, as suggested by Eq. (6), is the following: (23)where ĥ_x(ν_k) and ĥ_y(ν_k) are the smoothed/averaged PSD estimators of the two observed time series. The variance of the coherence estimator is given by the asymptotic formula (P81; Vaughan & Nowak 1997; Nowak et al. 1999; Bendat & Piersol 2011) (24)Since the intrinsic coherence, , is unknown, it is customary to replace it in Eqs. (21) and (24) by its estimate, , to obtain a numerical value.

6.1. Bias due to smoothing

To investigate whether smoothing or averaging increases the bias of the cross-spectral estimates, we smoothed the real and imaginary parts of the cross-periodograms estimated from the LS102.4 and OB light curves (for all numerical experiments) using Eqs. (16) and (17) for m = 20. We also used the cross-periodograms of the LS20.4/40.8 light curves, and averaged m = 20 of them to produce the corresponding averaged estimates. We then computed, in both cases and for all experiments, and using Eqs. (20) and (23), as well as their corresponding sample mean values, and . Figures D.5−D.8 show for the smoothed (top panels; open white circles and filled brown squares for the LS102.4/OB light curves) and averaged estimates (bottom panels; continuous black and red curves for the LS20.4/40.8 light curves).

In general, smoothing increases the time-lag bias at low frequencies. The top panel in Fig. 7 shows vs. for the smoothed estimates (in the case of the numerical experiments that do not exhibit phase-flipping in the sampled frequency range). The bottom panel shows the same plot for the averaged estimates. The blue dashed line shows the one-to-one relation in each case. Although the bias of the averaged estimates falls on the expected one-to-one relation line, this is not always the case for the smoothed estimates. Significant, additional bias appears in the case of experiments CD1, PLD4, and THRF1. The intrinsic CS in these cases exhibits a prominent non-linear variation with frequency. As a result, linear smoothing introduces a bias to the real and imaginary parts of the smoothed cross-periodogram. This is demonstrated in Fig. 8, which shows ⟨ ℜ [ ĥ_xy(ν_k) ] ⟩ (top left panel), ⟨ ℑ [ ĥ_xy(ν_k) ] ⟩ (top right panel), (bottom left panel), and (bottom right panel) for the m = 20 smoothed estimates of the LS102.4 light curves in experiment CD1. The intrinsic CS and time-lag spectrum are shown as blue dashed lines. At low frequencies (≲2 × 10^-4 Hz) where the intrinsic CS has the highest curvature (i.e. where its second derivative has the largest value), the bias of the smoothed cross-periodogram and time-lag estimate is maximised.

The mean sample coherence at all frequencies is ≈ 1 in most, but not all, cases. Figure D.9 shows the mean sample coherence for experiments CD2, CD3, and PLD5 (filled circles and solid lines indicate the m = 20 smoothed and averaged estimates, respectively). Despite the fact that, by construction, at all frequencies, the mean sample coherence is significantly smaller than unity in these cases. This is difficult to explain and predict. The coherence bias is a complicated function of the and bias, as well as the bias of the smoothed/averaged PSD estimates. For example, although the mean real part of the smoothed cross-periodogram is biased in experiment CD1 when m = 20 (see the top left panel in Fig. 8), the mean sample coherence is not, presumably because it is counterbalanced by the bias of the smoothed PSD estimates.

7. The effect of measurement errors

Every measured signal will inevitably be affected by experimental “noise”. To study the effect of measurement errors on the time-lag and coherence estimates, we considered the LS20.4 and LS40.8 light curve pairs in experiments CD1, PLD2, and THRF1. For these experiments, we increased the number of original realisations from 30 to 100. This increase in number of simulated light curves for each numerical experiment results in a significant increase in the computing time. For that reason, in this and the following section we consider only experiments CD1, PLD2, and THRF1. They do not exhibit phase-flipping effects, and are representative of the three categories of model phase-lag spectra we consider. In total, we thus ended up with 5 × 10⁴ LS20.4 and 2 × 10⁴ LS40.8 light curve pairs for each experiment.

Furthermore, we created five copies of each of these pairs to simulate the effects of different signal-to-noise ratio (S/N) combinations, { (S/N)_x,(S/N)_y }, for each experiment. The effects of measurement errors were simulated by adding a Gaussian random number with zero mean to each point of every light curve, which were uncorrelated both with each other, as well as with the intrinsic light curve values. The variance of the random numbers was chosen such that there were five different S/N combinations for each light curve pair in the aforementioned experiments: { (S/N)_x,(S/N)_y } = { 3,3 }, { 9,3 }, { 18,3 }, { 9,9 }, and { 18,9 }.

We then calculated the m = 10, 20, 30, and 40 averaged cross-periodogram to estimate the time-lags and coherence according to Eqs. (20) and (23), along with their error estimates according to Eqs. (22) and (24). The number of averaged time-lag and coherence estimates were thus 5000 (2000), 2500 (1000), 1666 (666), and 1250 (500) for the LS20.4 (LS40.8) light curves in each experiment and every S/N combination. Figures D.10−D.12 show (top row) and (bottom row) for the LS20.4/40.8 light curves (solid and dashed lines, respectively) for all S/N combinations, in each of the three experiments and for m = 20.

7.1. The effect of measurement errors on the time-lag bias

A comparison between, for example, the top panels in Fig. D.10 and the bottom left panel in Fig. D.5 reveals that the time-lag estimates become increasingly biased at high frequencies (in the sense that they converge to zero) with decreasing S/N. As the S/N of the light curves increases, becomes consistent with its corresponding value when there is no noise present. To illustrate the effect of noise on the time-lag bias, we computed as a function of for experiment THRF1 (this quantity is a measure of the combined S/N of both light curves). The results are shown in Fig. 9 for the LS20.4 and LS40.8 light curves. The horizontal dashed black and dotted-dashed red line correspond to the value when there is no noise (i.e. when (S/N)_xy → ∞) in the LS20.4 (filled black circles) and LS40.8 (open red squares) light curves, respectively. Clearly, the time-lag estimates are more biased when (S/N)_xy ≲ 3. Identical trends are observed in the case of experiments CD1 and PLD2 as well.

This is a rather unexpected result, since the noise introduced to one light curve is independent of the noise introduced to the other. We would therefore not expect the time-lag bias to be affected by noise (although we would expect the resulting estimates to have a larger scatter around the mean, which is indeed the case; compare, for example, the scatter of the time-lags in the top panels in Fig. D.10 and in the bottom left panel of Fig. D.5). As in the case of phase-flipping (see Sect. 5.3), this bias arises because the phase-lag estimates are defined on the interval (−π,π ]. As an example, in Fig. 10 we show the probability distribution of phase-lag estimates for the LS40.8 light curves with (S/N)_x = (S/N)_y = 9 in experiment THRF1 at three different frequencies: 6.1 × 10^-4 Hz, 1.2 × 10^-3 Hz, and 2.5 × 10^-3 Hz (top, middle, and bottom panel, respectively). Filled black bars and open red bars indicate the distribution of the m = 10 and m = 40 averaged phase-lag estimates, respectively. As the frequency increases, the mean sample coherence decreases (see the next section), and hence the scatter of the phase-lag estimate increases according to Eq. (21). When the magnitude of this scatter becomes sufficiently large, the “wings” of the distribution exceed the range of allowed values for the phase-lag estimate, and are thus folded back into the allowed range. Consequently, the distribution becomes increasingly uniform over the interval (−π,π ], and hence its mean converges to zero. This bias reduces as m increases, since, in this case, the scatter of the phase-lag estimates themselves, and hence their bias, decreases (see Eq. (21)). This is also evident from Fig. 10, as the “width” of the distributions is smaller in the case of larger m.

The vertical dashed lines in all panels of Figs. D.10−D.12 indicate the frequency at which the sample coherence becomes equal to 1.2/(1 + 0.2m) (this value is indicated by the horizontal dotted-dashed lines in the lower panels of the same figures). We refer to this frequency as the critical frequency, ν_crit, and in Sects. 8 and 9 we discuss its importance in detail. Interestingly, we found that, on average, the time-lag bias is similar to its value in the absence of measurement errors at frequencies below ν_crit. The same result holds for the time-lag estimates in the case of m = 10, 30, and 40 as well.

7.2. Effect of measurement errors on the sample coherence

Measurement errors have the intuitive effect of decreasing the intrinsic coherence between two time series. By construction, the intrinsic coherence of all light curves is equal to one at every frequency. As we show in Appendix C, noise decreases the intrinsic coherence at all frequencies. This decrease is minimal at frequencies where the amplitude of the measured signal’s intrinsic variations dominates over the amplitude of the noise variations. In the opposite case, the mean sample coherence will be significantly less than unity.

This is exactly what we see in the bottom panels of Figs. D.13−D.17; the sample coherence is always less than unity at high frequencies, even in the highest S/N cases. The decrease depends on the S/N; at a given frequency, the mean sample coherence decreases with decreasing S/N. By fitting various functions to the sample mean coherence we found that, for each numerical experiment and both light curve types, the following simple exponential function describes them well (see Appendix C for a theoretical justification of this fact): (25)This is demonstrated in Fig. 11, which shows the mean sample coherence for the LS40.8 light curves (m = 20) in experiment THRF1 for various S/N combinations (we get the same results in the case of experiments CD1 and PLD2 as well). The red dashed lines in this figure indicate the best-fit of the function defined above.

The blue horizontal dashed lines in the same figure indicate the value 1 /m. As we demonstrate in Appendix C, at frequencies where experimental noise dominates the intrinsic variations, the mean sample coherence is expected to be equal to that value. This is indeed what we observe in our simulations, as in the case of the lower S/N light curves the sample coherence indeed converges to the value of 1 /m (see the top panel in Fig. 11).

8. The error of the time-lag estimates

The top panels in Figs. D.13–D.17 show the mean sample coherence for the LS40.8 and LS20.4 light curves and various S/N combinations. Columns 1–4 show the results for m = 10, 20, 30, and 40, respectively. The plots show the mean sample coherence for all experiments where we considered the effects of noise (CD1, PLD2, and THRF1). They are difficult to identify, because their values are practically identical for all S/N combinations. In the middle panels, we plot the so-called error ratio, which we define as , where is the mean analytical error estimate (computed using Eq. (22) by substituting the value of the intrinsic coherence by its estimate), and is the observed 1σ scatter of the averaged time-lag estimates around their mean. The bottom panels in the same figures show the probability , i.e. the probability that the value of the time-lag estimate falls within the range (i.e. within 1σ of the sample mean).

The black, red, and green lines in the same figures indicate the results obtained from experiments CD1, PLD2, and THRF1, respectively. Continuous and dashed lines correspond to the LS40.8 and LS20.4 light curves. For a given m and S/N combination, the results are almost identical for the different experiments and light curve types. This indicates that the error of the time-lag estimates depends mainly on m and the S/N combination, irrespective of the intrinsic CS and light curve duration.

The ratio of the analytical error over the observed time-lag standard deviation is larger than ≈ 0.9 and approaches unity as m increases, at frequencies which are lower than a certain “critical” frequency, say ν_crit. This frequency is indicated by the vertical dashed lines in all panels of Figs. D.13−D.17. In the same frequency range, the probability is close to 0.68 (this value is indicated by the horizontal dotted lines in the bottom panels of the same figures). This is what we would expect if the distribution of the time-lag estimates is approximately Gaussian.

The critical frequency is independent of the intrinsic CS and on the light curve duration, and depends mainly on the light curve S/N and on m; ν_crit increases with increasing m and S/N. We found that ν_crit can be estimated by equating the mean sample coherence to the value 1.2/(1 + 0.2m), i.e. (26)This coherence value is indicated by the blue dotted-dashed horizontal line in the top panels of Figs. D.13−D.17.

At higher frequencies, the analytical error underestimates the true scatter around the mean of the time-lag estimates, and so the error ratio decreases. As we show in Appendix C, the mean sample coherence is always larger than the coherence of the “noisy” processes at all frequencies. Since the analytical time-lag error estimate depends on the sample coherence (see Eqs. (21) and (22) and the comment at the end of Sect. 6), it is not surprising that it underestimates the true scatter around the mean, and that the error ratio is smaller than unity. The mean sample coherence becomes significantly larger than the coherence of the “noisy” processes at frequencies ≳ ν_crit. Consequently, the analytical error estimate begins to increasingly underestimate the true scatter around the mean and hence the error ratio decreases at the same frequencies.

Interestingly, at frequencies ≲ν_crit the probability also remains constant and close to 0.68 (i.e. the value expected for a Gaussian distribution). At higher frequencies it begins to decrease, indicating that the distribution of the time-lag estimates develops longer “tails” than would be expected from a Gaussian distribution. We investigate the properties of the probability distribution of the time-lag estimates in more detail below.

9. The probability distribution of the time-lag estimates

The results presented above regarding the error of the time-lag estimates imply that, at frequencies below ν_crit, their distribution may be, approximately, Gaussian. To further investigate this issue, we estimated the excess kurtosis and skewness of the time-lag estimates distribution, for the three experiments and two light curve types. Figures D.18–D.22 show our results (upper and middle panels for the excess kurtosis and skewness, respectively). Black, red, and green lines indicate the results obtained from experiments CD1, PLD2, and THRF1, respectively, while continuous and dashed lines correspond to the LS40.8 and LS20.4 light curves. As before, the lines in each panel overlap, indicating that the distribution of the time-lag estimates depends mainly on m and the S/N, and not on the intrinsic CS or light curve duration.

The vertical lines in all panels indicate ν_crit, estimated as explained in the previous section. For a Gaussian distribution, the excess kurtosis and skewness should both be equal to zero. At frequencies lower than ν_crit, this is roughly the case for the distribution of the time-lag estimates when m ≳ 20. In fact, the skewness is ≈ 0 at all frequencies, although the scatter of the sample skewness increases at frequencies higher than ν_crit. This result indicates that the distribution of the time-lag estimates are symmetric around their mean. The excess kurtosis on the other hand is significantly different from zero at high frequencies, asymptotically reaching the value of −1.2. This is the excess kurtosis of the uniform distribution.

The lower panels in the same figures show the probability, p_KS(ν_k), that the time-lag distribution is Gaussian, with a mean and variance equal to the sample mean and variance of the distribution at each frequency. This probability was estimated using the Kolmogorov-Smirnov (KS) test. The dotted lines in the bottom panels show p_KS(ν_k) = 0.01. This is the typical threshold probability that would normally be considered if one wanted to reject the hypothesis of a Gaussian distribution for the time-lag estimates. The vertical line in the bottom left panels of Figs. D.18−D.22 show that this probability is higher than 0.01 at frequencies below the critical frequency, even in the case of m = 10.

As we discussed above, at frequencies ≲ν_crit the analytical time-lag error estimate differs by ≲10% from the true scatter around the mean. The reason for this effect can be explained by the properties of the mean sample coherence (as discussed in Appendix C). It is not easy to explain why the distribution of the time-lag estimates should approach a Gaussian in the same frequency range (i.e. below ν_crit), although our results indicate that this is the case.

Obviously, the time-lag distribution is not exactly Gaussian. In the case of the m = 10 and 20 averaged time-lag estimates, the excess kurtosis is certainly larger than zero at all frequencies, although not dramatically so at frequencies ≲ν_crit. Likewise, p_KS(ν_k) is rather low, but not less than 0.01 in the same frequency range. Perhaps the most interesting result for many practical applications is the fact that the ± 1 analytical error corresponds to 68% of the time-lag distribution when m ≳ 20.

10. Discussion

Our investigation was based on both analytical work and extensive numerical simulations. The simulations consisted of generating intrinsically coherent artificial light curve pairs with a prescribed model PSD and CS. The case of unevenly sampled light curves can be treated with maximum likelihood methods (e.g. Miller et al. 2010b,a; Zoghbi et al. 2013a), which we will address in a future work.

The artificial light curve characteristics (i.e. duration, sampling rate and time bin size) resembled those offered by AGN observations with present and past X-ray satellites. We assumed the same model PSD for all light curves, with a RMS amplitude of A = 0.01 and a power-law shape that smoothly “bends” from a slope of −1 to −2 after a “bend-frequency” ν_b = 2 × 10^-4 Hz. This particular shape is characteristic of AGN X-ray PSDs. The values of A = 0.01 and ν_b = 2 × 10^-4 Hz correspond to the mean values determined by fitting the X-ray PSDs of a sample of ≈ 100 nearby AGN (González-Martín & Vaughan 2012).

The phases of the model CS correspond to phase-lags that are commonly assumed between AGN X-ray light curves in different energy bands. They are divided into three main categories; constant delays of 10, 150, and 550 s, power-law time-lag spectra with an amplitude of 10^-3−10^-1 s and negative index 0.5−1.5, and time-lag spectra expected in a simple reverberation scenario where the light curves are related by a top-hat response function characterised by a centroid of t₀ = 200−2000 s, width Δ = 200−2000 s and f = 0.2. The constant delays correspond to a light-crossing time of ≈ 1−100r_g and ≈ 0.1−10r_g for a black hole of mass M_BH = 10⁶M_⊙ and 10⁷M_⊙, respectively (r_g = GM_BH/c³ ≈ 5(M_BH/ 10⁶M_⊙) s is the gravitational radius of a black hole of mass M_BH). Power-law delays with an amplitude of ≈ 10^-3 s and a negative index of ≈ 1 are those typically observed between AGN X-ray light curves in different energy bands (e.g. Papadakis et al. 2001; Emmanoulopoulos et al. 2011). In addition, modelling of observed time-lag spectra with a top-hat response function requires t₀ ≈ Δ ≈ 100−300 s (for AGN with M_BH ≈ 10⁶ M_⊙) and f ≈ 0.1−0.3 (e.g. Zoghbi et al. 2011; Emmanoulopoulos et al. 2011). Despite the fact that the light curve characteristics and model time-lag spectra that we considered are directly applicable to AGN studies, most of our results, which we summarise below, should be relevant to all time-lag studies where evenly sampled time series are considered.

10.1. Light curve sampling and binning effects

The time-lag estimates are affected by aliasing and binning of the observed time series. If the observed signal is the result of regularly recording the values of a continuous underlying processes at regular intervals Δt, then time-lags should be estimated at frequencies lower than ≈ ν_Nyq/ 5 (where ν_Nyq ≡ 1/2Δt is the Nyquist frequency). This is because aliasing has the effect of decreasing and increasing the imaginary and real parts of the cross-periodogram, respectively. As a result, their ratio (which determines the phase- and time-lag estimates) is decreased, on average, and the time-lag bias increases.

If the observed series are binned (instead of sampled) over time intervals Δt, then time-lags should be estimated at frequencies ≲ν_Nyq/ 2. This is because binning suppresses the aliasing effect on the cross-periodogram. Our results are consistent with the work of Crary et al. (1998), who reported similar findings in the context of XRB studies.

The model PSDs and CS we considered decrease with increasing frequency, and hence aliasing affects the time-lag bias more severely at increasingly higher frequencies. This is, however, expected to be true regardless of the intrinsic CS, as long as the measured signal is a (stationary) random process with finite variance (which is equal to the integrated PSD over all frequencies).

10.2. Phase-flipping effects

Phase-flipping occurs when the intrinsic phase-lag spectrum exceeds the boundaries of the interval (−π,π ]. In such cases, the phase-lag estimate will jump from π to −π (or vice-versa). This can severely alter the shape of the time-lag spectra; intrinsically constant or power-law-like time-lag spectra can show broad, prominent, “oscillatory” features. In theory, the transitions from −π to π, or vice versa, should be sharp, and hence easy to detect. However, in practice this is not the case, as the time-lag estimates will be significantly biased in the vicinity of frequencies where phase-flipping occurs. The origin of this bias can be traced to the fact that the “wings” of the phase-lag estimate’s distribution which exceed the interval (−π,π ] are folded back into the allowed range. This causes the mean of the distribution to shift towards zero. The magnitude of this bias cannot be predicted a priori, since it depends on both the intrinsic “width” of the phase-lag estimate’s distribution, as well as the unknown intrinsic phase-lag spectrum.

10.3. Smoothed vs. averaged estimates

It is customary to estimate the cross-periodogram from a single, long time series and then bin it over m neighbouring frequencies (a process called smoothing), or from m individual segments of identical duration and then bin the resulting cross-periodograms at each frequency. This process is necessary to reduce the error of the time-lag estimates, and is in fact necessary to predict their error. This requires an estimation of the coherence function itself (if we estimate the coherence from the “raw” cross-periodogram it will be equal to unity at all frequencies, irrespective of the intrinsic coherence; see P81).

Smoothing of the cross-periodogram can potentially introduce a serious bias to the time-lag estimates. For example, the model CS we considered have a power-law dependence on frequency. This can result in biased estimates of the real and imaginary part when we perform a linear type of smoothing. As a result, in some of the numerical experiments we performed, the time-lag estimates converge to zero at low (≲2 × 10^-4 Hz) frequencies. Since the time-lag bias owing to smoothing originates from the cross-periodogram bias, it can only be predicted by prescribing a model CS (and not just a model time-lag spectrum). We therefore suggest the cross-periodogram to be averaged over individual segments, rather than smoothed, before computing time-lag estimates.

10.4. Measurement error effects

In the presence of measurement errors, the coherence decreases at all frequencies. This effect becomes more severe at frequencies where the amplitude of noise variations dominates over the amplitude of the measured signal’s intrinsic variations, and the coherence tends to zero. In the case when the intrinsic coherence (i.e. the coherence in the absence of measurement errors) is unity, we found that the mean sample coherence is well fitted by Eq. (25). The parameters ν₀ and α depend only on the S/N of the light curves. The lower the S/N, the lower ν₀ and the steeper α will be. At frequencies ≳ ν_crit the sample coherence is biased, in the sense that its mean converges to the constant 1 /m, while the coherence of the noisy processes tends to zero.

Measurement errors cause the distribution of the time-lag estimates to become increasingly uniform at high frequencies. Their mean converges to zero, and their scatter is poorly approximated by standard analytical prescriptions. This is because, as the coherence decreases with increasing frequency, the “width” of the phase-lag distribution increases. When this width is sufficiently large, the wings of the distribution exceed the interval (−π,π ], and are thus folded back into the allowed range. This causes the distribution to become increasingly uniform, and its mean to shift towards zero.

We found that measurement errors have a minimal effect on the time-lag bias at frequencies where the sample coherence is ≳ 1.2/(1 + 0.2m). Furthermore, the analytical time-lag error estimate will differ from the true scatter by ≲10% in the same frequency range, as long as m ≳ 10. In addition, if m ≳ 20, the ± 1 analytical error corresponds to 68% of the time-lag distribution, and the probability distribution of the time-lag estimates approximates a Gaussian (see Sect. 9). These results hold regardless of the intrinsic PSDs and CS.

10.5. The effects of finite light curve duration

The results mentioned in Sects. 10.1–10.4 above should be applicable in all cases, i.e. they should be independent of the intrinsic PSDs and CS. However, one potentially important issue in correctly determining the intrinsic time-lag spectrum is the cross-periodogram bias. The mean cross-periodogram is not equal to the intrinsic CS, but to its convolution with the so-called Féjer kernel (see Eq. (13)). This function approaches the Dirac δ-function in the limit T → ∞, but when T is finite the bias is non-zero. Its magnitude is difficult to predict, as it depends on the shape of the unknown intrinsic CS.

Owing to the above fact, the time-lag estimates will also be biased. We found that, for the model CS we considered, the time-lag bias is positive, in the sense that the mean time-lag estimate at a given frequency has a smaller magnitude than its intrinsic value. We also found that the relative time-lag bias is ≲15% when T ≳ 20 ks, for all model time-lag spectra we considered.

In general, we found that the time-lag bias decreases with increasing light curve duration as . However, it is difficult to determine the normalisation of this relation. For a given intrinsic CS, it depends on the frequency range on which time-lags are estimated as well as the cross-periodogram bias, i.e. the shape and amplitude of the intrinsic CS. In other words, the naive expectation of the relative time-lag bias being negligible if the light curve duration is much larger than the magnitude of the time-lag estimates is generally incorrect. For example, in experiment CD1, where the intrinsic time-lag spectrum is equal to 10 s at all frequencies, the relative time-lag bias is ≈ 20% at the lowest sampled frequency, even if we use ≈ 100 ks light curves (see e.g. the continuous black curve curve in the top left panel of Fig. D.1).

To illustrate the complexity of the time-lag bias’ dependence on the intrinsic CS, we revisited experiment CD1. The solid line in Fig. 12 shows the relative time-lag bias determined from ≈ 100 ks light curves (this line is identical to the black line in the top left panel of Fig. D.1). The dashed red line shows the expected relative time-lag bias for the same light curves and intrinsic time-lag spectrum, and an intrinsic PSD whose amplitude is larger by a factor of five. We used Eq. (13) to estimate the mean cross-periodogram, E [ I_xy(ν_p) ], and then determined the expected relative time-lag bias by computing { 10 s−(2πν_p)^-1arg { E [ I_xy(ν_p) ] } } /10 s at each frequency. The relative bias does not change when compared to the original PSD amplitude. We then repeated the same calculation, but assumed a “bend-frequency” of 2 × 10^-5 Hz (i.e. 10 times lower than the original value). The expected relative bias (shown by the blue dashed-dotted line in the same figure) is now increased at each frequency, reaching a value of ≈ 30% at the lowest frequency.

These results can be understood by inspecting Eq. (13), which determines the cross-periodogram bias. In our simulations, the CS is given by h_XY(ν) = (μ_Y/μ_X)h_X(ν)e^iφ_XY(ν), where h_X(ν) and φ_XY(ν) are the intrinsic PSD and phase-lag spectrum, respectively. Since the intrinsic PSD is given by Eq. (14), a change in the PSD amplitude has the effect of altering the amplitude of the imaginary and real parts of the intrinsic CS by the same multiplicative factor at each frequency. This however does not change their ratio, which determines the mean of the phase- and time-lag estimates, and hence leaves the relative time-lag bias unchanged. On the other hand, a change in the bend-frequency alters the amplitude of the imaginary and real parts of the intrinsic CS differently at each frequency, resulting in a change of the relative time-lag bias.

We therefore suggest that it is important to always check the expected relative time-lag bias to validate the time-lag analysis in practice. This can be done by assuming a possible model CS using Eq. (13). If the intrinsic coherence of the light curves is close to unity, such a model CS may be determined by modelling the PSDs and considering various model time-lag spectra that may be potential candidates for the intrinsic time-lag spectrum.

11. Summary

We investigated the statistical properties of Fourier-based time-lag estimates. These estimates are based on the cross-periodogram, which is an estimator of the intrinsic CS between two random time series. Unlike the periodogram (the traditional PSD estimator) which is always positive-definite, the cross-periodogram is generally a complex number whose phase is an estimator of the phase-lag spectrum (the time-lag is equal to the phase-lag divided by the angular frequency). Our aims were to quantify the effects of light curve characteristics and measurement errors on the time-lag bias (i.e. the difference between their mean and intrinsic values), to study their distribution, and to investigate whether the standard analytical error prescriptions are accurate approximations of their intrinsic scatter around the mean.

Based on our results, we suggest the following steps when estimating time-lags:

• Estimate the averaged cross-periodogram (Eq. (18)), using at least m = 20 pairs of data segments.

• Use the averaged cross-periodogram to estimate the coherence and time-lags (Eqs. (23) and (20)), along with their errors (Eqs. (24) and (22)).

• Determine the frequency range for which the sample coherence is greater than 1.2/(1 + 0.2m). If the intrinsic coherence is close to unity⁶, then the coherence estimates should be well fitted by Eq. (25). By determining the best-fit ν₀ and α values of this function, the aforementioned frequency range can then be determined by solving (1 − m^-1)exp [ −(ν/ν₀)^α ] + m^-1 ≥ 1.2/(1 + 0.2m).

At the frequency range where the sample coherence is larger than 1.2/(1 + 0.2m), the bias of the averaged time-lag estimates will be similar to its value in the absence of measurement errors, the analytical prescription for the error of the time-lag estimates will be similar to the true scatter around the mean, and their distribution will be roughly approximated by a Gaussian, in the sense discussed in Sect. 9.

---

Acknowledgments

We thank the referee for his/her suggestions, which significantly improved the quality and clarity of the manuscript. This work was supported by the AGNQUEST project, which is implemented under the Aristeia II action of the Education and Lifelong Learning operational programme of the GSRT, Greece. It was also supported in part by the grant PIRSES-GA-2012-31578 EuroCal.

References

Appendix A: The cross-spectrum of two binned time series

Let us consider the continuous process { X(t),Y(t) } and the discrete process { X(t_s),Y(t_s) }, defined according to Eq. (7) in Sect. 3. It is straight-forward to show that the processes have the same mean values. However, this is not true for their CCFs and CS as well.

For notational convenience and without loss of generality, we will henceforth assume that the intrinsic process has a mean value of zero. By definition, the CCF of the discrete process is only defined at lags which are integer multiples of the sampling period, i.e. τ_l = lΔt_sam, where l = 0, ± 1, ± 2,..., and is given by (A.1)where R_XY(τ) ≡ E [ Y(t)X(t + τ) ] is, by definition, the CCF of the continuous process. According to Eq. (A.1) R_XY(τ_l) and R_XY(τ_l) are not equal at any value of τ_l. Given that , Eq. (A.1) becomes (A.2)Splitting the infinite frequency integration interval into an infinite number of segments with width 1/2Δt_sam, Eq. (A.2) becomes (A.3)since e^i2πkl = 1 for all k, l. The relation between the CCF and CS of the discrete process, h_XY(ν) (which is only defined on the interval | ν | ≤ 1/2Δt_sam), is (A.4)Comparing Eqs. (A.3) and (A.4), we arrive at the desired result: (A.5)

Appendix B: The mean of the cross-periodogram

Using Eqs. (9) and (10), we may write the cross-periodogram of the time series { x(t_r),y(t_r) } as follows: (B.1)We can now transform the indices r and k to s = r − k and l, where s goes from {−(N − 1) } to (N − 1), and, at fixed s, l goes from 1 to (N − | s |). With this transformation, Eq. (B.1) becomes (B.2)where is, by definition, CCF, , and τ_s = t_j − t_k = sΔt_sam is the lag. Therefore, Eq. (B.2) can be re-cast into the following form: (B.3)Applying the expectation operator on both sides of Eq. (B.3), we get (B.4)where we have ignored the effect of estimating the mean values of the observed time series to reach the last equation. R_XY(s) is the CCF of the discrete processes { X(t_s),Y(t_s) } that was defined in Sect. 3. Since the CCF is equal to the inverse Fourier transform of the CS, Eq. (B.4) becomes (B.5)where (B.6)is the so-called Fejér kernel. One of the properties of this function is that, as N → ∞, F_N(ν) → δ(ν) (the Dirac δ-function). Thus, as N → ∞ the right-hand side of Eq. (B.5) converges to the intrinsic CS (modified by the effects of discrete binning or sampling according to Eq. (A.5)) value at ν, h_XY(ν). Therefore, the cross-periodogram is an asymptotically unbiased estimator of the CS.

Appendix C: Effects of measurement errors on the coherence

Let us consider an intrinsically continuous random process { X(t),Y(t) }, whose mean, PSDs and CS are { μ_X,μ_Y }, { h_X(ν),h_Y(ν) } and h_XY(ν), respectively. Let us also assume that this process is continuously sampled at regular time intervals t_s = sΔt_sam, where s = 0, ± 1, ± 2,..., and binned over time bins of size Δt_bin. In the absence of measurement errors, the measured signal corresponds to a single realisation of a discrete version of the intrinsic process, { X(t_s),Y(t_s) }, as defined by Eq. (7). If there are measurement errors (as is always the case), the observed time series constitute a realisation of the “noisy” process (C.1)where ϵ_X and ϵ_Y are a “purely random” processes with zero mean, and variance and , respectively. In this case, where and are the PSDs of the purely random processes. Equations (C.2) and (C.3) describe the well-known fact that the addition of so-called white noise to a process has the effect of adding a constant amount of “power” to its intrinsic PSD. When the measurement errors in one light curve are uncorrelated with the errors in the other, then h_XY,n(ν) = h_XY(ν). As a result, (C.4)where and is the coherence of the noisy and the “intrinsic” process, respectively (here intrinsic refers to the discrete process { X(t_s),Y(t_s) }). According to Eq. (C.4), the presence of (intrinsically uncorrelated) measurement errors decreases the coherence of the processes under study. In fact, at high frequencies where the intrinsic PSDs are much smaller than h_ϵX and h_ϵY, the ratio will be equal to h_X(ν)h_Y(ν)/(h_ϵXh_ϵX). In the case of power-law like PSDs, which decrease with increasing frequency, as is typically the case in AGN, this ratio will then tend to zero at high frequencies.

In addition to decreasing the intrinsic coherence, measurement errors will introduce a bias in the sample coherence as well, i.e. the mean of sample coherence will not be equal to . This can be shown by following (Vaughan & Nowak 1997), who showed that (C.5)where is the CS estimate in the absence of measurement errors, and is a random complex number with zero mean and variance given by (C.6)Since is a random complex number, it may be represented by a vector in the complex plane with a random phase distributed uniformly over the interval (−π,π ]. As such, randomly “perturbs” both the magnitude as well as the direction of the vector . On average this perturbation will have no net effect on either the magnitude or direction of ĥ_xy(ν_k), hence , although it will obviously increase its variance. This implies that measurement errors will not bias the cross-periodogram (although they will bias the phase-lag estimates at high frequencies, as discussed in Sect. 7.1).

The mean sample coherence can now be determined. By first applying the expectation operator on both sides of Eq. (23), we get (C.7)Given that E [ ĥ_x(ν_k) ] ~ h_X,n(ν_k) and E [ ĥ_y(ν_k) ] ~ h_Y,n(ν_k), if we substitute Eqs. (C.2)−(C.4) and (C.6) into Eq. (C.7), we get (C.8)where is defined by Eq. (C.4). The above equation holds true for any intrinsic coherence. In our case where at all frequencies, Eq. (C.8) reduces to (C.9)This relation is identical to Eq. (25) when . We found that such a relation fits the mean sample coherence well in the case of experiments CD1, PLD2, and THRF1, and both light curve types we considered in Sects. 7–9. A comparison between this relation and the mean sample coherence is shown in Fig. 11 for the LS40.8 light curves in experiment THRF1. Equation (C.9) indicates that when m is large, the mean sample coherence will be equal to the coherence of the noisy process.

However, in general, at sufficiently high frequencies the ratio tends to zero because tends to zero and is always smaller than at all frequencies (Eq. (C.4)). In this case, the first term in the right-hand side of Eq. (C.8) will therefore converge to zero. Consequently, (C.10)

Equation (C.10) shows that, in the frequency range where the experimental noise variations dominate over the intrinsic ones (i.e. when ), the mean sample coherence will be equal to 1 /m.

Appendix D: Figures that show our results

All Figures