Statistical properties of Fourierbased timelag estimates
^{1}
Department of Physics and Institute of Theoretical and Computational
Physics, University of Crete, 71003
Heraklion, Greece
email:
epitrop@physics.uoc.gr
^{2}
IESL, Foundation for Research and TechnologyHellas,
71110
Heraklion, Crete, Greece
Received: 30 October 2015
Accepted: 2 April 2016
Context. The study of Xray timelag spectra in active galactic nuclei (AGN) is currently an active research area, since it has the potential to illuminate the physics and geometry of the innermost region (i.e. close to the putative supermassive black hole) in these objects. To obtain reliable information from these studies, the statistical properties of timelags estimated from data must be known as accurately as possible.
Aims. We investigated the statistical properties of Fourierbased timelag estimates (i.e. based on the crossperiodogram), using evenly sampled time series with no missing points. Our aim is to provide practical “guidelines” on estimating timelags that are minimally biased (i.e. whose mean is close to their intrinsic value) and have known errors.
Methods. Our investigation is based on both analytical work and extensive numerical simulations. The latter consisted of generating artificial time series with various signaltonoise ratios and sampling patterns/durations similar to those offered by AGN observations with present and past Xray satellites. We also considered a range of different model timelag spectra commonly assumed in Xray analyses of compact accreting systems.
Results. Discrete sampling, binning and finite light curve duration cause the mean of the timelag estimates to have a smaller magnitude than their intrinsic values. Smoothing (i.e. binning over consecutive frequencies) of the crossperiodogram can add extra bias at low frequencies. The use of light curves with low signaltonoise ratio reduces the intrinsic coherence, and can introduce a bias to the sample coherence, timelag estimates, and their predicted error.
Conclusions. Our results have direct implications for Xray timelag studies in AGN, but can also be applied to similar studies in other research fields. We find that: a) timelags should be estimated at frequencies lower than ≈ 1/2 the Nyquist frequency to minimise the effects of discrete binning of the observed time series; b) smoothing of the crossperiodogram should be avoided, as this may introduce significant bias to the timelag estimates, which can be taken into account by assuming a model crossspectrum (and not just a model timelag spectrum); c) timelags should be estimated by dividing observed time series into a number, say m, of shorter data segments and averaging the resulting crossperiodograms; d) if the data segments have a duration ≳ 20 ks, the timelag bias is ≲15% of its intrinsic value for the model crossspectra and powerspectra considered in this work. This bias should be estimated in practise (by considering possible intrinsic crossspectra that may be applicable to the timelag spectra at hand) to assess the reliability of any timelag analysis; e) the effects of experimental noise can be minimised by only estimating timelags in the frequency range where the sample coherence is larger than 1.2/(1 + 0.2m). In this range, the amplitude of noise variations caused by measurement errors is smaller than the amplitude of the signal’s intrinsic variations. As long as m ≳ 20, timelags estimated by averaging over individual data segments have analytical error estimates that are within 95% of the true scatter around their mean, and their distribution is similar, albeit not identical, to a Gaussian.
Key words: methods: statistical
© ESO, 2016
1. Introduction
The study of timelags as a function of temporal frequency between Xray 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 Xray binaries (XRBs; e.g. Miyamoto & Kitamoto 1989; Nowak & Vaughan 1996; Nowak et al. 1999). In the last few years, timelag studies have revealed that “soft” Xray variations in AGN are delayed with respect to “hard” Xray 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 timelags 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 Xray reverberation. Further credence to such a scenario came with recent discoveries of timelags between the Fe Kα emission line (≈ 5−7 keV) and the Xray 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 socalled Compton hump (≈ 10−30 keV) and the Xray continuum (e.g. Zoghbi et al. 2014; Kara et al. 2015) in a few AGN.
The study of Xray reverberation timelags can provide valuable geometrical and physical information on the Xray 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 timelags 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 timelag spectrum), error, and probability distribution. The later is necessary if one wishes to fit the observed timelag 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 Fourierbased timelag estimates, based on both analytical work and extensive numerical simulations. Our primary goals are: a) to investigate whether the frequently used timelag estimates are indeed reliable estimates of the intrinsic timelag 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 χ^{2} minimisation techniques. Our work should have direct impact on current timelag studies in the area of AGN Xray 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 timevarying 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 autocovariance functions (R_{X}(τ) and R_{Y}(τ); τ is the socalled lag) of the individual time series are finite and timeindependent (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 crosscovariance 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 crossspectrum (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 phaselag spectrum, which is defined as (4)The phaselag 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 timelag 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 timelag 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 socalled coherence function, (6)where h_{X}(ν) and h_{Y}(ν) denote the powerspectral density functions (PSDs) of the time series X(t) and Y(t), respectively (the PSD is defined as the Fourier transform of the autocovariance 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^{1} in astronomy are the result of binning the intrinsic signal over time bins of size Δt_{bin}^{2}. 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 straightforward 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 socalled 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^{2} function, whose argument depends on the bin size, Δt_{bin}. Since the sinc^{2} 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^{2} term in the righthand side of Eq. (8) vanishes, and hence the effects of aliasing on h_{XY}(ν) are maximised.
3. Estimation of crossspectra; the crossperiodogram
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 socalled 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 crossperiodogram of two time series is defined as (10)where the asterisk ^{∗} denotes complex conjugation. The crossperiodogram 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 timelag spectrum, respectively.
3.1. The bias of crossspectral 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 righthand side of Eq. (13) converges to h_{XY}(ν_{p}). The crossperiodogram 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 crossperiodogram 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 timelag estimates based on the crossperiodogram, 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 Xray 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)^{3} 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 countrate, μ_{X}. We assumed a model PSD of the form (14)This function describes a power law with a lowfrequency slope of −1 which smoothly “bends” to a slope of −2 at frequencies above the bendfrequency, ν_{b}, as in the case of inferred Xray 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^{6} and Δt_{sam} = 1 s. By construction, the intrinsic PSD of the generated light curves is discrete, and has nonzero 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 countrate μ_{Y}, assuming a specified model phaselag 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, nonzero only at the same frequencies where the PSD is nonzero as well. The light curve pairs generated following this procedure will henceforth be referred to as original realisations.
4.1. The model timelag spectra
We considered three different model timelag 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 powerlaw index of the corresponding model timelag 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.
Tophat response functions: in this case we assumed that φ_{XY}(ν) = arg [ 1 + fe^{i2πνt0}sinc(πνΔ) ]. This phaselag spectrum is expected if the intrinsic time series are related by the following equation: , where Ψ(t) (the socalled response function) is a simple socalled tophat function, i.e. it has a constant value of f/ Δ in the interval  t − t_{0}  ≤ Δ/2 and zero otherwise (Δ is the width of the tophat and t_{0} its centroid). We considered the model parameter values { f,t_{0},Δ } = { 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 timelag 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 Xray astronomy are provided by satellites in lowEarth 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 timescales. 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 lowEarth orbit satellites, XMMNewton 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 crossperiodogram and timelag estimates directly, or “chop” them into segments of shorter duration, calculate the crossperiodogram for each segment and then bin the resulting estimates at certain frequencies to estimate the timelags (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 timelags 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 timelag estimates, (henceforth, the timelag bias).
5. The bias of the timelag 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 Xray timing studies use this estimator. In addition, we computed the quantity to quantify the timelag bias in terms of its intrinsic value (henceforth, the relative bias). This allows us to directly compare the timelag 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 timelag spectra is that they decrease to zero at high frequencies, irrespective of the light curve sampling pattern, duration, and model timelag spectrum. This decrease is due to the effects of aliasing and lightcurve 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.42 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.43 light curves). Since the original realisations in this case have the same sampling rate as the LS102.43 light curves, they should not be affected by either binning or aliasing.
Fig. 1 Sample mean of the real and imaginary parts of the crossperiodogram (top left and right panels, respectively), the sample mean timelag spectrum and the relative timelag bias (bottom left and right panel, respectively) in experiment PLD2. The vertical black and red dashed lines indicate ν_{Nyq}/ 2 and ν_{Nyq}/ 5, respectively. Above these frequencies, the LS102.4 (black curve) and LS102.42 (red curve) relative timelag bias begins to noticeably increase (see text for more details). The horizontal dotted line in the bottom right panel, and in all subsequent plots, indicates the 0.1 (i.e. 10%) relative timelag bias. 

Open with DEXTER 
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.42, and LS102.43 light curves (black, red, and brown curves, respectively)^{4}. The sample mean timelag spectra for the LS102.42 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.43 light curves, as shown by the brown curves in the bottom panels of Fig. 1. The timelag 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 crossperiodogram 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 phaselag estimate) is decreased on average, hence the timelag bias increases. This increase is most severe at high frequencies, where the effects of aliasing and light curve binning on the crossperiodogram are maximised. The timelag bias owing to these effects is more pronounced for the LS102.42 light curves, since light curve binning suppresses the aliasing effect (see Sect. 3).
Fig. 2 Mean relative timelag bias over all frequencies below ν_{max}, plotted as a function of ν_{max}, for different light curve types in various numerical experiments. The left and right dashed vertical lines indicate ν_{Nyq}/ 2 for the OB and LStype light curves, respectively. 

Open with DEXTER 
Figures D.1–D.4 show that the onset of the highfrequency increase in the relative timelag bias occurs at ≈ ν_{Nyq}/ 2 for all light curve types, model timelag spectra, and light curve durations. To further illustrate this effect, we define the function as the mean relative timelag 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 timelag bias remains approximately constant, but increases markedly at higher frequencies (the same trend is observed in all numerical experiments we considered).
We conclude that timelags should be estimated at frequencies ≲ν_{Nyq}/ 2 to minimise the effects of lightcurve binning and aliasing on the timelag bias.
5.2. Finite light curve duration
The timelag estimates are biased even at frequencies ≲ν_{Nyq}/ 2. Figures D.1–D.4 show that the relative timelag 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 socalled phaseflipping (see Sect. 5.3 for details). The dashed red line indicates the bestfit relation to the LStype data: . This relation fits the LStype data quite well. We therefore conclude that the mean relative timelag 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 LStype light curves. This discrepancy arises because the frequency range probed by the OB light curves is significantly lower than by the LStype light curves (owing to their longer duration and larger time bin size), as well as the fact that the relative timelag 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 timelag 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 OBtype 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 OBtype 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 timelag bias depends on , although the normalisation of this relation depends on the frequency range on which the timelags are estimated; the lower the frequencies, the more biased the timelag estimates will be. For the model CS we considered, our results suggest that, to obtain timelag estimates with a relative bias ≲15%, the light curves must have T ≳ 20 ks.
Fig. 3 Mean relative timelag bias, plotted as function of light curve duration for LStype and OB light curves. 

Open with DEXTER 
Fig. 4 Relative timelag bias for the LS and OBtype light curves (top and bottom panel, respectively), for various durations in experiment CD1. Filled squares indicate the mean relative timelag bias over the full sampled frequency range, evaluated at the mean logarithmic frequency. The dashed vertical line in the top and bottom panel indicate ν_{Nyq}/ 2 for the LS and OBtype light curves, respectively. 

Open with DEXTER 
5.3. Phaseflipping
Since the phaselag 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 viceversa). This effect is known as phaseflipping (or phasewrapping).
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 timelag spectrum (black dashed line). The simulated light curves in this experiment are separated by a constant delay. As a result, the model phaselag 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 timelag spectrum will likewise jump from 550 s to −550 s. Subsequent phaseflips occur at frequencies ν = j/ 2d (j = 2,3,...), where the model timelag spectrum undergoes discontinuous jumps from positive to negative values with decreasing amplitude.
The mean sample timelag 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 phaselag 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 phaseflip occurs, while the other two are lower and higher than that frequency. The solid and dashed vertical lines in the topright and bottom panels of Fig. 5 show the model and sample mean phaselag value at these frequencies, respectively. Since the phaselag estimate is defined on the interval (−π,π ], the socalled 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 phaselag spectrum. This effect is most severe in the vicinity of frequencies where phaseflipping occurs, since, in this case, the mean of the phaselag estimate is close to zero.
Fig. 5 Top left panel: mean sample timelag spectrum for the LS102.4 light curves in experiment CD3 (continuous black line). The black dashed line indicates the model timelag spectrum. Top right and bottom panels: the probability distribution of the phaselag estimates at the frequencies indicated by the vertical dashed lines in the topleft panel. The solid and dashed vertical line in these panels indicate the model and sample mean phaselag value at these frequencies, respectively. 

Open with DEXTER 
Phaseflipping does not necessarily occur at increasingly higher frequencies. If, as in the case of experiments PLD4 and PLD5 for example, the intrinsic phaselag spectrum increases (in absolute magnitude) with decreasing frequency, phaseflipping 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 timelag spectrum for this experiment is τ_{XY}(ν) = 0.01ν^{1.5} (indicated by the continuous brown line). The corresponding phaselag spectrum is φ_{XY}(ν) = (2πν)0.01ν^{1.5}, and phaseflipping occurs at (4 × 10^{4}/j^{2}) Hz (j = 1,2,...). At frequencies close to 4 × 10^{4} Hz, where the first phaseflip occurs, the mean sample timelag spectrum exhibits a jump which is smoother than the abrupt jump of the model timelag 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 (bottomright 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 crossperiodogram (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 timelag 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}(ν) ] }.
Fig. 6 As in Fig. 1, for experiment PLD5. The blue dashed lines indicate the model CS (upper panels) and timelag spectrum (lower left panel). The solid brown line in the bottom left panel indicates the model timelag spectrum without taking the effects of phaseflipping into account. 

Open with DEXTER 
The timelag estimates at the two lowest frequencies are heavily biased owing to the bias of the real and imaginary parts of the crossperiodogram. This is reflected in the mean sample timelag 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 crossperiodogram to diverge from the model CS because of its rapid oscillatory behaviour at low frequencies.
6. Smoothed/averaged crossspectral estimators
The variance of the real and imaginary parts of the crossperiodogram 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 crossperiodogram 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^{5}. 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 crossperiodogram for each segment, (l = 1,2,...,m, and ν_{k} = k/ (T/m), where k = 1,2,...,N/m), and then average the different crossperiodogram 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 crossperiodogram (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 timelag 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 timelag estimator. The lower the coherence, the higher the variance of the phase and timelag 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 crossspectral estimates, we smoothed the real and imaginary parts of the crossperiodograms 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 crossperiodograms 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 timelag 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 phaseflipping in the sampled frequency range). The bottom panel shows the same plot for the averaged estimates. The blue dashed line shows the onetoone relation in each case. Although the bias of the averaged estimates falls on the expected onetoone 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 nonlinear variation with frequency. As a result, linear smoothing introduces a bias to the real and imaginary parts of the smoothed crossperiodogram. 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 timelag 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 crossperiodogram and timelag estimate is maximised.
Fig. 7 Mean relative bias of the m = 20 smoothed (top panel) and averaged (bottom panel) timelag estimates, plotted as a function of the relative bias of the nonsmoothed/averaged timelag estimates (the dashed lines show the onetoone relation line) in various numerical experiments. 

Open with DEXTER 
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 crossperiodogram 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 timelag 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 phaseflipping effects, and are representative of the three categories of model phaselag spectra we consider. In total, we thus ended up with 5 × 10^{4} LS20.4 and 2 × 10^{4} LS40.8 light curve pairs for each experiment.
Fig. 8 As in Figs. 1 and 6, for the m = 20 smoothed estimates obtained from the LS102.4 light curves in experiment CD1. The blue dashed lines indicate the model CS (upper panels) and timelag spectrum (lower left panel). 

Open with DEXTER 
Furthermore, we created five copies of each of these pairs to simulate the effects of different signaltonoise 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 crossperiodogram to estimate the timelags and coherence according to Eqs. (20) and (23), along with their error estimates according to Eqs. (22) and (24). The number of averaged timelag 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 timelag bias
Fig. 9 Relative bias of the averaged timelag estimates for the LS20.4 and LS40.8 light curves, plotted as a function of (S/N)_{xy} in experiment THRF1. The horizontal dashed lines correspond to the value of the relative bias when (S/N)_{xy} → ∞. 

Open with DEXTER 
A comparison between, for example, the top panels in Fig. D.10 and the bottom left panel in Fig. D.5 reveals that the timelag 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 timelag 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 dotteddashed 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 timelag 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 timelag 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 timelags in the top panels in Fig. D.10 and in the bottom left panel of Fig. D.5). As in the case of phaseflipping (see Sect. 5.3), this bias arises because the phaselag estimates are defined on the interval (−π,π ]. As an example, in Fig. 10 we show the probability distribution of phaselag 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 phaselag estimates, respectively. As the frequency increases, the mean sample coherence decreases (see the next section), and hence the scatter of the phaselag 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 phaselag 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 phaselag 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 dotteddashed 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 timelag bias is similar to its value in the absence of measurement errors at frequencies below ν_{crit}. The same result holds for the timelag estimates in the case of m = 10, 30, and 40 as well.
Fig. 10 Probability distribution of the averaged phaselag estimate for the LS40.8 light curves in experiment THRF1 at different frequencies. 

Open with DEXTER 
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 bestfit 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).
Fig. 11 Mean sample coherence (m = 20) obtained from the LS40.8 light curves in experiment THRF1 for different S/N combinations. The blue horizontal dashed lines indicate the value 1 /m (see Sect. 7.2 for details). 

Open with DEXTER 
8. The error of the timelag 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 socalled 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 timelag estimates around their mean. The bottom panels in the same figures show the probability , i.e. the probability that the value of the timelag 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 timelag 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 timelag 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 timelag 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 dotteddashed 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 timelag 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 timelag 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 timelag estimates develops longer “tails” than would be expected from a Gaussian distribution. We investigate the properties of the probability distribution of the timelag estimates in more detail below.
9. The probability distribution of the timelag estimates
The results presented above regarding the error of the timelag 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 timelag 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 timelag 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 timelag 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 timelag 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 timelag 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 KolmogorovSmirnov (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 timelag 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 timelag 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 timelag 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 timelag distribution is not exactly Gaussian. In the case of the m = 10 and 20 averaged timelag 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 timelag 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 Xray satellites. We assumed the same model PSD for all light curves, with a RMS amplitude of A = 0.01 and a powerlaw shape that smoothly “bends” from a slope of −1 to −2 after a “bendfrequency” ν_{b} = 2 × 10^{4} Hz. This particular shape is characteristic of AGN Xray PSDs. The values of A = 0.01 and ν_{b} = 2 × 10^{4} Hz correspond to the mean values determined by fitting the Xray PSDs of a sample of ≈ 100 nearby AGN (GonzálezMartín & Vaughan 2012).
The phases of the model CS correspond to phaselags that are commonly assumed between AGN Xray light curves in different energy bands. They are divided into three main categories; constant delays of 10, 150, and 550 s, powerlaw timelag spectra with an amplitude of 10^{3}−10^{1} s and negative index 0.5−1.5, and timelag spectra expected in a simple reverberation scenario where the light curves are related by a tophat response function characterised by a centroid of t_{0} = 200−2000 s, width Δ = 200−2000 s and f = 0.2. The constant delays correspond to a lightcrossing time of ≈ 1−100r_{g} and ≈ 0.1−10r_{g} for a black hole of mass M_{BH} = 10^{6}M_{⊙} and 10^{7}M_{⊙}, respectively (r_{g} = GM_{BH}/c^{3} ≈ 5(M_{BH}/ 10^{6}M_{⊙}) s is the gravitational radius of a black hole of mass M_{BH}). Powerlaw delays with an amplitude of ≈ 10^{3} s and a negative index of ≈ 1 are those typically observed between AGN Xray light curves in different energy bands (e.g. Papadakis et al. 2001; Emmanoulopoulos et al. 2011). In addition, modelling of observed timelag spectra with a tophat response function requires t_{0} ≈ Δ ≈ 100−300 s (for AGN with M_{BH} ≈ 10^{6} 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 timelag spectra that we considered are directly applicable to AGN studies, most of our results, which we summarise below, should be relevant to all timelag studies where evenly sampled time series are considered.
10.1. Light curve sampling and binning effects
The timelag 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 timelags 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 crossperiodogram, respectively. As a result, their ratio (which determines the phase and timelag estimates) is decreased, on average, and the timelag bias increases.
If the observed series are binned (instead of sampled) over time intervals Δt, then timelags should be estimated at frequencies ≲ν_{Nyq}/ 2. This is because binning suppresses the aliasing effect on the crossperiodogram. 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 timelag 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. Phaseflipping effects
Phaseflipping occurs when the intrinsic phaselag spectrum exceeds the boundaries of the interval (−π,π ]. In such cases, the phaselag estimate will jump from π to −π (or viceversa). This can severely alter the shape of the timelag spectra; intrinsically constant or powerlawlike timelag 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 timelag estimates will be significantly biased in the vicinity of frequencies where phaseflipping occurs. The origin of this bias can be traced to the fact that the “wings” of the phaselag 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 phaselag estimate’s distribution, as well as the unknown intrinsic phaselag spectrum.
10.3. Smoothed vs. averaged estimates
It is customary to estimate the crossperiodogram 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 crossperiodograms at each frequency. This process is necessary to reduce the error of the timelag 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” crossperiodogram it will be equal to unity at all frequencies, irrespective of the intrinsic coherence; see P81).
Smoothing of the crossperiodogram can potentially introduce a serious bias to the timelag estimates. For example, the model CS we considered have a powerlaw 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 timelag estimates converge to zero at low (≲2 × 10^{4} Hz) frequencies. Since the timelag bias owing to smoothing originates from the crossperiodogram bias, it can only be predicted by prescribing a model CS (and not just a model timelag spectrum). We therefore suggest the crossperiodogram to be averaged over individual segments, rather than smoothed, before computing timelag 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 ν_{0} and α depend only on the S/N of the light curves. The lower the S/N, the lower ν_{0} 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 timelag 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 phaselag 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 timelag bias at frequencies where the sample coherence is ≳ 1.2/(1 + 0.2m). Furthermore, the analytical timelag 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 timelag distribution, and the probability distribution of the timelag 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 timelag spectrum is the crossperiodogram bias. The mean crossperiodogram is not equal to the intrinsic CS, but to its convolution with the socalled Féjer kernel (see Eq. (13)). This function approaches the Dirac δfunction in the limit T → ∞, but when T is finite the bias is nonzero. Its magnitude is difficult to predict, as it depends on the shape of the unknown intrinsic CS.
Owing to the above fact, the timelag estimates will also be biased. We found that, for the model CS we considered, the timelag bias is positive, in the sense that the mean timelag estimate at a given frequency has a smaller magnitude than its intrinsic value. We also found that the relative timelag bias is ≲15% when T ≳ 20 ks, for all model timelag spectra we considered.
In general, we found that the timelag 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 timelags are estimated as well as the crossperiodogram bias, i.e. the shape and amplitude of the intrinsic CS. In other words, the naive expectation of the relative timelag bias being negligible if the light curve duration is much larger than the magnitude of the timelag estimates is generally incorrect. For example, in experiment CD1, where the intrinsic timelag spectrum is equal to 10 s at all frequencies, the relative timelag 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 timelag bias’ dependence on the intrinsic CS, we revisited experiment CD1. The solid line in Fig. 12 shows the relative timelag 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 timelag bias for the same light curves and intrinsic timelag spectrum, and an intrinsic PSD whose amplitude is larger by a factor of five. We used Eq. (13) to estimate the mean crossperiodogram, E [ I_{xy}(ν_{p}) ], and then determined the expected relative timelag bias by computing { 10 s−(2πν_{p})^{1}arg { 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 “bendfrequency” of 2 × 10^{5} Hz (i.e. 10 times lower than the original value). The expected relative bias (shown by the blue dasheddotted 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 crossperiodogram 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 phaselag 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 timelag estimates, and hence leaves the relative timelag bias unchanged. On the other hand, a change in the bendfrequency alters the amplitude of the imaginary and real parts of the intrinsic CS differently at each frequency, resulting in a change of the relative timelag bias.
We therefore suggest that it is important to always check the expected relative timelag bias to validate the timelag 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 timelag spectra that may be potential candidates for the intrinsic timelag spectrum.
Fig. 12 Relative timelag bias for the LS102.4 light curves and different values of the model PSD parameters, { A,ν_{b} }, in experiment CD1. 

Open with DEXTER 
11. Summary
We investigated the statistical properties of Fourierbased timelag estimates. These estimates are based on the crossperiodogram, which is an estimator of the intrinsic CS between two random time series. Unlike the periodogram (the traditional PSD estimator) which is always positivedefinite, the crossperiodogram is generally a complex number whose phase is an estimator of the phaselag spectrum (the timelag is equal to the phaselag divided by the angular frequency). Our aims were to quantify the effects of light curve characteristics and measurement errors on the timelag 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 timelags:

Estimate the averaged crossperiodogram (Eq. (18)), using at least m = 20 pairs of data segments.

Use the averaged crossperiodogram to estimate the coherence and timelags (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^{6}, then the coherence estimates should be well fitted by Eq. (25). By determining the bestfit ν_{0} and α values of this function, the aforementioned frequency range can then be determined by solving (1 − m^{1})exp [ −(ν/ν_{0})^{α} ] + 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 timelag estimates will be similar to its value in the absence of measurement errors, the analytical prescription for the error of the timelag 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.
We note that Emmanoulopoulos et al. (2013) have more recently presented an improved method for producing artificial light curves. Since we are not interested in the probability distribution of the synthetic time series in our study, we used the Timmer & Koenig method to minimise the time needed to perform the simulations.
The sample mean of the real and imaginary parts of the crossperiodograms shown in Fig. 1, and all similar subsequent figures, have been normalised by the factor μ_{X}μ_{Y} (which is the common normalisation factor of the intrinsic CS used in our numerical experiments; see Sect. 4) to make them dimensionless and independent of the specific countrates assumed to construct the simulated light curve pairs.
There are many smoothing schemes for the crossperiodogram. The one defined by Eq. (15) corresponds to a particular type of smoothing “window” called the Daniell, or rectangular, spectral window (see P81 for more details).
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 PIRSESGA201231578 EuroCal.
References
 Alston, W. N., Vaughan, S., & Uttley, P. 2013, MNRAS, 435, 1511 [NASA ADS] [CrossRef] [Google Scholar]
 Arévalo, P., Papadakis, I. E., Uttley, P., McHardy, I. M., & Brinkmann, W. 2006, MNRAS, 372, 401 [NASA ADS] [CrossRef] [Google Scholar]
 Arévalo, P., McHardy, I. M., & Summons, D. P. 2008, MNRAS, 388, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Bendat, J., & Piersol, A. 2011, Random Data: Analysis and Measurement Procedures (New York: Wiley) [Google Scholar]
 Crary, D. J., Finger, M. H., Kouveliotou, C., et al. 1998, ApJ, 493, L71 [NASA ADS] [CrossRef] [Google Scholar]
 De Marco, B., Ponti, G., Cappi, M., et al. 2013, MNRAS, 431, 2441 [NASA ADS] [CrossRef] [Google Scholar]
 Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2011, MNRAS, 416, L94 [NASA ADS] [CrossRef] [Google Scholar]
 Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Fabian, A. C., Zoghbi, A., Ross, R. R., et al. 2009, Nature, 459, 540 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 GonzálezMartín, O., & Vaughan, S. 2012, A&A, 544, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kara, E., Fabian, A. C., Cackett, E. M., Miniutti, G., & Uttley, P. 2013a, MNRAS, 430, 1408 [NASA ADS] [CrossRef] [Google Scholar]
 Kara, E., Fabian, A. C., Cackett, E. M., et al. 2013b, MNRAS, 428, 2795 [NASA ADS] [CrossRef] [Google Scholar]
 Kara, E., Fabian, A. C., Cackett, E. M., et al. 2013c, MNRAS, 434, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Kara, E., Zoghbi, A., Marinucci, A., et al. 2015, MNRAS, 446, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Marinucci, A., Matt, G., Kara, E., et al. 2014, MNRAS, 440, 2347 [NASA ADS] [CrossRef] [Google Scholar]
 McHardy, I. M., Papadakis, I. E., Uttley, P., Page, M. J., & Mason, K. O. 2004, MNRAS, 348, 783 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, L., Turner, T. J., Reeves, J. N., & Braito, V. 2010a, MNRAS, 408, 1928 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, L., Turner, T. J., Reeves, J. N., et al. 2010b, MNRAS, 403, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, S., & Kitamoto, S. 1989, Nature, 342, 773 [NASA ADS] [CrossRef] [Google Scholar]
 Nowak, M. A., & Vaughan, B. A. 1996, MNRAS, 280, 227 [NASA ADS] [Google Scholar]
 Nowak, M. A., Vaughan, B. A., Wilms, J., Dove, J. B., & Begelman, M. C. 1999, ApJ, 510, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Papadakis, I. E., Nandra, K., & Kazanas, D. 2001, ApJ, 554, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Priestley, M. B. 1981, Spectral Analysis and Time Series (London: Academic Press) [Google Scholar]
 Sriram, K., Agrawal, V. K., & Rao, A. R. 2009, ApJ, 700, 1042 [NASA ADS] [CrossRef] [Google Scholar]
 Timmer, J., & Koenig, M. 1995, A&A, 300, 707 [NASA ADS] [Google Scholar]
 Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014, A&A Rev., 22, 72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaughan, B. A., & Nowak, M. A. 1997, ApJ, 474, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Fabian, A. C., Uttley, P., et al. 2010, MNRAS, 401, 2419 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Uttley, P., & Fabian, A. C. 2011, MNRAS, 412, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Fabian, A. C., Reynolds, C. S., & Cackett, E. M. 2012, MNRAS, 422, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Reynolds, C., & Cackett, E. M. 2013a, ApJ, 777, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Reynolds, C., Cackett, E. M., et al. 2013b, ApJ, 767, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Zoghbi, A., Cackett, E. M., Reynolds, C., et al. 2014, ApJ, 789, 56 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: The crossspectrum 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 straightforward 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 crossperiodogram
Using Eqs. (9) and (10), we may write the crossperiodogram 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 recast 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 socalled Fejér kernel. One of the properties of this function is that, as N → ∞, F_{N}(ν) → δ(ν) (the Dirac δfunction). Thus, as N → ∞ the righthand 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 crossperiodogram 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 wellknown fact that the addition of socalled 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_{ϵX}h_{ϵX}). In the case of powerlaw 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 crossperiodogram (although they will bias the phaselag 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 righthand 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
Fig. D.1 Relative timelag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative timelag bias for the LS102.4 (black curves), LS3.2 (red curves), and OB (brown curves) light curves. Bottom row: relative timelag bias for the LS40.8 (black curves), LS20.4 (red curves), LS10.2 (brown curves), and LS5.1 (green curves) light curves. The horizontal dotted lines in this, and all similar subsequent figures, indicate the 0.1 relative timelag bias. 

Open with DEXTER 
Fig. D.2 As in Fig. D.1, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). 

Open with DEXTER 
Fig. D.3 As in Fig. D.1, for experiments PLD4 (first column) and PLD5 (second column). 

Open with DEXTER 
Fig. D.4 As in Fig. D.1, for experiments THRF1 (first column) and THRF2 (second column). 

Open with DEXTER 
Fig. D.5 Relative timelag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative timelag bias for the LS102.4 light curves with a smoothing parameter m = 20 (open black circles), and OB light curves with a smoothing parameter m = 13 (filled brown squares). Bottom row: relative timelag bias for the LS40.8 (black curves) and LS20.4 (red curves) light curves in the case of m = 20 segments. 

Open with DEXTER 
Fig. D.6 As in Fig. D.5, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). 

Open with DEXTER 
Fig. D.7 As in Fig. D.5, for experiments PLD4 (first column) and PLD5 (second column). 

Open with DEXTER 
Fig. D.8 As in Fig. D.5, for experiments THRF1 (first column) and THRF1 (second column). 

Open with DEXTER 
Fig. D.9 Mean sample coherence for experiments CD2 (first column), CD3 (second column), and PLD5 (third column). Top row: mean sample coherence for the LS102.4 light curves with a smoothing parameter m = 20 (open black circles), and OB light curves with a smoothing parameter m = 13 (filled brown squares). Bottom row: mean sample coherence for the LS40.8 (black curves) and LS20.4 (red curves) light curves estimated from m = 20 segments. 

Open with DEXTER 
Fig. D.10 Top row: relative timelag bias for the LS20.4 (continuous black curves) and LS40.8 (dashed red curves) light curves in experiment CD1. Bottom row: mean sample coherence for the LS20.4 (continuous black curves) and LS40.8 (dashed red curves) light curves in experiment CD1. Different columns correspond to different { (S/N)_{x},(S/N)_{y} } combinations, while the estimates were determined from m = 20 segments in each case. The vertical dashed lines in all panels indicate ν_{crit}, while the blue dotteddashed lines in the lower panels indicate the 1.2/(1 + 0.2m) mean sample coherence value (see Sects. 8 and 9 for details). 

Open with DEXTER 
Fig. D.11 As in Fig. D.10, for experiment PLD2. 

Open with DEXTER 
Fig. D.12 As in Fig. D.10, for experiment THRF1. 

Open with DEXTER 
Fig. D.13 The mean sample coherence, timelag error ratio and probability that the value of the timelag estimate is within 1σ of the sample mean (top, middle, and bottom row, respectively) for different values of m and light curve types. In each case, (S/N)_{x} = (S/N)_{y} = 3. The vertical dashed lines indicate ν_{crit} (see Sect. 8 for details). 

Open with DEXTER 
Fig. D.14 As in Fig. D.13, for (S/N)_{x} = 9 and (S/N)_{y} = 3. 

Open with DEXTER 
Fig. D.15 As in Fig. D.13, for (S/N)_{x} = 9 and (S/N)_{y} = 9. 

Open with DEXTER 
Fig. D.16 As in Fig. D.13, for (S/N)_{x} = 18 and (S/N)_{y} = 3. 

Open with DEXTER 
Fig. D.17 As in Fig. D.13, for (S/N)_{x} = 18 and (S/N)_{y} = 9. 

Open with DEXTER 
Fig. D.18 The excess kurtosis, skewness and probability that the averaged timelag estimate is Gaussiandistributed according to the KS test (top, middle, and bottom row, respectively) for different values of m and light curve types. In each case, (S/N)_{x} = (S/N)_{y} = 3. The vertical dashed lines indicate ν_{crit} (see Sect. 9 for details). 

Open with DEXTER 
Fig. D.19 As in Fig. D.18, for (S/N)_{x} = 9 and (S/N)_{y} = 3. 

Open with DEXTER 
Fig. D.20 As in Fig. D.18, for (S/N)_{x} = 9 and (S/N)_{y} = 9. 

Open with DEXTER 
Fig. D.21 As in Fig. D.18, for (S/N)_{x} = 18 and (S/N)_{y} = 3. 

Open with DEXTER 
Fig. D.22 As in Fig. D.18, for (S/N)_{x} = 18 and (S/N)_{y} = 9. 

Open with DEXTER 
All Figures
Fig. 1 Sample mean of the real and imaginary parts of the crossperiodogram (top left and right panels, respectively), the sample mean timelag spectrum and the relative timelag bias (bottom left and right panel, respectively) in experiment PLD2. The vertical black and red dashed lines indicate ν_{Nyq}/ 2 and ν_{Nyq}/ 5, respectively. Above these frequencies, the LS102.4 (black curve) and LS102.42 (red curve) relative timelag bias begins to noticeably increase (see text for more details). The horizontal dotted line in the bottom right panel, and in all subsequent plots, indicates the 0.1 (i.e. 10%) relative timelag bias. 

Open with DEXTER  
In the text 
Fig. 2 Mean relative timelag bias over all frequencies below ν_{max}, plotted as a function of ν_{max}, for different light curve types in various numerical experiments. The left and right dashed vertical lines indicate ν_{Nyq}/ 2 for the OB and LStype light curves, respectively. 

Open with DEXTER  
In the text 
Fig. 3 Mean relative timelag bias, plotted as function of light curve duration for LStype and OB light curves. 

Open with DEXTER  
In the text 
Fig. 4 Relative timelag bias for the LS and OBtype light curves (top and bottom panel, respectively), for various durations in experiment CD1. Filled squares indicate the mean relative timelag bias over the full sampled frequency range, evaluated at the mean logarithmic frequency. The dashed vertical line in the top and bottom panel indicate ν_{Nyq}/ 2 for the LS and OBtype light curves, respectively. 

Open with DEXTER  
In the text 
Fig. 5 Top left panel: mean sample timelag spectrum for the LS102.4 light curves in experiment CD3 (continuous black line). The black dashed line indicates the model timelag spectrum. Top right and bottom panels: the probability distribution of the phaselag estimates at the frequencies indicated by the vertical dashed lines in the topleft panel. The solid and dashed vertical line in these panels indicate the model and sample mean phaselag value at these frequencies, respectively. 

Open with DEXTER  
In the text 
Fig. 6 As in Fig. 1, for experiment PLD5. The blue dashed lines indicate the model CS (upper panels) and timelag spectrum (lower left panel). The solid brown line in the bottom left panel indicates the model timelag spectrum without taking the effects of phaseflipping into account. 

Open with DEXTER  
In the text 
Fig. 7 Mean relative bias of the m = 20 smoothed (top panel) and averaged (bottom panel) timelag estimates, plotted as a function of the relative bias of the nonsmoothed/averaged timelag estimates (the dashed lines show the onetoone relation line) in various numerical experiments. 

Open with DEXTER  
In the text 
Fig. 8 As in Figs. 1 and 6, for the m = 20 smoothed estimates obtained from the LS102.4 light curves in experiment CD1. The blue dashed lines indicate the model CS (upper panels) and timelag spectrum (lower left panel). 

Open with DEXTER  
In the text 
Fig. 9 Relative bias of the averaged timelag estimates for the LS20.4 and LS40.8 light curves, plotted as a function of (S/N)_{xy} in experiment THRF1. The horizontal dashed lines correspond to the value of the relative bias when (S/N)_{xy} → ∞. 

Open with DEXTER  
In the text 
Fig. 10 Probability distribution of the averaged phaselag estimate for the LS40.8 light curves in experiment THRF1 at different frequencies. 

Open with DEXTER  
In the text 
Fig. 11 Mean sample coherence (m = 20) obtained from the LS40.8 light curves in experiment THRF1 for different S/N combinations. The blue horizontal dashed lines indicate the value 1 /m (see Sect. 7.2 for details). 

Open with DEXTER  
In the text 
Fig. 12 Relative timelag bias for the LS102.4 light curves and different values of the model PSD parameters, { A,ν_{b} }, in experiment CD1. 

Open with DEXTER  
In the text 
Fig. D.1 Relative timelag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative timelag bias for the LS102.4 (black curves), LS3.2 (red curves), and OB (brown curves) light curves. Bottom row: relative timelag bias for the LS40.8 (black curves), LS20.4 (red curves), LS10.2 (brown curves), and LS5.1 (green curves) light curves. The horizontal dotted lines in this, and all similar subsequent figures, indicate the 0.1 relative timelag bias. 

Open with DEXTER  
In the text 
Fig. D.2 As in Fig. D.1, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). 

Open with DEXTER  
In the text 
Fig. D.3 As in Fig. D.1, for experiments PLD4 (first column) and PLD5 (second column). 

Open with DEXTER  
In the text 
Fig. D.4 As in Fig. D.1, for experiments THRF1 (first column) and THRF2 (second column). 

Open with DEXTER  
In the text 
Fig. D.5 Relative timelag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative timelag bias for the LS102.4 light curves with a smoothing parameter m = 20 (open black circles), and OB light curves with a smoothing parameter m = 13 (filled brown squares). Bottom row: relative timelag bias for the LS40.8 (black curves) and LS20.4 (red curves) light curves in the case of m = 20 segments. 

Open with DEXTER  
In the text 
Fig. D.6 As in Fig. D.5, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). 

Open with DEXTER  
In the text 
Fig. D.7 As in Fig. D.5, for experiments PLD4 (first column) and PLD5 (second column). 

Open with DEXTER  
In the text 
Fig. D.8 As in Fig. D.5, for experiments THRF1 (first column) and THRF1 (second column). 

Open with DEXTER  
In the text 
Fig. D.9 Mean sample coherence for experiments CD2 (first column), CD3 (second column), and PLD5 (third column). Top row: mean sample coherence for the LS102.4 light curves with a smoothing parameter m = 20 (open black circles), and OB light curves with a smoothing parameter m = 13 (filled brown squares). Bottom row: mean sample coherence for the LS40.8 (black curves) and LS20.4 (red curves) light curves estimated from m = 20 segments. 

Open with DEXTER  
In the text 
Fig. D.10 Top row: relative timelag bias for the LS20.4 (continuous black curves) and LS40.8 (dashed red curves) light curves in experiment CD1. Bottom row: mean sample coherence for the LS20.4 (continuous black curves) and LS40.8 (dashed red curves) light curves in experiment CD1. Different columns correspond to different { (S/N)_{x},(S/N)_{y} } combinations, while the estimates were determined from m = 20 segments in each case. The vertical dashed lines in all panels indicate ν_{crit}, while the blue dotteddashed lines in the lower panels indicate the 1.2/(1 + 0.2m) mean sample coherence value (see Sects. 8 and 9 for details). 

Open with DEXTER  
In the text 
Fig. D.11 As in Fig. D.10, for experiment PLD2. 

Open with DEXTER  
In the text 
Fig. D.12 As in Fig. D.10, for experiment THRF1. 

Open with DEXTER  
In the text 
Fig. D.13 The mean sample coherence, timelag error ratio and probability that the value of the timelag estimate is within 1σ of the sample mean (top, middle, and bottom row, respectively) for different values of m and light curve types. In each case, (S/N)_{x} = (S/N)_{y} = 3. The vertical dashed lines indicate ν_{crit} (see Sect. 8 for details). 

Open with DEXTER  
In the text 
Fig. D.14 As in Fig. D.13, for (S/N)_{x} = 9 and (S/N)_{y} = 3. 

Open with DEXTER  
In the text 
Fig. D.15 As in Fig. D.13, for (S/N)_{x} = 9 and (S/N)_{y} = 9. 

Open with DEXTER  
In the text 
Fig. D.16 As in Fig. D.13, for (S/N)_{x} = 18 and (S/N)_{y} = 3. 

Open with DEXTER  
In the text 
Fig. D.17 As in Fig. D.13, for (S/N)_{x} = 18 and (S/N)_{y} = 9. 

Open with DEXTER  
In the text 
Fig. D.18 The excess kurtosis, skewness and probability that the averaged timelag estimate is Gaussiandistributed according to the KS test (top, middle, and bottom row, respectively) for different values of m and light curve types. In each case, (S/N)_{x} = (S/N)_{y} = 3. The vertical dashed lines indicate ν_{crit} (see Sect. 9 for details). 

Open with DEXTER  
In the text 
Fig. D.19 As in Fig. D.18, for (S/N)_{x} = 9 and (S/N)_{y} = 3. 

Open with DEXTER  
In the text 
Fig. D.20 As in Fig. D.18, for (S/N)_{x} = 9 and (S/N)_{y} = 9. 

Open with DEXTER  
In the text 
Fig. D.21 As in Fig. D.18, for (S/N)_{x} = 18 and (S/N)_{y} = 3. 

Open with DEXTER  
In the text 
Fig. D.22 As in Fig. D.18, for (S/N)_{x} = 18 and (S/N)_{y} = 9. 

Open with DEXTER  
In the text 