1 Introduction

Calculating Time Lags From Unevenly-Sampled Light Curves


Timing techniques offer powerful tools to study dynamical astrophysical phenomena. In the X-ray band, they offer the potential of probing accretion physics down to the event horizon. Recent work has used frequency and energy-dependent time lags as a tool for studying relativistic reverberation around the black holes in several Seyfert galaxies. This was achieved thanks to the evenly-sampled light curves obtained using XMM-Newton. Continuous-sampled data is however not always available and standard Fourier techniques are not applicable. Here, building on the work of Miller et al. (2010), we discuss and use a maximum likelihood method to obtain frequency-dependent lags that takes into account light curve gaps. Instead of calculating the lag directly, the method estimates the most likely lag values at a particular frequency given two observed light curves. We use Monte Carlo simulations to assess the method’s applicability, and use it to obtain lag-energy spectra from Suzaku data for two objects, NGC 4151 and MCG-5-23-16, that had previously shown signatures of iron K reverberation. The lags obtained are consistent with those calculated using standard methods using XMM-Newton data.

Subject headings:
AGN etc

1. Introduction

Variability is ubiquitous in astrophysical phenomena. The observed fluxes are seen to vary on different time-scales for the different classes of objects. From milliseconds in neutron stars, stellar mass black holes and Gamma Ray bursts, to progressively increasing time-scales for the different phases of planetary, stellar and galaxy evolution. Variability has been invaluable in understanding the dynamics of many systems that is otherwise observationally inaccessible.

Although a lot of the methodology relies on time-domain analysis, frequency-domain techniques remain the standard tool in characterizing time-scale dependency and (quasi-)periodicity when good data coverage is available. Stellar mass black holes and neutron star LMXBs in particular show a rich phenomenology in the frequency domain, that is very tied to state transitions that have distinct spectroscopic signatures (van der Klis, 2000; Remillard and McClintock, 2006).

In addition to estimating the power spectral density (PSD) in the frequency domain for a single light curve, other measures exist when multiple, simultaneous (e.g. at different bands or energies) light curves exist. The cross spectrum gives a measure of the combined variability power in two light curves, the coherence measures the fraction of one light curve that can be predicted from the other (Vaughan and Nowak, 1997), and the phase lag gives the relative delay (in units of radians, between and ) between the two light curves as a function of frequency (Miyamoto and Kitamoto, 1989). The phase lag is converted into a time lag (in units of seconds) by dividing the phase lag by the angular frequency of the measurement.

In AGN X-ray studies that motivated this work, frequency-domain techniques have been used extensively to characterize the broadband variability (Papadakis and McHardy, 1995; Uttley et al., 2002; Vaughan et al., 2003; McHardy et al., 2004). Periodicity is generally not seen (see Gierliński et al., 2008, for the exception), and the power spectra are characterized with a power-law of index at high frequencies that breaks or bends to at a characteristic frequency that appears to scale with mass (McHardy et al., 2006). The standard tools in this case is the Fast Fourier Transform (FFT), which, for a continuously-sampled time series, gives a set of complex numbers at specific frequencies. The periodogram, which estimates the power spectrum, is the squared amplitude of these complex transforms.

Inter-band time delays in the standard X-ray bandpass ( keV) in AGN have just started to be explored in detail. Low frequency hard lags have been seen in early XMM-Newton observations (Papadakis et al., 2001; McHardy et al., 2007; Arévalo et al., 2008), where hard bands lag softer bands and the lag magnitude depends on the separation of energies, similar to that in stellar mass black holes (Nowak et al., 1999; Kotov et al., 2001). More recently, high frequency lags have also been seen (e.g. Fabian et al., 2009; Zoghbi et al., 2010; De Marco et al., 2013; Cackett et al., 2013). In this case, the lag can be soft or hard depending on the selected energies, but it is distinguished from the low frequency lags by its energy dependence. The shape of the lag-energy spectra, which gives a measure of the inter-band delays as a function of energy at a particular Fourier frequency, appears to be closely related to the spectroscopic components in a standard spectrum (Zoghbi et al., 2011; Kara et al., 2013; Zoghbi et al., 2013). The fact that the keV band leads both the keV and keV at high frequencies points to a reflection origin, with relativistic reflection being the most plausible explanation. In this case, the reflection spectrum, matched by the lag-energy spectrum, is produced within gravitational radii from the event horizon of the black hole, and reverberation is produced when the reflecting medium responds to the fast variations of the illuminating source, providing a powerful tool to probe these environments (Reynolds et al., 1999).

Similar to power spectra, time lags in these cases are calculated from the FFT (see sec. 2) for continuously-sampled light curves. The statistical properties of lag measurements in this case are discussed in Nowak et al. (1999). Extending reverberation studies beyond XMM-Newton data is not possible using the standard Fourier techniques because of the inherent non-continuous sampling forced by the low earth orbits of other observatories like Suzaku, NuSTAR and AstroH. The lowest frequency probed with the standard Fourier techniques are those associated with the orbital period of the satellite, which for the case of low-earth orbits is higher than the frequency of the interesting reverberation in AGN. Miller et al. (2010) introduced a method based on likelihood maximization that directly fits for frequency-dependent time lag (along with the power spectrum, based on the work of Bond et al., 1998). In this work, we explicitly discuss it in detail, assessing its applicability using Monte Carlo simulations. Then, we apply it to Suzaku observations of two objects, NGC 4151 and MCG-5-23-16, that had previously shown relativistic reverberation delays in the iron K band. We start section 2 of this article by reviewing the standard Fourier techniques for both power spectra and time lags. In section 3, we describe the formalism of the likelihood method. Section 4 discusses the detailed Monte Carlo simulations of the applicability of the method. The application of the method to Suzaku archival observations of NGC 4151 and MCG-5-23-16 is presented in section 5.

2. Standard Fourier Techniques

In this section, we briefly review the standard techniques based on the Fourier transform, which are commonly used with evenly-sampled light curves. The power spectral density is a property of the stochastic process producing the variability, and it gives a measure of the variability power as a function of temporal frequency . It is estimated by calculating the periodogram . If the observed data is in the form of a vector of length that gives the count rates at times , where takes the integral values and is the time bin size, the periodogram is given by the squared amplitude of Discrete Fourier Transform of :


where with . A is a normalization factor, which we take in this work to be (Vaughan et al., 2003). The periodogram itself is an inconsistent estimator of , where its standard deviation at frequency is equal to its value (Priestley, 1981). The variance is reduced significantly if several frequencies are grouped together (e.g. Papadakis and Lawrence, 1993).

Let us consider a second light curve which gives the count rate at the same time intervals but in another energy band. The cross spectrum can be estimated as , where and are the Fourier transform of and respectively, and the is the complex conjugate of . The cross spectrum is a complex number. Its amplitude is usually expressed in the form of the coherence function (Vaughan and Nowak, 1997), where the angle brackets denote averaging. The phase of the complex cross spectrum gives the phase lag between the two light curves (Miyamoto and Kitamoto, 1989; Nowak et al., 1999):


The time lag is then obtained by dividing by , so that . gives measure of time delay between and as a function of frequency (or variability time-scale).

The above calculations require the light curve to be evenly-sampled so the Fourier transform can be utilized. If this is not the case, other techniques are needed. The following section discusses the method of using the likelihood function to directly fit for the best estimates for the power and cross spectra as well as the phase/time lags directly.

3. Likelihood Analysis

The principle idea behind the method was first presented in the context of X-ray light curves by Miller et al. (2010). Here, we expand it and show explicitly how the method works and perform Monte Carlo simulation to assess its applicability. The method fits for the most likely variability powers and time lags given the observed data. Starting with a model for the power and time lags (which can be of a functional form such as a power-law, or parameterized with the values of power and time lags in pre-defined frequency bins), a likelihood function that compares the model to the data is constructed (by comparing the auto- and cross-correlations of the data with those expected from the model), and the best estimates are obtained by maximizing this likelihood function. We start in section 3.1 by applying it to estimating the power spectrum, then in section 3.2, with a simple extension, we use the method to estimate the cross power spectrum and the time lag as a function of Fourier frequency.

3.1. Power Spectrum Estimate

As before, the light curve is taken to be with values for , but now . Following Bond et al. (1998) (although that work is for 2d CMB data), each is the sum of the contribution from the signal and noise . The noise is assumed Gaussian, which is almost always the case given that each results from binning measurements obtained at a sampling smaller than . So the observed light curve is:


with a correlation matrix given by:


where is the source signal correlation matrix and is the noise correlation matrix. The angle brackets indicate ensemble average and we have assumed that the source noise components are independent. Because the measurement errors in light curves are independent, is diagonal with entries , . In general, if the observations have correlated noise, they can be easily incorporated here by adding non-diagonal elements to .

is unknown, and its components defined by are related to the underlying power spectrum through the autocorrelation function .


using the relation that the autocorrelation is the Fourier transform of the power spectrum, and for real functions only the cosine term is included. Now starting from that depends on a number of parameters (of length say) to be found, we construct and calculate the likelihood functions for those parameter:


where the dependence on is in through its dependence on , and is the transpose of . Thus, is calculated from the model and is the data vector. The procedure now is to select a model and fit for the parameters that maximize the likelihood function in equation 6. can be taken to be a power-law or a broken power-law function of . Alternatively, we can fit the band powers directly, taking the powers in pre-defined frequency bands as the parameters . This is the best option when the intrinsic shape is unknown, which might not be the case for , but is certainly the case for frequency-dependent lags .

The standard is to maximize instead of , and because of the functional form of the likelihood, the gradient and the second derivatives of the likelihood can be calculated (Bond et al., 1998). An iterative quadratic approximation is then used to find the maximum likelihood. The structure of the log-likelihood function is relatively smooth to converge within a few iterations.

3.2. Time Lag Estimate

Extending the previous formalism to include time lags is straight forward. Now we have another light curve , that represents the count rates in a different band for example. can be appended to the vector to give an augmented data vector (Rybicki and Press, 1992). The covariance matrix of the new data vector is:


where is the covariance matrix of the second light curve defined in a similar way to equation 4, with being in general different from . The matrix is the cross-covariance matrix of . The noise components of and are assumed independent because the two light curves are produced by independent events in the two bands. The noise components are also independent of the two light curves so that , etc. Therefore .

The cross-covariance is related to the cross-power spectrum and the phase lag through the cross-correlation function :


The parameters we are interested in are now in the matrix . The likelihood equation is similar to equation 6, replacing and with and respectively. Again, we use a pre-selected model (e.g. power-law), or choose the powers and lags in pre-defined frequency bins as the parameters of interest. In this work, we choose the latter parameterization. If is the number of frequency bins, then we have parameters: the powers for each light curve, the cross powers and the phase lags. For this case, the maximization procedure starts with obtaining the PSD values for individual light curves first, then for the cross power and phase lags.

In practice, there are also other effects that need to be considered. Aliasing is a consequence of the fact that power cannot be calculated beyond the Nyquist frequency . The result is that the measured powers at a frequency contain also contribution from its aliases above . Fortunately however, X-ray light curves generally have power-law PSDs, so the power above is small, and also the fact is a width of a bin not the actual sampling time. The binning process is equivalent to convolving the light curve with a binning window for and otherwise. The result is that the is multiplied by the Fourier transform of , which is: (e.g. van der Klis, 1989). Red noise leak is another problem and it is the result of the finite length of the observation. in this case is convolved with the Fourier transform of the window function, and mainly causes power below the lowest measured frequency (, where is the length of the observation) to leak into frequencies above . One can explicitly include the convolution of the window in equations 5 and 8. However, we found that it is computationally easier to include the additional power below in the fit by extending the lowest boundary of the lowest frequency bin to values smaller than , this was found to correct for the power biases (see sections 4.1). Extending the first bin to frequencies less than assumes that the power below does not change significantly, which is a reasonable assumption given that the psd in almost all cases is a smooth power-law.

3.3. Estimating Uncertainties

As discussed in Miller et al. (2010), the uncertainties can be estimated, as part of the fitting procedure, by calculating the Fisher matrix, which is related to the second derivative of the log-likelihood (see detailed related discussion in Tegmark et al., 1997). The Fisher matrix basically measures how fast on average does the likelihood function fall off around the its maximum. When the best fit is found, an estimate of the covariance matrix of the parameters is given by the inverse of the Fisher matrix. The variance of the estimates parameters are the diagonal elements of this covariance matrix. These estimates are however only a lower limit on the uncertainties when the off-diagonal values are not small (i.e. parameters are correlated).

The alternative is to step through the parameters, taking the uncertainty as the value that changes by 1 (Miller et al., 2010). Another approach involves using Monte Carlo Markov Chain (MCMC) to map the full probability space, obtaining probability distributions for the parameters directly. The uncertainties quoted in this work, unless stated otherwise, are the result of stepping through each parameter individually, allowing the rest to change, and taking the error as the value that changes the value of by 1. This choice works when the number of parameters to be fitted is small ( , so stepping through parameters is computationally feasible relatively quickly). If the number is large, the best option is to use MCMC to obtain uncertainties.

4. Simulations

In order to test the above method, we simulate light curves with known underlying power-spectra and time delays, introduce gaps, and explore how well they can be recovered. Starting with a functional form for the , we randomize the amplitude and the phase then inverse Fourier transform to obtain one light curve realization (Timmer and Koenig, 1995). When a second light curve is needed, we shift the phase by the desired amount before performing the inverse Fourier transform. This assumes unity coherence. When fitting real data, the coherence can be estimated from the cross spectrum and the individual power spectra. Poisson noise is added to all light curves.

Figure 1.— The two cases of base power-spectra used in the simulations. Solid lines represent the underlying generating spectrum. Dashed lines include the effect of Poisson noise. Case 1 corresponds to a bright variable ( rms) source and case 2 corresponds to a relatively fainter and less variable source ( rms).
Figure 2.— Typical light curve realizations from cases 1 (panel a) and 2 (panel b) in Fig.1. In each case, the second light curve is lagged by 1 radians with respect to first. Panel c shows typical light curves with gaps for the two cases. The y-axis is similar to panels a and b.
Figure 3.— Left: The average estimated PSD for the high (a) and low-rms (b) cases without including gaps. Right: The histogram of the values for two selected frequency bands (marked with vertical lines in the left panel ). For each of the high and low-rms cases, more than 2000 separate light curves are simulated. For each one, the power is estimated at nine frequency bins. The means of the resulting distributions are plotted in the left panel. Their errors bars represent the standard deviation of the distribution. The average error from the 2000 estimates are plotted as the dotted lines above and below the best estimate. The solid line is the the value of (i.e. input) at that frequency. The frequency error bars representing the width of the bin are omitted for clarity. The horizontal lines in the histogram plots represent the mean (solid) and input model (dotted).

In this work, we take the input power spectrum to be a broken power-law of the form:


where below some break frequency , and above it, and is a normalization factor. We take Hz. This break frequency is consistent with a black hole of , which is typical of many Seyfert galaxies (McHardy et al., 2006). The lag in the simulated light curves is taken to be constant in phase (at 1 radians), so that the lag scales with . For the gaps, we tried different patterns, as will be discussed. The most relevant given X-ray observations are those which have a period of hours, typical of low earth orbit observations. Throughout the following simulations, we study two cases: and , representing high and low power cases respectively. The count rates for the two cases are 5 and 1 respectively. These are chosen as typical values for a bright variable, and relatively faint, less variable sources, corresponding to rms variability of 35 and 10 percent in each case. For each of these two cases, we run simulations with, and without gaps. The simulations without gaps are used for comparison and consistency check. All the simulated light curves are equivalent to an exposure of 200 ks. In simulations with gaps, we discuss both on-source and total exposures of 200 ks as detailed below.

4.1. Power spectra

First, we discuss estimating the power spectrum, starting with a simple, high power, high signal to noise case without including gaps to use as a proof of concept. We simulated more than 2000 light curve realizations from case-1 PSD with 1 second sampling rate, we add Poisson noise then binned the light curves to 512 second bins. We used frequency bins that give in the case of even sampling at least 10 Fourier frequency points per bin. The same experiment is repeated for the low-rms case. The model PSDs and typical light curve realizations for these two cases are shown in Fig. 1 and Fig. 2.

The result is summarized in Fig. 3, where we show the ensemble-averaged measured power at nine frequency bins, along with histogram distributions for two selected frequency bins. Each simulated light curve gives an estimate of the power spectrum at the nine frequencies and their uncertainties. The left panel of Fig. 3 shows the means of these estimates (points). The errors on those points are taken to be the standard deviation of estimates around the mean. The average of the measured uncertainties are also plotted as the dotted lines around the best estimates. The right panel shows the distribution of the 2000 values for two selected frequency bins (1st and 5th bins).

There are several points to note from Fig. 3. The power spectrum is well-recovered both for the high and low rms cases. Even in noise-dominated bands in the low rms case ( Hz, see Fig. 1 ). The noise in the light curves is accounted for automatically ( and in equations 3 and 4), so that the measured PSD is the underlying, noise-less .The estimates are nearly Gaussian, particularly at intermediate frequencies. The distribution tails at the lowest frequencies are consistent with expectations from standard Fourier analysis for bins with a small number of averaged frequencies (e.g. Papadakis and Lawrence, 1993). The shape of the distribution depends essentially on the effective number of independent frequencies present in the light curve. The central limit theorem ensures that when a relatively large number is averaged, the distribution is Gaussian. If the number is small, an estimate is obtained, the errors may not be Gaussian, but the formalism presented here allow us to also estimate probability distribution of the lags using either direct evaluation of the likelihood function, or through Monte Carlo Markov Chains.

It is also clear that the method gives unbiased, consistent estimates of the power. The plot also shows that the uncertainty estimates (dotted envelopes), discussed in sec. 3.3 and taken here as the average of individual uncertainties, are very consistent with the standard deviation of an ensemble of estimates. In fact, for the cases of Fig. 3 where no gaps are included, the frequencies are independent, and so the uncertainties taken directly from the Fisher matrix and those estimated by stepping through the parameter space are the same.

Figure 4.— Left: The average estimated power spectra (similar to Fig. 3 ) including light curves with gaps for high (Top) and low (bottom) rms cases. The points represent the distribution means for the case of no gaps (similar to Fig. 3) and the error bars are the standard deviation around the mean. The red dotted envelope is the standard deviation of the estimates including gaps with the light curves of length 200 ks (data + gaps, G1). The green dashed envelope is the standard deviation of the estimates for light curves with gaps, but with on-source exposure of 200 ks (data only, G2). Right A histogram distribution of the estimates at the 5th frequency bin ( Hz) which corresponds roughly to the periodicity of the gaps. The blue line is for light curves without gaps. The red line is for the case of total exposure of 200 ks, and the green line is for the case of on-source exposure of 200 ks.

Similar simulations were performed for the low and high rms cases (as defined in section 4 and Fig. 1) considering light curves with gaps. For comparison, we simulate light curves where the length of observation is 200 ks, and also light curves where the on-source exposure is 200 ks. The gaps are generated randomly assuming that both the length of data stream and gaps are Gaussian random variables with means 5700 and 4000 seconds respectively, and a standard deviation of 100 seconds. These gap patterns roughly resemble those usually encountered in Suzaku observations and relevant to NuSTAR and AstroH. Other gap patterns have also been explored, and the conclusion are in general the same (with the obvious change of the frequencies affected.). The result is plotted in Fig. 4.

High and low rms cases are plotted in the top and bottom panels respectively. In each case, the left plot is similar to that in Fig. 3. The points and the errors bars are for the case with no gaps for comparison. The red-dotted and green-dashed lines are the envelope of the standard deviation of the PSD estimates for light curves with gaps and light curve lengths of 200 ks (hereafter case G1) and on-source exposure of 200 ks (hereafter case G2) respectively.

The gaps have several effects compared to the continuous case. The errors are in general larger because there is less data on the whole, except for the very lowest frequencies where the errors in G2 are smaller than the no-gaps case because the requirement of a on-source exposure of 200 ks means there is more low frequency data. Also, the errors for both G1 and G2 are larger for frequencies close to the gap periodicity. The reason is that information on those frequencies are missing because of the gaps. This is a general result that we found throughout the simulations, and it shows that the periodic gaps cause the uncertainties at the frequency corresponding to the gaps periodicity ( Hz). G1 has about 60-70% less exposure and its errors are slightly larger then G2 (the difference for a single frequency band is not huge, but all frequency bins are affected). The distribution histogram for the frequency bin closes to the gaps frequencies are also plotted in Fig. 4.

4.2. Time lag

Analysis similar to that presented in section 4.1 was extended to include time lags. Fig. 2 show typical light curves pairs for the high and low rms cases defined in section 4, where for each pair, the second light curve is shifted with a phase of 1 radian. The results are presented in Fig. 5. Although the presented simulations are for the case of constant phase lag at 1 radians, we tested for other forms (e.g. constant time lag, a time lag that has functional dependence on etc.), and the results are not different from those discussed here.

Figure 5.— Similar to Fig. 3 but now showing lags instead of psd. High and low rms cases are shown in panels a (top) and b (bottom) respectively for light curves without gaps. The average estimated lag is shown as points. The standard deviation around the mean is shown in as error pars. The envelope dotted line shows the average estimated uncertainties. The right panels in each case show the (normalized) number of values histogram for the 1st (un-shaded) and 5th (shaded) frequency bins, marked with vertical lines in the left panels..

The figure, analogous to Fig. 3, shows the averaged lag calculated from an ensemble of 2000 light curve realization for the high and low rms cases without gaps. The lags are well-recovered for all frequencies for both cases. The distributions of the estimates (shown in Fig. 5) are almost perfect Gaussians. The plot also shows that the uncertainty estimates (dotted envelopes), discussed in sec. 3.3 and taken here as the ensemble average of individual uncertainties, are also consistent with the standard deviation of an ensemble of estimates. The slight difference at the noise-dominated frequencies (highest frequencies at panel b in Fig. 5) is an artifact of the simulation, where the noise-dominated parameters sometimes fail to converge, and it is therefore hard to obtain uncertainties and those are removed when estimating the average uncertainties. In practical data analysis, one would reduce the number of frequency bins to improve the signal to noise ratio.

Figure 6.— Similar to Fig. 5 but now including light curves with gaps for the high rms case. The average estimated lag for light curves without gaps is shown as points. The standard deviations around the mean are shown as error bars. The envelope dotted line (red) shows the standard deviation for light curves with gaps and length of 200 ks (G1). The envelope dashed line (green) is for the case of light curves with gaps and an on-source exposure of 200 ks (G2). The right panels shows the corresponding (normalized) number of values histogram for the 5th frequency bins, which corresponds roughly to the frequency of the periodic gaps.

Extending the analysis to light curves with gaps is again straight forward. Fig. 6. As in the case of power spectra, the errors are larger for light curves with gaps because less data is used. The lowest frequencies are not affected much because the long time-scale trends in the light curves are not affected if there are gaps on smaller time-scale. The periodic gaps have the effect of increasing the uncertainty of the measured lags at frequencies close to the gaps frequency and also its harmonics where information is missing. This, combined with the gap randomness (i.e. it is not a single frequency) and frequency binning produces the fluctuations seen in Fig. 6. The results for the low rms case is very similar. The low signal to noise ratio however means the errors are larger, and sometimes simulations not constrained. Better estimates are obtained when using less frequency bins (i.e. improving the signal per bin), and in this case, the results are similar to those of the high rms case.

Figure 7.— Similar to Fig. 6 but with a different gap pattern. The gaps now have a periodicity that corresponds to a frequency of Hz. The histograms are for the second frequency point.

This increased uncertainty at the gaps periodicity is further illustrated in Fig. 7, which is similar to Fig. 6 but for a different gap pattern. Here, the gaps have a periodicity corresponding to a frequency of Hz. Again, the effect of the gaps is that less information is available to the gap frequency, and therefore the uncertainty is larger. The distribution of the estimates is Gaussian or very close to Gaussian in most cases. The power of the likelihood method presented here is that, even in frequency bands where the effective number of independent frequencies is small, one can obtain a direct measure of the probability distribution whether by stepping through the likelihood function, or more efficiently by using MCMC.

5. Applications

In this section, we discuss the application of the above method to calculate time lags in Suzaku observations of two sources: NGC 4151 and MCG-5-23-16. Time delays have been seen in this two objects using XMM-Newton and standard lag calculation methods (Zoghbi et al., 2012, 2013). The observations used in the following discussion are summarized in table 1.

Object Obs. ID Exposure (ks) Date
NGC 4151 701034010 125 08-2008
906006010 60 04-2011
906006020 60 04-2001
MCG-5-23-16 700002010 95 12-2005
Table 1Observations of NGC 4151 and MCG-5-23-16 from the Suzaku archive used in this work

Data were retrieved from the archives and reduced using heasoft 6.13 and the latest calibration files (caldb version 20130305). Cleaned events files for all XIS detectors operational during the observations were produced following the Suzaku user guide. Then xselect was used to extract source and background light curves with time bins of 512 seconds. The source region in each case was circular with radius of 3.5 arcmin, and background regions were selected from source-free regions on the CCD. In order to study time lags, we extracted light curves in eight energy bins between keV in steps of 1 keV. Background light curves were than scaled to match the area of the source region before subtracting them from the source light curves. The XIS0 and XIS3 counts were combined to produce a total front-illuminated light curve, while the XIS1 gives a back-illuminated light curve. The two light curves can then be fitted simultaneously using the formalism discussed in section 3.

Figure 8.— Lag-energy plot for NGC 4151 (left) and MCG-5-23-16 (right). The points represent the estimated lag at frequencies Hz (NGC 4151) and Hz (MCG-5-23-16) between the marked energy band and the the whole 2-10 keV band taken as a reference (after removing the band of interest to keep the noise uncorrelated. See Zoghbi et al. 2012 for details on lag-energy plots). The shaded plots are the lag-energy plots from the XMM-Newton data (Zoghbi et al., 2012, 2013)

Fig. 8 shows the lag-energy plots for NGC 4151 (left) and MCG-5-23-16 (right), along with plots from previously published XMM-Newton data (Zoghbi et al., 2012, 2013). Each point in the plots is a result of maximizing the likelihood function for the power spectra and phase lags between individual light curves and the total keV light curves, excluding the current energies (see Zoghbi et al., 2012, for details on lag-energy plots). The plotted frequencies are Hz and Hz respectively. Although the uncertainties at the highest energies are relatively large, it is clear that there is a structure at keV consistent that seen in the XMM-Newton data. The match between the Suzaku and XMM-Newton plots in the case of MCG-5-23-16 (Fig. 8-right) is remarkable. For the case of NGC 4151, although the shapes are statistically consistent, the lag-energy shape in this source is known to be flux- and frequency-dependent (Zoghbi et al., 2012). The length and quality of the Suzaku observations do not allow direct comparison at the same exact frequencies, but the fact that there is a peak at keV adds further evidence that the iron line is responsible for these lags. A further test is achieved by adding artificial gaps to the XMM-Newton and calculating lags. This however reduces the amount of available data and smears the signals out.

Figure 9.— Probability densities for the MCG-5-23-16 two lag points at and keV shown in Fig. 8 estimated using Monte Carlo Markov Chain.

The likelihood method allows us to obtain full probability densities for lag estimates and hence quantify directly the significance of any lag detection. For example, Fig. 9 shows the probability density of the estimated lag values at and keV for the case of MCG-5-23-16, plotted in the right panel of Fig. 8. After the best estimates are obtained by maximizing the likelihood function, the multi-dimensional parameter space is mapped out using MCMC. The chains were generated with an affine-invariant ensemble sampler (Goodman and Weare, 2010; Foreman-Mackey et al., 2013), and the result is a probability density for each of the estimated values.

6. Summary

We presented and discussed a method to calculate frequency-dependent power spectra and time lags for unevenly-sampled data. The method, first introduced by Miller et al. (2010), relies on likelihood maximization, and gives the most likely power and lag estimates given the data. We tested the method using Monte Carlo simulations, and showed that the main effect of periodic gaps, typical of low-earth orbit X-ray observations, is to give unconstrained estimates at the frequency corresponding to the gap periodicity, while information at other frequencies is recovered. We applied the method to Suzaku archival observations of NGC 4151 and MCG-5-23-16, and showed that their lag-energy spectra are consistent with those observed using XMM-Newton, giving further support to their interpretations as being due to relativistic reverberation close to the black holes.


  1. Arévalo, P., McHardy, I. M., and Summons, D. P.: 2008, mnras 388, 211
  2. Bond, J. R., Jaffe, A. H., and Knox, L.: 1998, prd 57, 2117
  3. Cackett, E. M., Fabian, A. C., Zogbhi, A., Kara, E., Reynolds, C., and Uttley, P.: 2013, apjl 764, L9
  4. De Marco, B., Ponti, G., Cappi, M., Dadina, M., Uttley, P., Cackett, E. M., Fabian, A. C., and Miniutti, G.: 2013, mnras 431, 2441
  5. Fabian, A. C., Zoghbi, A., Ross, R. R., Uttley, P., Gallo, L. C., Brandt, W. N., Blustin, A. J., Boller, T., Caballero-Garcia, M. D., Larsson, J., Miller, J. M., Miniutti, G., Ponti, G., Reis, R. C., Reynolds, C. S., Tanaka, Y., and Young, A. J.: 2009, nat 459, 540
  6. Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J.: 2013, pasp 125, 306
  7. Gierliński, M., Middleton, M., Ward, M., and Done, C.: 2008, nat 455, 369
  8. Goodman, J. and Weare, J.: 2010, Communications in Applied Mathematics and Computational Science 5(1), 65
  9. Kara, E., Fabian, A. C., Cackett, E. M., Miniutti, G., and Uttley, P.: 2013, mnras 430, 1408
  10. Kotov, O., Churazov, E., and Gilfanov, M.: 2001, mnras 327, 799
  11. McHardy, I. M., Arévalo, P., Uttley, P., Papadakis, I. E., Summons, D. P., Brinkmann, W., and Page, M. J.: 2007, mnras 382, 985
  12. McHardy, I. M., Koerding, E., Knigge, C., Uttley, P., and Fender, R. P.: 2006, nat 444, 730
  13. McHardy, I. M., Papadakis, I. E., Uttley, P., Page, M. J., and Mason, K. O.: 2004, mnras 348, 783
  14. Miller, L., Turner, T. J., Reeves, J. N., Lobban, A., Kraemer, S. B., and Crenshaw, D. M.: 2010, mnras 403, 196
  15. Miyamoto, S. and Kitamoto, S.: 1989, nat 342, 773
  16. Nowak, M. A., Vaughan, B. A., Wilms, J., Dove, J. B., and Begelman, M. C.: 1999, apj 510, 874
  17. Papadakis, I. E. and Lawrence, A.: 1993, mnras 261, 612
  18. Papadakis, I. E. and McHardy, I. M.: 1995, mnras 273, 923
  19. Papadakis, I. E., Nandra, K., and Kazanas, D.: 2001, apjl 554, L133
  20. Priestley, M. B.: 1981, Spectral analysis and time series / M.B. Priestley, Academic Press, London ; New York :
  21. Remillard, R. A. and McClintock, J. E.: 2006, araa 44, 49
  22. Reynolds, C. S., Young, A. J., Begelman, M. C., and Fabian, A. C.: 1999, apj 514, 164
  23. Rybicki, G. B. and Press, W. H.: 1992, apj 398, 169
  24. Tegmark, M., Taylor, A. N., and Heavens, A. F.: 1997, apj 480, 22
  25. Timmer, J. and Koenig, M.: 1995, aap 300, 707
  26. Uttley, P., McHardy, I. M., and Papadakis, I. E.: 2002, mnras 332, 231
  27. van der Klis, M.: 1989, in H. Ögelman and E. P. J. van den Heuvel (eds.), Timing Neutron Stars, p. 27
  28. van der Klis, M.: 2000, araa 38, 717
  29. Vaughan, B. A. and Nowak, M. A.: 1997, apjl 474, L43
  30. Vaughan, S., Edelson, R., Warwick, R. S., and Uttley, P.: 2003, mnras 345, 1271
  31. Zoghbi, A., Fabian, A. C., Reynolds, C. S., and Cackett, E. M.: 2012, mnras 422, 129
  32. Zoghbi, A., Fabian, A. C., Uttley, P., Miniutti, G., Gallo, L. C., Reynolds, C. S., Miller, J. M., and Ponti, G.: 2010, mnras 401, 2419
  33. Zoghbi, A., Reynolds, C., Cackett, E. M., Miniutti, G., Kara, E., and Fabian, A. C.: 2013, apj 767, 121
  34. Zoghbi, A., Uttley, P., and Fabian, A. C.: 2011, mnras 412, 59
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description