Issue |
A&A
Volume 591, July 2016
|
|
---|---|---|
Article Number | A113 | |
Number of page(s) | 29 | |
Section | Numerical methods and codes | |
DOI | https://doi.org/10.1051/0004-6361/201527665 | |
Published online | 23 June 2016 |
Statistical properties of Fourier-based time-lag estimates
1
Department of Physics and Institute of Theoretical and Computational
Physics, University of Crete, 71003
Heraklion, Greece
e-mail:
epitrop@physics.uoc.gr
2
IESL, Foundation for Research and Technology-Hellas,
71110
Heraklion, Crete, Greece
Received: 30 October 2015
Accepted: 2 April 2016
Context. The study of X-ray time-lag 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 super-massive black hole) in these objects. To obtain reliable information from these studies, the statistical properties of time-lags estimated from data must be known as accurately as possible.
Aims. We investigated the statistical properties of Fourier-based time-lag estimates (i.e. based on the cross-periodogram), using evenly sampled time series with no missing points. Our aim is to provide practical “guidelines” on estimating time-lags 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 signal-to-noise ratios and sampling patterns/durations similar to those offered by AGN observations with present and past X-ray satellites. We also considered a range of different model time-lag spectra commonly assumed in X-ray analyses of compact accreting systems.
Results. Discrete sampling, binning and finite light curve duration cause the mean of the time-lag estimates to have a smaller magnitude than their intrinsic values. Smoothing (i.e. binning over consecutive frequencies) of the cross-periodogram can add extra bias at low frequencies. The use of light curves with low signal-to-noise ratio reduces the intrinsic coherence, and can introduce a bias to the sample coherence, time-lag estimates, and their predicted error.
Conclusions. Our results have direct implications for X-ray time-lag studies in AGN, but can also be applied to similar studies in other research fields. We find that: a) time-lags 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 cross-periodogram should be avoided, as this may introduce significant bias to the time-lag estimates, which can be taken into account by assuming a model cross-spectrum (and not just a model time-lag spectrum); c) time-lags should be estimated by dividing observed time series into a number, say m, of shorter data segments and averaging the resulting cross-periodograms; d) if the data segments have a duration ≳ 20 ks, the time-lag bias is ≲15% of its intrinsic value for the model cross-spectra and power-spectra considered in this work. This bias should be estimated in practise (by considering possible intrinsic cross-spectra that may be applicable to the time-lag spectra at hand) to assess the reliability of any time-lag analysis; e) the effects of experimental noise can be minimised by only estimating time-lags 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, time-lags 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 time-lags as a function of temporal frequency between X-ray light curves in different energy bands has frequently been used in the last two decades to probe the emission mechanism and geometry of the innermost region in active galactic nuclei (AGN; e.g. Papadakis et al. 2001; McHardy et al. 2004; Arévalo et al. 2006, 2008; Sriram et al. 2009) and X-ray binaries (XRBs; e.g. Miyamoto & Kitamoto 1989; Nowak & Vaughan 1996; Nowak et al. 1999). In the last few years, time-lag studies have revealed that “soft” X-ray variations in AGN are delayed with respect to “hard” X-ray variations at frequencies higher than ≈ 10-4 Hz (e.g. Fabian et al. 2009; Zoghbi et al. 2010, 2011; Emmanoulopoulos et al. 2011; De Marco et al. 2013). These time-lags are commonly referred to as soft lags, and have been observed in ≈ 20 sources so far. They have attracted considerable attention, since they are thought to be direct evidence of X-ray reverberation. Further credence to such a scenario came with recent discoveries of time-lags between the Fe Kα emission line (≈ 5−7 keV) and the X-ray continuum (e.g. Zoghbi et al. 2012; Kara et al. 2013b,a,c; Zoghbi et al. 2013b; Marinucci et al. 2014), as well as between the so-called Compton hump (≈ 10−30 keV) and the X-ray continuum (e.g. Zoghbi et al. 2014; Kara et al. 2015) in a few AGN.
The study of X-ray reverberation time-lags can provide valuable geometrical and physical information on the X-ray source and reflector, since they should depend on, for example, their typical size, relative distance, proximity to the central black hole (BH), the mass and spin of the BH, and the inclination of the system. To obtain this information, the statistical properties of time-lags estimated from observed light curves must be known as accurately as possible. For example, one must know their bias (i.e. how accurately their mean approximates the intrinsic time-lag spectrum), error, and probability distribution. The later is necessary if one wishes to fit the observed time-lag spectra with theoretical models. To the best of our knowledge, such a detailed investigation has not yet been performed.
Results from preliminary studies along these lines have been presented by Alston et al. (2013) and Uttley et al. (2014). In this paper we report the results of a more detailed study regarding the statistical properties of Fourier-based time-lag estimates, based on both analytical work and extensive numerical simulations. Our primary goals are: a) to investigate whether the frequently used time-lag estimates are indeed reliable estimates of the intrinsic time-lag spectrum; b) to study the effects of light curve sampling patterns and duration, as well as measurement errors, on the statistical properties of these estimates; and c) to provide practical guidelines which ensure estimates that are minimally biased, have know errors, and follow a Gaussian distribution. The latter property would be desirable, for example, in the case of model fitting using traditional χ2 minimisation techniques. Our work should have direct impact on current time-lag studies in the area of AGN X-ray timing analyses. We believe that our results can apply equally well to all scientific fields where similar techniques are employed to search for delays between two observed time-varying signals.
2. Definitions
Consider a continuous bivariate random time series { X(t),Y(t) }. We assume that the mean values (μX and μY), as well as the auto-covariance functions (RX(τ) and RY(τ); τ is the so-called lag) of the individual time series are finite and time-independent (i.e. they are stationary random processes). A random function that is frequently used to quantify the correlation between two time series in the time domain is the cross-covariance function (CCF), (1)where E is the expectation operator. We assume that the CCF depends on τ only, i.e. the CCF does not vary with time. Defined as above, RXY(τ) is the CCF with Y(t) leading X(t). The Fourier transform of the CCF,
(2)defines the cross-spectrum (CS) of the bivariate processes. The CCF is not necessarily symmetric about τ = 0, and hence the CS is generally a complex number that can be written as
(3)The real functions cXY(ν) and {−qXY(ν) } represent the real, ℜ [ hXY(ν) ], and imaginary, ℑ [ hXY(ν) ], parts of the CS, respectively. The function cXY(ν) is an even function of ν, while qXY(ν) is an odd function of ν (Priestley 1981, P81 hereafter). The quantity | hXY(ν) | is the CS amplitude, and φXY(ν) is the phase-lag spectrum, which is defined as
(4)The phase-lag spectrum represents the average phase shift between sinusoidal components of the two time series with frequency ν. Since it is defined modulo 2π, it is customary to define its principal value in the interval (−π,π ], which is the convention adopted in the present work as well. The time-lag spectrum is defined as
(5)and represents the average temporal delay between sinusoidal components of the two time series with frequency ν. Given the definition of the CCF by Eq. (1), a positive time-lag value at a frequency ν indicates that, on average, the sinusoidal component of X(t) lags behind the respective component of Y(t).
Another statistic that is often used to study correlations between two random processes in Fourier space is the so-called coherence function, (6)where hX(ν) and hY(ν) denote the power-spectral density functions (PSDs) of the time series X(t) and Y(t), respectively (the PSD is defined as the Fourier transform of the auto-covariance function of a random process). The coherence function is interpreted as the correlation coefficient between two processes at frequency ν, as it measures the degree of linear correlation between them at each frequency. As with ordinary correlation coefficients,
, where
indicates a perfect (linear) correlation, and
implies the absence of any correlation at frequency ν.
2.1. The effects of binning and discrete sampling
In practice, the data correspond to values of a single realisation of a random process that are recorded over a finite time interval, T (the duration). The recording is typically performed at regular time intervals, Δtsam (the sampling period).
Consider a pair of observed time series { x(tr),y(tr) }, where tr = rΔtsam, 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(tr),y(tr) } correspond to a particular, finite realisation of a discrete version of the intrinsic process, which we denote by { X(ts),Y(ts) }, where ts = sΔtsam and s = 0, ± 1, ± 2,.... In addition, light curves1 in astronomy are the result of binning the intrinsic signal over time bins of size Δtbin2. In this case, the relation between X(ts) and X(t) is given by (7)with an identical relation holding between Y(ts) and Y(t).
It is rather straight-forward to show that { X(ts),Y(ts) } 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, hXY(ν), which is defined only at frequencies | ν | ≤ 1/2Δtsam, is given by (8)In other words, at each frequency ν within the aforementioned interval, hXY(ν) is the superposition of the intrinsic CS values, hXY(ν), at frequencies ν,ν ± (1/Δtsam),ν ± (2/Δtsam),.... This is entirely analogous to the so-called aliasing effect in the case of the PSD. However, although PSDs are always positive, hence aliasing in this case always implies transfer of “power” from higher to the sampled frequencies, the effects of aliasing on hXY(ν) cannot be predicted without a priori knowledge of hXY(ν). Aliasing affects both cXY(ν) and qXY(ν), which are not necessarily positive at all frequencies. The situation is further complicated by the fact that while cXY(ν) is an even function of ν, qXY(ν) 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 hXY(ν) that are aliased in the sampled frequency range are suppressed by the sinc2 function, whose argument depends on the bin size, Δtbin. Since the sinc2 function asymptotically approaches unity in the limit whereby its argument goes to zero, if the observed time series are not binned (i.e. when Δtbin → 0) the sinc2 term in the right-hand side of Eq. (8) vanishes, and hence the effects of aliasing on hXY(ν) are maximised.
3. Estimation of cross-spectra; the cross-periodogram
The discrete Fourier transform (DFT) of a time series x(tr) is defined as (9)where
is the sample mean of x(tr). It is customary to estimate DFTs at the following set of frequencies: νp = p/NΔtsam, where p = 1,2,...,N/ 2, and νNyq = 1/2Δtsam is the so-called Nyquist frequency, which corresponds to the highest frequency that can be probed for a given sampling period, Δtsam (the symbol νp will henceforth always stand for this particular set of frequencies). The cross-periodogram of two time series is defined as
(10)where the asterisk ∗ denotes complex conjugation. The cross-periodogram is used in practice as an estimator of the CS. Given Eqs. (4) and (5), it seems reasonable to use the real and imaginary parts of Ixy(νp) and accept
as estimators of the phase- and time-lag spectrum, respectively.
3.1. The bias of cross-spectral estimators
As we show in Appendix B, (13)where hXY(ν) is the CS of { X(ts),Y(ts) } (as defined by Eq. (8)). The function FN(ν′ − νp), which is called the Fejér kernel, has a large peak at ν′ = νp of magnitude NΔtsam, and decays to zero as | ν′ | → ∞. Subsequent peaks appear at ν′ ≈ νp ± 3/2Δtsam,νp ± 5/2Δtsam,..., while the zeros of the function occur at ν′ = νp ± 1 /NΔtsam,νp ± 2 /NΔtsam,.... As N → ∞, FN(ν′ − νp) → δ(ν′ − νp) (the Dirac δ-function). Therefore, in the limit N → ∞ the right-hand side of Eq. (13) converges to hXY(νp). The cross-periodogram is therefore an asymptotically (when N → ∞) unbiased estimate of the modified (owing to the effects of binning and discrete sampling) intrinsic CS.
Since hXY(ν) is not equal to hXY(ν), 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Δtsam of a given pair of observed time series will be sufficiently long for the effects of the convolution of hXY(ν) with the Fejér kernel (Eq. (13)) to be significant of not, unless hXY(ν) is known a priori. We define the “bias” of the cross-periodogram as bI(νp) ≡ E [ Ixy(νp) ] − hXY(νp) (henceforth, the term bias for a statistical estimator will always refer to the difference between its mean and intrinsic value). Even if bI(νp) was known, and hence could be used to “correct” Ixy(νp), it would still not be straightforward to determine the bias of or
, since E { arg [ Ixy(νp) ] } is not necessarily equal to arg { E [ Ixy(νp) ] }.
To quantify the bias of time-lag estimates based on the cross-periodogram, we performed an extensive number of simulations, which we describe below. The characteristics of the simulated time series (i.e. time bin size, sampling rate, duration, intrinsic CS, and PSDs) are representative of those observed in AGN X-ray light curves. However, most of our results should apply to any time series observed in practice (see discussion in Sect. 10).
4. Simulating correlated random processes
The simulations we performed are based on the procedure outlined by Timmer & Koenig (1995)3 to generate artificial realisations of a discrete stationary random process with a specified model PSD, hX(ν), number of points, N, sampling period, Δtsam, and mean count-rate, μX. We assumed a model PSD of the form (14)This function describes a power law with a low-frequency slope of −1 which smoothly “bends” to a slope of −2 at frequencies above the bend-frequency, νb, as in the case of inferred X-ray PSDs of most AGN. We chose an amplitude value of A = 0.01, and assumed that νb = 2 × 10-4 Hz. Furthermore, we set N = 10.24 × 106 and Δtsam = 1 s. By construction, the intrinsic PSD of the generated light curves is discrete, and has non-zero values only at frequencies νj = j/NΔtsam, where j = ± 1, ± 2,..., ± N/ 2.
For each simulated light curve we also created a corresponding “partner” with a mean count-rate μY, assuming a specified model phase-lag spectrum, φXY(ν). This was achieved by multiplying the DFT of the first light curve by the factor (μY/μX)e− iφXY(ν). The light curves thus have the same PSD shape, and can be considered as realisations of a discrete process whose intrinsic CS is equal to hXY(ν) = (μY/μX)hX(ν)eiφXY(ν). Since the CS, in effect, represents the average product of the Fourier component amplitudes of each light curve, and this amplitude is proportional to the PSD, the CS is, by construction, non-zero only at the same frequencies where the PSD is non-zero as well. The light curve pairs generated following this procedure will henceforth be referred to as original realisations.
4.1. The model time-lag spectra
We considered three different model time-lag spectra:
-
1.
Constant delays: in this case we assumed that φXY(ν) = 2πνd, where d = 10 s, 150 s and 550 s is aconstant delay between the light curve pairs (henceforth,experiments CD1, CD2, and CD3, respectively).
-
2.
Power law delays: in this case we assumed that φXY(ν) = (2πν)Bν− β, where B is the normalisation and β the power-law index of the corresponding model time-lag spectrum. We considered the values { B,β } = { 0.001,1 }, { 0.01,1 }, { 0.1,1 }, { 0.01,0.5 }, { 0.01,1.5 } (henceforth, experiments PLD1, PLD2, PLD3, PLD4, and PLD5, respectively).
-
3.
Top-hat response functions: in this case we assumed that φXY(ν) = arg [ 1 + fei2πνt0sinc(πνΔ) ]. This phase-lag spectrum is expected if the intrinsic time series are related by the following equation:
, where Ψ(t) (the so-called response function) is a simple so-called top-hat function, i.e. it has a constant value of f/ Δ in the interval | t − t0 | ≤ Δ/2 and zero otherwise (Δ is the width of the top-hat and t0 its centroid). We considered the model parameter values { f,t0,Δ } = { 0.2,200 s,200 s } and { 0.2,2000 s,2000 s } (henceforth, experiments THRF1 and THRF2, respectively).
4.2. The sampling patterns of the simulated light curves
For each model time-lag spectrum listed above, we created 30 light curve pairs as per the specifications given in Sect. 4. To simulate light curves encountered in practice, the original realisations need to be properly “chopped”, binned, and sampled.
Most of the data in X-ray astronomy are provided by satellites in low-Earth orbit with a typical orbital period of ≈ 96 min and bin size of 16 s, such as ASCA, RXTE, Suzaku and NuSTAR. Light curves obtained from observations with these satellites are “affected” by periodic Earth occultations of a target during every orbit. As a result, they contain gaps that are typically ≈ 1−3 ks long. They are hence not appropriate for Fourier analysis using the “standard” techniques we considered in this work. One can bin the data at one orbital period to acquire evenly sampled light curves with no missing points. Such light curves can be used to probe variability on long time-scales. On the other hand, the data usually contain an appreciable number of continuously sampled segments, with a duration typically ≲3 ks. These segments can, in principle, be used to probe variability on short timescales.
Contrary to low-Earth orbit satellites, XMM-Newton observations result in continuously sampled light curves up to ≈ 120 ks long owing to its highly elliptical orbit. We can use such a pair of light curves to compute the cross-periodogram and time-lag estimates directly, or “chop” them into segments of shorter duration, calculate the cross-periodogram for each segment and then bin the resulting estimates at certain frequencies to estimate the time-lags (see Sect. 6 for a more detailed discussion on this issue).
The purpose of the simulations we performed is to study the sampling properties of time-lags when estimated using light curves with durations and sampling patterns similar to those described above. To this end, we adopted the following strategy:
-
1.
Chop each of the original realisations into100 parts and bin them at 100 s. This process generates3000 light curve pairs with T = 102.4 ks and Δtsam = Δtbin = 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 Δtsam = Δtbin = 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 Δtsam = Δtbin = 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 Δtsam = Δtbin = 5760 s (OB light curves, hereafter).
The reason for the original realisations being longer and with a finer sampling rate than the “final” light curves is to simulate the effects of binning and finite light curve duration on the bias of the time-lag estimates, (henceforth, the time-lag bias).
5. The bias of the time-lag estimates in practice
To quantify , we calculated
using Eq. (12) for the 10 numerical experiments and each light curve type described in Sects. 4.1 and 4.2. We then computed the sample mean at each frequency,
(quantities in angle brackets will hereafter denote their sample mean). We chose to study the properties of
mainly, since the majority of works concerned with AGN X-ray timing studies use this estimator. In addition, we computed the quantity
to quantify the time-lag bias in terms of its intrinsic value (henceforth, the relative bias). This allows us to directly compare the time-lag bias obtained from different light curve types in the various numerical experiments we considered.
Figures D.1–D.4 show for the LS102.4/LS3.2/OB (top rows; continuous black, red, and brown curves, respectively) and LS40.8/20.4/10.2/5.1 (bottom rows; continuous black, red, brown, and green curves, respectively) light curves. Each column in these figures corresponds to a different numerical experiment. The main results of the simulations we performed are summarised below.
5.1. Light curve binning and aliasing
A common feature in all mean sample time-lag spectra is that they decrease to zero at high frequencies, irrespective of the light curve sampling pattern, duration, and model time-lag spectrum. This decrease is due to the effects of aliasing and light-curve binning.
To demonstrate this issue, we considered experiment PLD2 (the results are similar for the other numerical experiments as well). We generated two additional ensembles of light curve pairs that have the same length as the LS102.4 light curves. The first ensemble was constructed using the original realisations of experiment PLD2, which were chopped into 3000 parts of length 102.4 ks each, and sampled (not binned) every 100 s (henceforth, the LS102.4-2 light curves). These light curves are still affected by aliasing, since the original realisations have a finer sampling rate of 1 s. For the second ensemble, we constructed 30 new original realisations with the same length as the rest of the numerical experiments (i.e. 10.24 Ms) and a sampling period of 100 s. We subsequently chopped them into 3000 parts of length 102.4 ks each (henceforth, the LS102.4-3 light curves). Since the original realisations in this case have the same sampling rate as the LS102.4-3 light curves, they should not be affected by either binning or aliasing.
![]() |
Fig. 1 Sample mean of the real and imaginary parts of the cross-periodogram (top left and right panels, respectively), the sample mean time-lag spectrum and the relative time-lag 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.4-2 (red curve) relative time-lag bias begins to noticeably increase (see text for more details). The horizontal dotted line in the bottom right panel, and in all subsequent |
We then calculated and
for these two additional light curve types. Figure 1 shows ⟨ ℜ [ Ixy(νp) ] ⟩, ⟨ ℑ [ Ixy(νp) ] ⟩,
, and
(top left and right, bottom left and right panels, respectively) for the LS102.4, LS102.4-2, and LS102.4-3 light curves (black, red, and brown curves, respectively)4. The sample mean time-lag spectra for the LS102.4-2 and LS102.4 light curves decrease to zero at frequencies higher than ≈ 10-3 Hz (= νNyq/ 5) and ≈ 2.5 × 10-3 Hz (= νNyq/ 2), respectively. These frequencies are indicated by the vertical dashed lines in the same figure. This is not the case for the LS102.4-3 light curves, as shown by the brown curves in the bottom panels of Fig. 1. The time-lag bias of the sampled and binned light curves can be understood by the plots in the top panels of Fig. 1. On average, the gradients of the imaginary and real parts of the cross-periodogram enter steeper rates of decrease and increase, respectively, compared to their intrinsic values, around the frequencies indicated by the vertical dashed lines in Fig. 1. As a result, their ratio (which determines the phase-lag estimate) is decreased on average, hence the time-lag bias increases. This increase is most severe at high frequencies, where the effects of aliasing and light curve binning on the cross-periodogram are maximised. The time-lag bias owing to these effects is more pronounced for the LS102.4-2 light curves, since light curve binning suppresses the aliasing effect (see Sect. 3).
![]() |
Fig. 2 Mean relative time-lag 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 LS-type light curves, respectively. |
Figures D.1–D.4 show that the onset of the high-frequency increase in the relative time-lag bias occurs at ≈ νNyq/ 2 for all light curve types, model time-lag spectra, and light curve durations. To further illustrate this effect, we define the function as the mean relative time-lag bias at all frequencies below νmax. Figure 2 shows
, evaluated at νmax = νNyq/ 3, νNyq/ 2.5, νNyq/ 2, and νNyq/ 1.5, for the OB (diamonds), LS40.8 (squares), and LS102.4 (circles) light curves in experiments CD1 (black), PLD2 (red), and THRF1 (blue). The vertical red and black dashed lines indicate νNyq/ 2 for the OB and LS102.4/40.8 light curves, respectively. At frequencies below νNyq/ 2, the relative time-lag bias remains approximately constant, but increases markedly at higher frequencies (the same trend is observed in all numerical experiments we considered).
We conclude that time-lags should be estimated at frequencies ≲νNyq/ 2 to minimise the effects of light-curve binning and aliasing on the time-lag bias.
5.2. Finite light curve duration
The time-lag estimates are biased even at frequencies ≲νNyq/ 2. Figures D.1–D.4 show that the relative time-lag bias is generally ≲15% at frequencies ≲νNyq/ 2 for the LS102.4/40.8/20.4 and OB light curves, and larger for the rest. It is usually positive, in the sense that the estimates are, on average, smaller in absolute magnitude than their corresponding intrinsic values. The bias at intermediate/low frequencies is not due to light curve binning and/or aliasing (we note, for example, that at frequencies ≲νNyq/ 5 in Fig. 1 is identical for all three light curve types, which have the same length but different time bin sizes). This bias is due to the finite duration, T, of the light curves, or, technically speaking, due to the convolution of hXY(νp) with the Fejér kernel (see Eq. (13)).
To quantify the dependence of this bias on T, Fig. 3 shows as a function of T for the LS5.1/10.2/20.4/40.8/102.4 and OB light curves (filled and open points, respectively) in all experiments that do not exhibit so-called phase-flipping (see Sect. 5.3 for details). The dashed red line indicates the best-fit relation to the LS-type data:
. This relation fits the LS-type data quite well. We therefore conclude that the mean relative time-lag bias decreases with increasing light curve duration as
. However, this relation is not consistent with the points corresponding to the OB light curves. This is unexpected at first, since their duration (≈ 300 ks) is larger than the duration of all LS-type light curves. This discrepancy arises because the frequency range probed by the OB light curves is significantly lower than by the LS-type light curves (owing to their longer duration and larger time bin size), as well as the fact that the relative time-lag bias increases at lower frequencies.
To illustrate this effect, we generated additional ensembles of light curve pairs that have the same time bin size (100 s) and a duration two, three, four, and five times longer than LS102.4 light curves in the case of experiment CD1 (the results are identical for the other numerical experiments as well). The solid lines in the top panel of Fig. 4 show the relative time-lag bias at all frequencies below νNyq/ 2 for all these light curves (their duration is listed next to each curve). As expected, the relative bias at each frequency decreases with increasing T. For a fixed T, the relative bias increases with decreasing frequency. The filled boxes in the same plot indicate the mean relative bias over the whole frequency range, plotted at the mean logarithmic frequency. The mean relative bias, over the full frequency range, decreases with increasing T. In the bottom panel of the same figure we plot the same results in the case of OB-type light curves with a duration of ≈ 50−500 ks. The relative bias at each frequency is almost identical to the previous case, despite the larger time bin size of the OB-type light curves. Although the mean relative bias decreases with increasing T as before, owing to the larger time bin size in this case, we sample a lower frequency range, and hence the mean relative bias is larger than before (for the same T). We therefore conclude that the relative time-lag bias depends on , although the normalisation of this relation depends on the frequency range on which the time-lags are estimated; the lower the frequencies, the more biased the time-lag estimates will be. For the model CS we considered, our results suggest that, to obtain time-lag estimates with a relative bias ≲15%, the light curves must have T ≳ 20 ks.
![]() |
Fig. 3 Mean relative time-lag bias, plotted as function of light curve duration for LS-type and OB light curves. |
![]() |
Fig. 4 Relative time-lag bias for the LS- and OB-type light curves (top and bottom panel, respectively), for various durations in experiment CD1. Filled squares indicate the mean relative time-lag 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 OB-type light curves, respectively. |
5.3. Phase-flipping
Since the phase-lag spectrum is defined on the interval (−π,π ], there can be frequencies whereby this function will exceed the value of π and “flip back” to the value of −π (and vice-versa). This effect is known as phase-flipping (or phase-wrapping).
These kinds of events take place in experiments CD2, CD3, and PLD5 (see Figs. D.1 and D.3). The top left panel of Fig. 5 shows for the LS102.4 light curves in experiment CD3 (continuous black line), along with the corresponding model time-lag spectrum (black dashed line). The simulated light curves in this experiment are separated by a constant delay. As a result, the model phase-lag spectrum is a linearly increasing function of frequency, i.e. φXY(ν) = 2πνd. At ν = 1/2d ≈ 9.1 × 10-4 Hz we get φXY(ν) = π and then, by definition, φXY(ν) “jumps” to the value −π. Since τXY(ν) = φXY(ν)/2πν, the model time-lag spectrum will likewise jump from 550 s to −550 s. Subsequent phase-flips occur at frequencies ν = j/ 2d (j = 2,3,...), where the model time-lag spectrum undergoes discontinuous jumps from positive to negative values with decreasing amplitude.
The mean sample time-lag spectrum also fluctuates from positive to negative values at these frequencies, although the transition is much smoother. To illustrate the reason for this behaviour, in the top right and bottom panels of Fig. 5, we plot the probability distribution of the phase-lag estimate at the following frequencies: 8.1 × 10-4 Hz (top right panel), 9.1 × 10-4 Hz and 1.0 × 10-3 Hz (bottom left and right panel, respectively). These frequencies are indicated by the vertical dashed lines in the top left panel of the same figure. The middle frequency is very close to the frequency where the first phase-flip occurs, while the other two are lower and higher than that frequency. The solid and dashed vertical lines in the top-right and bottom panels of Fig. 5 show the model and sample mean phase-lag value at these frequencies, respectively. Since the phase-lag estimate is defined on the interval (−π,π ], the so-called wings of its distribution that exceed these boundaries (indicated by the “empty” histograms in the plots) are “folded back” into the allowed range. This causes the mean of the distribution to shift towards a value lower than that of the model phase-lag spectrum. This effect is most severe in the vicinity of frequencies where phase-flipping occurs, since, in this case, the mean of the phase-lag estimate is close to zero.
![]() |
Fig. 5 Top left panel: mean sample time-lag spectrum for the LS102.4 light curves in experiment CD3 (continuous black line). The black dashed line indicates the model time-lag spectrum. Top right and bottom panels: the probability distribution of the phase-lag estimates at the frequencies indicated by the vertical dashed lines in the top-left panel. The solid and dashed vertical line in these panels indicate the model and sample mean phase-lag value at these frequencies, respectively. |
Phase-flipping does not necessarily occur at increasingly higher frequencies. If, as in the case of experiments PLD4 and PLD5 for example, the intrinsic phase-lag spectrum increases (in absolute magnitude) with decreasing frequency, phase-flipping takes place at increasingly lower frequencies. This is seen in the bottom left panel of Fig. 6, where we plot for the LS102.4 light curves in experiment PLD5 (open black circles). The model time-lag spectrum for this experiment is τXY(ν) = 0.01ν-1.5 (indicated by the continuous brown line). The corresponding phase-lag spectrum is φXY(ν) = (2πν)0.01ν-1.5, and phase-flipping occurs at (4 × 10-4/j2) Hz (j = 1,2,...). At frequencies close to 4 × 10-4 Hz, where the first phase-flip occurs, the mean sample time-lag spectrum exhibits a jump which is smoother than the abrupt jump of the model time-lag spectrum, for exactly the same reasons we discussed above. This is not very clear in the plot shown in the bottom left panel of Fig. 6, but is evident in the plot of
(bottom-right panel in the same figure).
The dashed blue line in the top panels of Fig. 6 indicate the model real and imaginary parts of the CS. Interestingly, the sample mean of the real and imaginary parts of the cross-periodogram (indicated by the open circles in the top panels) are not biased at frequencies around ≈ 4 × 10-4 Hz. In fact, the argument of their ratio divided by the angular frequency, i.e. arg { E [ ĥXY(νp) ] } /2πνp, is very similar to the model time-lag spectrum, and exhibits a sharp jump at this frequency. This is a case where, for the reasons explained above, arg { E [ ĥXY(ν) ] } is not equal to E { arg [ ĥXY(ν) ] }.
![]() |
Fig. 6 As in Fig. 1, for experiment PLD5. The blue dashed lines indicate the model CS (upper panels) and time-lag spectrum (lower left panel). The solid brown line in the bottom left panel indicates the model time-lag spectrum without taking the effects of phase-flipping into account. |
The time-lag estimates at the two lowest frequencies are heavily biased owing to the bias of the real and imaginary parts of the cross-periodogram. This is reflected in the mean sample time-lag spectrum, which has a large relative bias at those frequencies (bottom right panel in Fig. 6). This bias originates from the convolution of the intrinsic CS with the Fejér kernel, which causes the mean real and imaginary parts of the cross-periodogram to diverge from the model CS because of its rapid oscillatory behaviour at low frequencies.
6. Smoothed/averaged cross-spectral estimators
The variance of the real and imaginary parts of the cross-periodogram is and
, respectively (P81). They are unknown, and do not depend on the number of points in the observed time series (i.e. they will not decrease if we increase their duration). In practice, we usually average a certain number, say m, of consecutive cross-periodogram estimates, i.e.
(15)and accept ĥxy(νk) as the CS estimator at the frequencies νk = (1 /m) ∑ pνp (k = 1,2,...,N/m). The process of binning over consecutive frequencies is called smoothing5. The real and imaginary parts of the smoothed CS estimator are correspondingly given by
An alternative procedure is to partition the available time series into m shorter segments of duration T/m, compute the cross-periodogram for each segment,
(l = 1,2,...,m, and νk = k/ (T/m), where k = 1,2,...,N/m), and then average the different cross-periodogram values at each νk:
(18)with the real and imaginary parts estimated as in the case of the smoothed estimates (to differentiate between estimators that are smoothed or averaged over individual segments, we henceforth refer to the former as smoothed and the latter as averaged estimates).
One can show that the variance of the real and imaginary parts of the smoothed/averaged cross-periodogram (as well as their covariance) is inversely proportional to m. As the duration of the observed light curves increases (i.e. N increases for a fixed Δtsam), m can be increased proportionally without degrading the frequency resolution. In this sense, the variance of the smoothed/averaged estimates decreases with increasing N (in the limit N → ∞, it approaches zero).
We can use the real and imaginary parts of the smoothed/averaged CS estimates to construct the following estimators of the phase- and time-lag spectrum: Their variance is given by the following asymptotic formulae (P81; Nowak et al. 1999; Bendat & Piersol 2011)
Equations (21) and (22) highlight the importance of
in constructing a reliable time-lag estimator. The lower the coherence, the higher the variance of the phase- and time-lag estimates will be (this is expected, as it must be more “difficult” to detect delays when the two light curves are highly incoherent). The “natural” choice for the coherence estimator, as suggested by Eq. (6), is the following:
(23)where ĥx(νk) and ĥy(νk) are the smoothed/averaged PSD estimators of the two observed time series. The variance of the coherence estimator is given by the asymptotic formula (P81; Vaughan & Nowak 1997; Nowak et al. 1999; Bendat & Piersol 2011)
(24)Since the intrinsic coherence,
, is unknown, it is customary to replace it in Eqs. (21) and (24) by its estimate,
, to obtain a numerical value.
6.1. Bias due to smoothing
To investigate whether smoothing or averaging increases the bias of the cross-spectral estimates, we smoothed the real and imaginary parts of the cross-periodograms estimated from the LS102.4 and OB light curves (for all numerical experiments) using Eqs. (16) and (17) for m = 20. We also used the cross-periodograms of the LS20.4/40.8 light curves, and averaged m = 20 of them to produce the corresponding averaged estimates. We then computed, in both cases and for all experiments, and
using Eqs. (20) and (23), as well as their corresponding sample mean values, and
. Figures D.5−D.8 show
for the smoothed (top panels; open white circles and filled brown squares for the LS102.4/OB light curves) and averaged estimates (bottom panels; continuous black and red curves for the LS20.4/40.8 light curves).
In general, smoothing increases the time-lag bias at low frequencies. The top panel in Fig. 7 shows vs.
for the smoothed estimates (in the case of the numerical experiments that do not exhibit phase-flipping in the sampled frequency range). The bottom panel shows the same plot for the averaged estimates. The blue dashed line shows the one-to-one relation in each case. Although the bias of the averaged estimates falls on the expected one-to-one relation line, this is not always the case for the smoothed estimates. Significant, additional bias appears in the case of experiments CD1, PLD4, and THRF1. The intrinsic CS in these cases exhibits a prominent non-linear variation with frequency. As a result, linear smoothing introduces a bias to the real and imaginary parts of the smoothed cross-periodogram. This is demonstrated in Fig. 8, which shows ⟨ ℜ [ ĥxy(νk) ] ⟩ (top left panel), ⟨ ℑ [ ĥxy(νk) ] ⟩ (top right panel),
(bottom left panel), and
(bottom right panel) for the m = 20 smoothed estimates of the LS102.4 light curves in experiment CD1. The intrinsic CS and time-lag spectrum are shown as blue dashed lines. At low frequencies (≲2 × 10-4 Hz) where the intrinsic CS has the highest curvature (i.e. where its second derivative has the largest value), the bias of the smoothed cross-periodogram and time-lag estimate is maximised.
![]() |
Fig. 7 Mean relative bias of the m = 20 smoothed (top panel) and averaged (bottom panel) time-lag estimates, plotted as a function of the relative bias of the non-smoothed/averaged time-lag estimates (the dashed lines show the one-to-one relation line) in various numerical experiments. |
The mean sample coherence at all frequencies is ≈ 1 in most, but not all, cases. Figure D.9 shows the mean sample coherence for experiments CD2, CD3, and PLD5 (filled circles and solid lines indicate the m = 20 smoothed and averaged estimates, respectively). Despite the fact that, by construction, at all frequencies, the mean sample coherence is significantly smaller than unity in these cases. This is difficult to explain and predict. The coherence bias is a complicated function of the
and
bias, as well as the bias of the smoothed/averaged PSD estimates. For example, although the mean real part of the smoothed cross-periodogram is biased in experiment CD1 when m = 20 (see the top left panel in Fig. 8), the mean sample coherence is not, presumably because it is counterbalanced by the bias of the smoothed PSD estimates.
7. The effect of measurement errors
Every measured signal will inevitably be affected by experimental “noise”. To study the effect of measurement errors on the time-lag and coherence estimates, we considered the LS20.4 and LS40.8 light curve pairs in experiments CD1, PLD2, and THRF1. For these experiments, we increased the number of original realisations from 30 to 100. This increase in number of simulated light curves for each numerical experiment results in a significant increase in the computing time. For that reason, in this and the following section we consider only experiments CD1, PLD2, and THRF1. They do not exhibit phase-flipping effects, and are representative of the three categories of model phase-lag spectra we consider. In total, we thus ended up with 5 × 104 LS20.4 and 2 × 104 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 time-lag spectrum (lower left panel). |
Furthermore, we created five copies of each of these pairs to simulate the effects of different signal-to-noise ratio (S/N) combinations, { (S/N)x,(S/N)y }, for each experiment. The effects of measurement errors were simulated by adding a Gaussian random number with zero mean to each point of every light curve, which were uncorrelated both with each other, as well as with the intrinsic light curve values. The variance of the random numbers was chosen such that there were five different S/N combinations for each light curve pair in the aforementioned experiments: { (S/N)x,(S/N)y } = { 3,3 }, { 9,3 }, { 18,3 }, { 9,9 }, and { 18,9 }.
We then calculated the m = 10, 20, 30, and 40 averaged cross-periodogram to estimate the time-lags and coherence according to Eqs. (20) and (23), along with their error estimates according to Eqs. (22) and (24). The number of averaged time-lag and coherence estimates were thus 5000 (2000), 2500 (1000), 1666 (666), and 1250 (500) for the LS20.4 (LS40.8) light curves in each experiment and every S/N combination. Figures D.10−D.12 show (top row) and
(bottom row) for the LS20.4/40.8 light curves (solid and dashed lines, respectively) for all S/N combinations, in each of the three experiments and for m = 20.
7.1. The effect of measurement errors on the time-lag bias
![]() |
Fig. 9 Relative bias of the averaged time-lag 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 → ∞. |
A comparison between, for example, the top panels in Fig. D.10 and the bottom left panel in Fig. D.5 reveals that the time-lag estimates become increasingly biased at high frequencies (in the sense that they converge to zero) with decreasing S/N. As the S/N of the light curves increases, becomes consistent with its corresponding value when there is no noise present. To illustrate the effect of noise on the time-lag bias, we computed
as a function of
for experiment THRF1 (this quantity is a measure of the combined S/N of both light curves). The results are shown in Fig. 9 for the LS20.4 and LS40.8 light curves. The horizontal dashed black and dotted-dashed red line correspond to the
value when there is no noise (i.e. when (S/N)xy → ∞) in the LS20.4 (filled black circles) and LS40.8 (open red squares) light curves, respectively. Clearly, the time-lag estimates are more biased when (S/N)xy ≲ 3. Identical trends are observed in the case of experiments CD1 and PLD2 as well.
This is a rather unexpected result, since the noise introduced to one light curve is independent of the noise introduced to the other. We would therefore not expect the time-lag bias to be affected by noise (although we would expect the resulting estimates to have a larger scatter around the mean, which is indeed the case; compare, for example, the scatter of the time-lags in the top panels in Fig. D.10 and in the bottom left panel of Fig. D.5). As in the case of phase-flipping (see Sect. 5.3), this bias arises because the phase-lag estimates are defined on the interval (−π,π ]. As an example, in Fig. 10 we show the probability distribution of phase-lag estimates for the LS40.8 light curves with (S/N)x = (S/N)y = 9 in experiment THRF1 at three different frequencies: 6.1 × 10-4 Hz, 1.2 × 10-3 Hz, and 2.5 × 10-3 Hz (top, middle, and bottom panel, respectively). Filled black bars and open red bars indicate the distribution of the m = 10 and m = 40 averaged phase-lag estimates, respectively. As the frequency increases, the mean sample coherence decreases (see the next section), and hence the scatter of the phase-lag estimate increases according to Eq. (21). When the magnitude of this scatter becomes sufficiently large, the “wings” of the distribution exceed the range of allowed values for the phase-lag estimate, and are thus folded back into the allowed range. Consequently, the distribution becomes increasingly uniform over the interval (−π,π ], and hence its mean converges to zero. This bias reduces as m increases, since, in this case, the scatter of the phase-lag estimates themselves, and hence their bias, decreases (see Eq. (21)). This is also evident from Fig. 10, as the “width” of the distributions is smaller in the case of larger m.
The vertical dashed lines in all panels of Figs. D.10−D.12 indicate the frequency at which the sample coherence becomes equal to 1.2/(1 + 0.2m) (this value is indicated by the horizontal dotted-dashed lines in the lower panels of the same figures). We refer to this frequency as the critical frequency, νcrit, and in Sects. 8 and 9 we discuss its importance in detail. Interestingly, we found that, on average, the time-lag bias is similar to its value in the absence of measurement errors at frequencies below νcrit. The same result holds for the time-lag estimates in the case of m = 10, 30, and 40 as well.
![]() |
Fig. 10 Probability distribution of the averaged phase-lag estimate for the LS40.8 light curves in experiment THRF1 at different frequencies. |
7.2. Effect of measurement errors on the sample coherence
Measurement errors have the intuitive effect of decreasing the intrinsic coherence between two time series. By construction, the intrinsic coherence of all light curves is equal to one at every frequency. As we show in Appendix C, noise decreases the intrinsic coherence at all frequencies. This decrease is minimal at frequencies where the amplitude of the measured signal’s intrinsic variations dominates over the amplitude of the noise variations. In the opposite case, the mean sample coherence will be significantly less than unity.
This is exactly what we see in the bottom panels of Figs. D.13−D.17; the sample coherence is always less than unity at high frequencies, even in the highest S/N cases. The decrease depends on the S/N; at a given frequency, the mean sample coherence decreases with decreasing S/N. By fitting various functions to the sample mean coherence we found that, for each numerical experiment and both light curve types, the following simple exponential function describes them well (see Appendix C for a theoretical justification of this fact): (25)This is demonstrated in Fig. 11, which shows the mean sample coherence for the LS40.8 light curves (m = 20) in experiment THRF1 for various S/N combinations (we get the same results in the case of experiments CD1 and PLD2 as well). The red dashed lines in this figure indicate the best-fit of the function defined above.
The blue horizontal dashed lines in the same figure indicate the value 1 /m. As we demonstrate in Appendix C, at frequencies where experimental noise dominates the intrinsic variations, the mean sample coherence is expected to be equal to that value. This is indeed what we observe in our simulations, as in the case of the lower S/N light curves the sample coherence indeed converges to the value of 1 /m (see the top panel in Fig. 11).
![]() |
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). |
8. The error of the time-lag estimates
The top panels in Figs. D.13–D.17 show the mean sample coherence for the LS40.8 and LS20.4 light curves and various S/N combinations. Columns 1–4 show the results for m = 10, 20, 30, and 40, respectively. The plots show the mean sample coherence for all experiments where we considered the effects of noise (CD1, PLD2, and THRF1). They are difficult to identify, because their values are practically identical for all S/N combinations. In the middle panels, we plot the so-called error ratio, which we define as , where
is the mean analytical error estimate (computed using Eq. (22) by substituting the value of the intrinsic coherence by its estimate), and
is the observed 1σ scatter of the averaged time-lag estimates around their mean. The bottom panels in the same figures show the probability
, i.e. the probability that the value of the time-lag estimate falls within the range
(i.e. within 1σ of the sample mean).
The black, red, and green lines in the same figures indicate the results obtained from experiments CD1, PLD2, and THRF1, respectively. Continuous and dashed lines correspond to the LS40.8 and LS20.4 light curves. For a given m and S/N combination, the results are almost identical for the different experiments and light curve types. This indicates that the error of the time-lag estimates depends mainly on m and the S/N combination, irrespective of the intrinsic CS and light curve duration.
The ratio of the analytical error over the observed time-lag standard deviation is larger than ≈ 0.9 and approaches unity as m increases, at frequencies which are lower than a certain “critical” frequency, say νcrit. This frequency is indicated by the vertical dashed lines in all panels of Figs. D.13−D.17. In the same frequency range, the probability is close to 0.68 (this value is indicated by the horizontal dotted lines in the bottom panels of the same figures). This is what we would expect if the distribution of the time-lag estimates is approximately Gaussian.
The critical frequency is independent of the intrinsic CS and on the light curve duration, and depends mainly on the light curve S/N and on m; νcrit increases with increasing m and S/N. We found that νcrit can be estimated by equating the mean sample coherence to the value 1.2/(1 + 0.2m), i.e. (26)This coherence value is indicated by the blue dotted-dashed horizontal line in the top panels of Figs. D.13−D.17.
At higher frequencies, the analytical error underestimates the true scatter around the mean of the time-lag estimates, and so the error ratio decreases. As we show in Appendix C, the mean sample coherence is always larger than the coherence of the “noisy” processes at all frequencies. Since the analytical time-lag error estimate depends on the sample coherence (see Eqs. (21) and (22) and the comment at the end of Sect. 6), it is not surprising that it underestimates the true scatter around the mean, and that the error ratio is smaller than unity. The mean sample coherence becomes significantly larger than the coherence of the “noisy” processes at frequencies ≳ νcrit. Consequently, the analytical error estimate begins to increasingly underestimate the true scatter around the mean and hence the error ratio decreases at the same frequencies.
Interestingly, at frequencies ≲νcrit the probability also remains constant and close to 0.68 (i.e. the value expected for a Gaussian distribution). At higher frequencies it begins to decrease, indicating that the distribution of the time-lag estimates develops longer “tails” than would be expected from a Gaussian distribution. We investigate the properties of the probability distribution of the time-lag estimates in more detail below.
9. The probability distribution of the time-lag estimates
The results presented above regarding the error of the time-lag estimates imply that, at frequencies below νcrit, their distribution may be, approximately, Gaussian. To further investigate this issue, we estimated the excess kurtosis and skewness of the time-lag estimates distribution, for the three experiments and two light curve types. Figures D.18–D.22 show our results (upper and middle panels for the excess kurtosis and skewness, respectively). Black, red, and green lines indicate the results obtained from experiments CD1, PLD2, and THRF1, respectively, while continuous and dashed lines correspond to the LS40.8 and LS20.4 light curves. As before, the lines in each panel overlap, indicating that the distribution of the time-lag estimates depends mainly on m and the S/N, and not on the intrinsic CS or light curve duration.
The vertical lines in all panels indicate νcrit, estimated as explained in the previous section. For a Gaussian distribution, the excess kurtosis and skewness should both be equal to zero. At frequencies lower than νcrit, this is roughly the case for the distribution of the time-lag estimates when m ≳ 20. In fact, the skewness is ≈ 0 at all frequencies, although the scatter of the sample skewness increases at frequencies higher than νcrit. This result indicates that the distribution of the time-lag estimates are symmetric around their mean. The excess kurtosis on the other hand is significantly different from zero at high frequencies, asymptotically reaching the value of −1.2. This is the excess kurtosis of the uniform distribution.
The lower panels in the same figures show the probability, pKS(νk), that the time-lag distribution is Gaussian, with a mean and variance equal to the sample mean and variance of the distribution at each frequency. This probability was estimated using the Kolmogorov-Smirnov (KS) test. The dotted lines in the bottom panels show pKS(νk) = 0.01. This is the typical threshold probability that would normally be considered if one wanted to reject the hypothesis of a Gaussian distribution for the time-lag estimates. The vertical line in the bottom left panels of Figs. D.18−D.22 show that this probability is higher than 0.01 at frequencies below the critical frequency, even in the case of m = 10.
As we discussed above, at frequencies ≲νcrit the analytical time-lag error estimate differs by ≲10% from the true scatter around the mean. The reason for this effect can be explained by the properties of the mean sample coherence (as discussed in Appendix C). It is not easy to explain why the distribution of the time-lag estimates should approach a Gaussian in the same frequency range (i.e. below νcrit), although our results indicate that this is the case.
Obviously, the time-lag distribution is not exactly Gaussian. In the case of the m = 10 and 20 averaged time-lag estimates, the excess kurtosis is certainly larger than zero at all frequencies, although not dramatically so at frequencies ≲νcrit. Likewise, pKS(νk) is rather low, but not less than 0.01 in the same frequency range. Perhaps the most interesting result for many practical applications is the fact that the ± 1 analytical error corresponds to 68% of the time-lag distribution when m ≳ 20.
10. Discussion
Our investigation was based on both analytical work and extensive numerical simulations. The simulations consisted of generating intrinsically coherent artificial light curve pairs with a prescribed model PSD and CS. The case of unevenly sampled light curves can be treated with maximum likelihood methods (e.g. Miller et al. 2010b,a; Zoghbi et al. 2013a), which we will address in a future work.
The artificial light curve characteristics (i.e. duration, sampling rate and time bin size) resembled those offered by AGN observations with present and past X-ray satellites. We assumed the same model PSD for all light curves, with a RMS amplitude of A = 0.01 and a power-law shape that smoothly “bends” from a slope of −1 to −2 after a “bend-frequency” νb = 2 × 10-4 Hz. This particular shape is characteristic of AGN X-ray PSDs. The values of A = 0.01 and νb = 2 × 10-4 Hz correspond to the mean values determined by fitting the X-ray PSDs of a sample of ≈ 100 nearby AGN (González-Martín & Vaughan 2012).
The phases of the model CS correspond to phase-lags that are commonly assumed between AGN X-ray light curves in different energy bands. They are divided into three main categories; constant delays of 10, 150, and 550 s, power-law time-lag spectra with an amplitude of 10-3−10-1 s and negative index 0.5−1.5, and time-lag spectra expected in a simple reverberation scenario where the light curves are related by a top-hat response function characterised by a centroid of t0 = 200−2000 s, width Δ = 200−2000 s and f = 0.2. The constant delays correspond to a light-crossing time of ≈ 1−100rg and ≈ 0.1−10rg for a black hole of mass MBH = 106M⊙ and 107M⊙, respectively (rg = GMBH/c3 ≈ 5(MBH/ 106M⊙) s is the gravitational radius of a black hole of mass MBH). Power-law delays with an amplitude of ≈ 10-3 s and a negative index of ≈ 1 are those typically observed between AGN X-ray light curves in different energy bands (e.g. Papadakis et al. 2001; Emmanoulopoulos et al. 2011). In addition, modelling of observed time-lag spectra with a top-hat response function requires t0 ≈ Δ ≈ 100−300 s (for AGN with MBH ≈ 106 M⊙) and f ≈ 0.1−0.3 (e.g. Zoghbi et al. 2011; Emmanoulopoulos et al. 2011). Despite the fact that the light curve characteristics and model time-lag spectra that we considered are directly applicable to AGN studies, most of our results, which we summarise below, should be relevant to all time-lag studies where evenly sampled time series are considered.
10.1. Light curve sampling and binning effects
The time-lag estimates are affected by aliasing and binning of the observed time series. If the observed signal is the result of regularly recording the values of a continuous underlying processes at regular intervals Δt, then time-lags should be estimated at frequencies lower than ≈ νNyq/ 5 (where νNyq ≡ 1/2Δt is the Nyquist frequency). This is because aliasing has the effect of decreasing and increasing the imaginary and real parts of the cross-periodogram, respectively. As a result, their ratio (which determines the phase- and time-lag estimates) is decreased, on average, and the time-lag bias increases.
If the observed series are binned (instead of sampled) over time intervals Δt, then time-lags should be estimated at frequencies ≲νNyq/ 2. This is because binning suppresses the aliasing effect on the cross-periodogram. Our results are consistent with the work of Crary et al. (1998), who reported similar findings in the context of XRB studies.
The model PSDs and CS we considered decrease with increasing frequency, and hence aliasing affects the time-lag bias more severely at increasingly higher frequencies. This is, however, expected to be true regardless of the intrinsic CS, as long as the measured signal is a (stationary) random process with finite variance (which is equal to the integrated PSD over all frequencies).
10.2. Phase-flipping effects
Phase-flipping occurs when the intrinsic phase-lag spectrum exceeds the boundaries of the interval (−π,π ]. In such cases, the phase-lag estimate will jump from π to −π (or vice-versa). This can severely alter the shape of the time-lag spectra; intrinsically constant or power-law-like time-lag spectra can show broad, prominent, “oscillatory” features. In theory, the transitions from −π to π, or vice versa, should be sharp, and hence easy to detect. However, in practice this is not the case, as the time-lag estimates will be significantly biased in the vicinity of frequencies where phase-flipping occurs. The origin of this bias can be traced to the fact that the “wings” of the phase-lag estimate’s distribution which exceed the interval (−π,π ] are folded back into the allowed range. This causes the mean of the distribution to shift towards zero. The magnitude of this bias cannot be predicted a priori, since it depends on both the intrinsic “width” of the phase-lag estimate’s distribution, as well as the unknown intrinsic phase-lag spectrum.
10.3. Smoothed vs. averaged estimates
It is customary to estimate the cross-periodogram from a single, long time series and then bin it over m neighbouring frequencies (a process called smoothing), or from m individual segments of identical duration and then bin the resulting cross-periodograms at each frequency. This process is necessary to reduce the error of the time-lag estimates, and is in fact necessary to predict their error. This requires an estimation of the coherence function itself (if we estimate the coherence from the “raw” cross-periodogram it will be equal to unity at all frequencies, irrespective of the intrinsic coherence; see P81).
Smoothing of the cross-periodogram can potentially introduce a serious bias to the time-lag estimates. For example, the model CS we considered have a power-law dependence on frequency. This can result in biased estimates of the real and imaginary part when we perform a linear type of smoothing. As a result, in some of the numerical experiments we performed, the time-lag estimates converge to zero at low (≲2 × 10-4 Hz) frequencies. Since the time-lag bias owing to smoothing originates from the cross-periodogram bias, it can only be predicted by prescribing a model CS (and not just a model time-lag spectrum). We therefore suggest the cross-periodogram to be averaged over individual segments, rather than smoothed, before computing time-lag estimates.
10.4. Measurement error effects
In the presence of measurement errors, the coherence decreases at all frequencies. This effect becomes more severe at frequencies where the amplitude of noise variations dominates over the amplitude of the measured signal’s intrinsic variations, and the coherence tends to zero. In the case when the intrinsic coherence (i.e. the coherence in the absence of measurement errors) is unity, we found that the mean sample coherence is well fitted by Eq. (25). The parameters ν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 time-lag estimates to become increasingly uniform at high frequencies. Their mean converges to zero, and their scatter is poorly approximated by standard analytical prescriptions. This is because, as the coherence decreases with increasing frequency, the “width” of the phase-lag distribution increases. When this width is sufficiently large, the wings of the distribution exceed the interval (−π,π ], and are thus folded back into the allowed range. This causes the distribution to become increasingly uniform, and its mean to shift towards zero.
We found that measurement errors have a minimal effect on the time-lag bias at frequencies where the sample coherence is ≳ 1.2/(1 + 0.2m). Furthermore, the analytical time-lag error estimate will differ from the true scatter by ≲10% in the same frequency range, as long as m ≳ 10. In addition, if m ≳ 20, the ± 1 analytical error corresponds to 68% of the time-lag distribution, and the probability distribution of the time-lag estimates approximates a Gaussian (see Sect. 9). These results hold regardless of the intrinsic PSDs and CS.
10.5. The effects of finite light curve duration
The results mentioned in Sects. 10.1–10.4 above should be applicable in all cases, i.e. they should be independent of the intrinsic PSDs and CS. However, one potentially important issue in correctly determining the intrinsic time-lag spectrum is the cross-periodogram bias. The mean cross-periodogram is not equal to the intrinsic CS, but to its convolution with the so-called Féjer kernel (see Eq. (13)). This function approaches the Dirac δ-function in the limit T → ∞, but when T is finite the bias is non-zero. Its magnitude is difficult to predict, as it depends on the shape of the unknown intrinsic CS.
Owing to the above fact, the time-lag estimates will also be biased. We found that, for the model CS we considered, the time-lag bias is positive, in the sense that the mean time-lag estimate at a given frequency has a smaller magnitude than its intrinsic value. We also found that the relative time-lag bias is ≲15% when T ≳ 20 ks, for all model time-lag spectra we considered.
In general, we found that the time-lag bias decreases with increasing light curve duration as . However, it is difficult to determine the normalisation of this relation. For a given intrinsic CS, it depends on the frequency range on which time-lags are estimated as well as the cross-periodogram bias, i.e. the shape and amplitude of the intrinsic CS. In other words, the naive expectation of the relative time-lag bias being negligible if the light curve duration is much larger than the magnitude of the time-lag estimates is generally incorrect. For example, in experiment CD1, where the intrinsic time-lag spectrum is equal to 10 s at all frequencies, the relative time-lag bias is ≈ 20% at the lowest sampled frequency, even if we use ≈ 100 ks light curves (see e.g. the continuous black curve curve in the top left panel of Fig. D.1).
To illustrate the complexity of the time-lag bias’ dependence on the intrinsic CS, we revisited experiment CD1. The solid line in Fig. 12 shows the relative time-lag bias determined from ≈ 100 ks light curves (this line is identical to the black line in the top left panel of Fig. D.1). The dashed red line shows the expected relative time-lag bias for the same light curves and intrinsic time-lag spectrum, and an intrinsic PSD whose amplitude is larger by a factor of five. We used Eq. (13) to estimate the mean cross-periodogram, E [ Ixy(νp) ], and then determined the expected relative time-lag bias by computing { 10 s−(2πνp)-1arg { E [ Ixy(νp) ] } } /10 s at each frequency. The relative bias does not change when compared to the original PSD amplitude. We then repeated the same calculation, but assumed a “bend-frequency” of 2 × 10-5 Hz (i.e. 10 times lower than the original value). The expected relative bias (shown by the blue dashed-dotted line in the same figure) is now increased at each frequency, reaching a value of ≈ 30% at the lowest frequency.
These results can be understood by inspecting Eq. (13), which determines the cross-periodogram bias. In our simulations, the CS is given by hXY(ν) = (μY/μX)hX(ν)eiφXY(ν), where hX(ν) and φXY(ν) are the intrinsic PSD and phase-lag spectrum, respectively. Since the intrinsic PSD is given by Eq. (14), a change in the PSD amplitude has the effect of altering the amplitude of the imaginary and real parts of the intrinsic CS by the same multiplicative factor at each frequency. This however does not change their ratio, which determines the mean of the phase- and time-lag estimates, and hence leaves the relative time-lag bias unchanged. On the other hand, a change in the bend-frequency alters the amplitude of the imaginary and real parts of the intrinsic CS differently at each frequency, resulting in a change of the relative time-lag bias.
We therefore suggest that it is important to always check the expected relative time-lag bias to validate the time-lag analysis in practice. This can be done by assuming a possible model CS using Eq. (13). If the intrinsic coherence of the light curves is close to unity, such a model CS may be determined by modelling the PSDs and considering various model time-lag spectra that may be potential candidates for the intrinsic time-lag spectrum.
![]() |
Fig. 12 Relative time-lag bias for the LS102.4 light curves and different values of the model PSD parameters, { A,νb }, in experiment CD1. |
11. Summary
We investigated the statistical properties of Fourier-based time-lag estimates. These estimates are based on the cross-periodogram, which is an estimator of the intrinsic CS between two random time series. Unlike the periodogram (the traditional PSD estimator) which is always positive-definite, the cross-periodogram is generally a complex number whose phase is an estimator of the phase-lag spectrum (the time-lag is equal to the phase-lag divided by the angular frequency). Our aims were to quantify the effects of light curve characteristics and measurement errors on the time-lag bias (i.e. the difference between their mean and intrinsic values), to study their distribution, and to investigate whether the standard analytical error prescriptions are accurate approximations of their intrinsic scatter around the mean.
Based on our results, we suggest the following steps when estimating time-lags:
-
Estimate the averaged cross-periodogram (Eq. (18)), using at least m = 20 pairs of data segments.
-
Use the averaged cross-periodogram to estimate the coherence and time-lags (Eqs. (23) and (20)), along with their errors (Eqs. (24) and (22)).
-
Determine the frequency range for which the sample coherence is greater than 1.2/(1 + 0.2m). If the intrinsic coherence is close to unity6, then the coherence estimates should be well fitted by Eq. (25). By determining the best-fit ν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 time-lag estimates will be similar to its value in the absence of measurement errors, the analytical prescription for the error of the time-lag estimates will be similar to the true scatter around the mean, and their distribution will be roughly approximated by a Gaussian, in the sense discussed in Sect. 9.
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 cross-periodograms 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 count-rates assumed to construct the simulated light curve pairs.
There are many smoothing schemes for the cross-periodogram. 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 PIRSES-GA-2012-31578 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ález-Martí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 [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 cross-spectrum of two binned time series
Let us consider the continuous process { X(t),Y(t) } and the discrete process { X(ts),Y(ts) }, defined according to Eq. (7) in Sect. 3. It is straight-forward to show that the processes have the same mean values. However, this is not true for their CCFs and CS as well.
For notational convenience and without loss of generality, we will henceforth assume that the intrinsic process has a mean value of zero. By definition, the CCF of the discrete process is only defined at lags which are integer multiples of the sampling period, i.e. τl = lΔtsam, where l = 0, ± 1, ± 2,..., and is given by (A.1)where RXY(τ) ≡ E [ Y(t)X(t + τ) ] is, by definition, the CCF of the continuous process. According to Eq. (A.1) RXY(τl) and RXY(τ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Δtsam, Eq. (A.2) becomes
(A.3)since ei2πkl = 1 for all k, l. The relation between the CCF and CS of the discrete process, hXY(ν) (which is only defined on the interval | ν | ≤ 1/2Δtsam), is
(A.4)Comparing Eqs. (A.3) and (A.4), we arrive at the desired result:
(A.5)
Appendix B: The mean of the cross-periodogram
Using Eqs. (9) and (10), we may write the cross-periodogram of the time series { x(tr),y(tr) } 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 = tj − tk = sΔtsam is the lag. Therefore, Eq. (B.2) can be re-cast into the following form:
(B.3)Applying the expectation operator on both sides of Eq. (B.3), we get
(B.4)where we have ignored the effect of estimating the mean values of the observed time series to reach the last equation. RXY(s) is the CCF of the discrete processes { X(ts),Y(ts) } that was defined in Sect. 3. Since the CCF is equal to the inverse Fourier transform of the CS, Eq. (B.4) becomes
(B.5)where
(B.6)is the so-called Fejér kernel. One of the properties of this function is that, as N → ∞, FN(ν) → δ(ν) (the Dirac δ-function). Thus, as N → ∞ the right-hand side of Eq. (B.5) converges to the intrinsic CS (modified by the effects of discrete binning or sampling according to Eq. (A.5)) value at ν, hXY(ν). Therefore, the cross-periodogram is an asymptotically unbiased estimator of the CS.
Appendix C: Effects of measurement errors on the coherence
Let us consider an intrinsically continuous random process { X(t),Y(t) }, whose mean, PSDs and CS are { μX,μY }, { hX(ν),hY(ν) } and hXY(ν), respectively. Let us also assume that this process is continuously sampled at regular time intervals ts = sΔtsam, where s = 0, ± 1, ± 2,..., and binned over time bins of size Δtbin. In the absence of measurement errors, the measured signal corresponds to a single realisation of a discrete version of the intrinsic process, { X(ts),Y(ts) }, as defined by Eq. (7). If there are measurement errors (as is always the case), the observed time series constitute a realisation of the “noisy” process (C.1)where ϵX and ϵY are a “purely random” processes with zero mean, and variance
and
, respectively. In this case,
where
and
are the PSDs of the purely random processes. Equations (C.2) and (C.3) describe the well-known fact that the addition of so-called white noise to a process has the effect of adding a constant amount of “power” to its intrinsic PSD. When the measurement errors in one light curve are uncorrelated with the errors in the other, then hXY,n(ν) = hXY(ν). 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(ts),Y(ts) }). 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 hX(ν)hY(ν)/(hϵXhϵX). In the case of power-law like PSDs, which decrease with increasing frequency, as is typically the case in AGN, this ratio will then tend to zero at high frequencies.
In addition to decreasing the intrinsic coherence, measurement errors will introduce a bias in the sample coherence as well, i.e. the mean of sample coherence will not be equal to . This can be shown by following (Vaughan & Nowak 1997), who showed that
(C.5)where
is the CS estimate in the absence of measurement errors, and
is a random complex number with zero mean and variance given by
(C.6)Since
is a random complex number, it may be represented by a vector in the complex plane with a random phase distributed uniformly over the interval (−π,π ]. As such,
randomly “perturbs” both the magnitude as well as the direction of the vector
. On average this perturbation will have no net effect on either the magnitude or direction of ĥxy(νk), hence
, although it will obviously increase its variance. This implies that measurement errors will not bias the cross-periodogram (although they will bias the phase-lag estimates at high frequencies, as discussed in Sect. 7.1).
The mean sample coherence can now be determined. By first applying the expectation operator on both sides of Eq. (23), we get (C.7)Given that E [ ĥx(νk) ] ~ hX,n(νk) and E [ ĥy(νk) ] ~ hY,n(νk), if we substitute Eqs. (C.2)−(C.4) and (C.6) into Eq. (C.7), we get
(C.8)where
is defined by Eq. (C.4). The above equation holds true for any intrinsic coherence. In our case where
at all frequencies, Eq. (C.8) reduces to
(C.9)This relation is identical to Eq. (25) when
. We found that such a relation fits the mean sample coherence well in the case of experiments CD1, PLD2, and THRF1, and both light curve types we considered in Sects. 7–9. A comparison between this relation and the mean sample coherence is shown in Fig. 11 for the LS40.8 light curves in experiment THRF1. Equation (C.9) indicates that when m is large, the mean sample coherence will be equal to the coherence of the noisy process.
However, in general, at sufficiently high frequencies the ratio tends to zero because
tends to zero and is always smaller than
at all frequencies (Eq. (C.4)). In this case, the first term in the right-hand side of Eq. (C.8) will therefore converge to zero. Consequently,
(C.10)
Equation (C.10) shows that, in the frequency range where the experimental noise variations dominate over the intrinsic ones (i.e. when ), the mean sample coherence will be equal to 1 /m.
Appendix D: Figures that show our results
![]() |
Fig. D.1 Relative time-lag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative time-lag bias for the LS102.4 (black curves), LS3.2 (red curves), and OB (brown curves) light curves. Bottom row: relative time-lag 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 time-lag bias. |
![]() |
Fig. D.2 As in Fig. D.1, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). |
![]() |
Fig. D.5 Relative time-lag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative time-lag 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 time-lag bias for the LS40.8 (black curves) and LS20.4 (red curves) light curves in the case of m = 20 segments. |
![]() |
Fig. D.6 As in Fig. D.5, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). |
![]() |
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. |
![]() |
Fig. D.10 Top row: relative time-lag 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 dotted-dashed lines in the lower panels indicate the 1.2/(1 + 0.2m) mean sample coherence value (see Sects. 8 and 9 for details). |
![]() |
Fig. D.13 The mean sample coherence, time-lag error ratio and probability that the value of the time-lag 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). |
![]() |
Fig. D.18 The excess kurtosis, skewness and probability that the averaged time-lag estimate is Gaussian-distributed 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). |
All Figures
![]() |
Fig. 1 Sample mean of the real and imaginary parts of the cross-periodogram (top left and right panels, respectively), the sample mean time-lag spectrum and the relative time-lag 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.4-2 (red curve) relative time-lag bias begins to noticeably increase (see text for more details). The horizontal dotted line in the bottom right panel, and in all subsequent |
In the text |
![]() |
Fig. 2 Mean relative time-lag 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 LS-type light curves, respectively. |
In the text |
![]() |
Fig. 3 Mean relative time-lag bias, plotted as function of light curve duration for LS-type and OB light curves. |
In the text |
![]() |
Fig. 4 Relative time-lag bias for the LS- and OB-type light curves (top and bottom panel, respectively), for various durations in experiment CD1. Filled squares indicate the mean relative time-lag 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 OB-type light curves, respectively. |
In the text |
![]() |
Fig. 5 Top left panel: mean sample time-lag spectrum for the LS102.4 light curves in experiment CD3 (continuous black line). The black dashed line indicates the model time-lag spectrum. Top right and bottom panels: the probability distribution of the phase-lag estimates at the frequencies indicated by the vertical dashed lines in the top-left panel. The solid and dashed vertical line in these panels indicate the model and sample mean phase-lag value at these frequencies, respectively. |
In the text |
![]() |
Fig. 6 As in Fig. 1, for experiment PLD5. The blue dashed lines indicate the model CS (upper panels) and time-lag spectrum (lower left panel). The solid brown line in the bottom left panel indicates the model time-lag spectrum without taking the effects of phase-flipping into account. |
In the text |
![]() |
Fig. 7 Mean relative bias of the m = 20 smoothed (top panel) and averaged (bottom panel) time-lag estimates, plotted as a function of the relative bias of the non-smoothed/averaged time-lag estimates (the dashed lines show the one-to-one relation line) in various numerical experiments. |
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 time-lag spectrum (lower left panel). |
In the text |
![]() |
Fig. 9 Relative bias of the averaged time-lag 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 → ∞. |
In the text |
![]() |
Fig. 10 Probability distribution of the averaged phase-lag estimate for the LS40.8 light curves in experiment THRF1 at different frequencies. |
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). |
In the text |
![]() |
Fig. 12 Relative time-lag bias for the LS102.4 light curves and different values of the model PSD parameters, { A,νb }, in experiment CD1. |
In the text |
![]() |
Fig. D.1 Relative time-lag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative time-lag bias for the LS102.4 (black curves), LS3.2 (red curves), and OB (brown curves) light curves. Bottom row: relative time-lag 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 time-lag bias. |
In the text |
![]() |
Fig. D.2 As in Fig. D.1, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). |
In the text |
![]() |
Fig. D.3 As in Fig. D.1, for experiments PLD4 (first column) and PLD5 (second column). |
In the text |
![]() |
Fig. D.4 As in Fig. D.1, for experiments THRF1 (first column) and THRF2 (second column). |
In the text |
![]() |
Fig. D.5 Relative time-lag bias for experiments CD1 (first column), CD2 (second column), and CD3 (third column). Top row: relative time-lag 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 time-lag bias for the LS40.8 (black curves) and LS20.4 (red curves) light curves in the case of m = 20 segments. |
In the text |
![]() |
Fig. D.6 As in Fig. D.5, for experiments PLD1 (first column), PLD2 (second column), and PLD3 (third column). |
In the text |
![]() |
Fig. D.7 As in Fig. D.5, for experiments PLD4 (first column) and PLD5 (second column). |
In the text |
![]() |
Fig. D.8 As in Fig. D.5, for experiments THRF1 (first column) and THRF1 (second column). |
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. |
In the text |
![]() |
Fig. D.10 Top row: relative time-lag 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 dotted-dashed lines in the lower panels indicate the 1.2/(1 + 0.2m) mean sample coherence value (see Sects. 8 and 9 for details). |
In the text |
![]() |
Fig. D.11 As in Fig. D.10, for experiment PLD2. |
In the text |
![]() |
Fig. D.12 As in Fig. D.10, for experiment THRF1. |
In the text |
![]() |
Fig. D.13 The mean sample coherence, time-lag error ratio and probability that the value of the time-lag 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). |
In the text |
![]() |
Fig. D.14 As in Fig. D.13, for (S/N)x = 9 and (S/N)y = 3. |
In the text |
![]() |
Fig. D.15 As in Fig. D.13, for (S/N)x = 9 and (S/N)y = 9. |
In the text |
![]() |
Fig. D.16 As in Fig. D.13, for (S/N)x = 18 and (S/N)y = 3. |
In the text |
![]() |
Fig. D.17 As in Fig. D.13, for (S/N)x = 18 and (S/N)y = 9. |
In the text |
![]() |
Fig. D.18 The excess kurtosis, skewness and probability that the averaged time-lag estimate is Gaussian-distributed 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). |
In the text |
![]() |
Fig. D.19 As in Fig. D.18, for (S/N)x = 9 and (S/N)y = 3. |
In the text |
![]() |
Fig. D.20 As in Fig. D.18, for (S/N)x = 9 and (S/N)y = 9. |
In the text |
![]() |
Fig. D.21 As in Fig. D.18, for (S/N)x = 18 and (S/N)y = 3. |
In the text |
![]() |
Fig. D.22 As in Fig. D.18, for (S/N)x = 18 and (S/N)y = 9. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.