## Abstract

Progress in astronomy comes from interpreting the signals encoded in the light received from distant objects: the distribution of light over the sky (images), over photon wavelength (spectrum), over polarization angle and over time (usually called light curves by astronomers). In the time domain, we see transient events such as supernovae, gamma-ray bursts and other powerful explosions; we see periodic phenomena such as the orbits of planets around nearby stars, radio pulsars and pulsations of stars in nearby galaxies; and we see persistent aperiodic variations (‘noise’) from powerful systems such as accreting black holes. I review just a few of the recent and future challenges in the burgeoning area of time domain astrophysics, with particular attention to persistently variable sources, the recovery of reliable noise power spectra from sparsely sampled time series, higher order properties of accreting black holes, and time delays and correlations in multi-variate time series.

## 1. Introduction

Pick up any good book on time-series analysis and you will probably find that the topics covered include: trends and seasonal components, autoregressive and moving average processes (ARMAs), forecasting, Kalman filters, state–space models, spectral methods and so on. Many of these will sound unfamiliar to astronomers, even those who regularly examine time series in the form of light curves of variable astronomical objects. From the astronomer’s perspective time-series analysis, as a subject itself, appears to be a branch of statistics but with its own specialist jargon and approach, developed in no small part from applications in signal-processing engineering and econometrics. Astronomers generally do not have much formal training in statistics (if any), tend to have their own jargon, and read papers and textbooks only by other astronomers. The exchange of ideas between time-series experts and astronomers with data to analyse is not as effective as it could be. The result is, at best, inefficiency; astronomers spend time struggling with the details of unfamiliar methods or re-inventing techniques that have already been discussed extensively elsewhere. Even between different sub-branches of astronomy (e.g. between X-ray astronomy, optical and radio astronomy; or stellar and extragalactic astronomy) there may be differences of convention and practice (in data collection, curation, analysis and interpretation) that hamper effective communication and transfer of ideas.

As well as the ‘sociological’ barriers that divide astronomers from time-series analysts elsewhere, there are more practical reasons that may contribute to their separation. Astronomers often have to deal with limitations that may not trouble data collection in other fields. The most obvious problem is breaks in observations owing to observing conditions related to, for example, visibility of the target on the sky, bad weather, or telescope scheduling constraints, resulting in discontinuous and often uneven time sampling. There may be a variable background, measurement errors or instrumental effects (such as ‘deadtime’ or readout cycles) that are beyond experimental control, but must be understood through careful instrumental calibration. Another difference between much of conventional time-series analysis and astronomy is their objectives. Much of the time-series literature is concerned with modelling time series for the purpose of forecasting (e.g. in meteorology, econometrics and finance, engineering applications). Astronomers are (usually) more interested in developing and testing physical models, to better understand the physical principles at work in a variable system, rather than in predicting future values of some time-varying function for their own worth.

The main purpose of this article is, in the spirit of this Theme Issue, to facilitate discussion and the exchange of ideas between those who work on astronomical problems involving variable objects and those who work on methods for analysing time-series data. If, as a result, a few astronomers go on to consider alternative methods for analysing their data, or a few time-series specialists from the statistics or signal-processing communities get involved with astronomical data, I will consider this to have been a success.

## 2. Why time-series analysis matters to astronomers?

Astronomy is an observational science. Unlike experimental sciences, we are not able to interact with the objects we study, alter experimental conditions or repeat experiments under the same conditions to obtain more data. The information we can gather from astronomical objects (planets, stars, galaxies, etc.) is usually limited to the light we receive from them. (I will not deal here with other types of radiation, e.g. neutrino, cosmic ray and gravitational wave astronomy. The latter topic is discussed elsewhere in this Theme Issue.) In principle, we are able to examine how this light is distributed in terms of position on the sky (images), wavelength or photon energy (spectroscopy), polarization, and time. Early examples of astronomical time series include observations of the positions of the stars and planets, the phases of the Moon and planets, and sunspot numbers. The Wolf^{1} sunspot data are a staple of introductory time-series analysis books.

Most of the variable sources are effectively point-like on the sky, and most detectors are not sensitive to polarization, and so we will be concerned here with datasets that have at most two dimensions, characterizing the brightness (often called *flux* by astronomers) against time and wavelength. Astronomers usually use the term *light curve* to mean estimates of the brightness of a target at different times, i.e. a univariate time series of the target’s brightness, although other kinds of time series are also common in astronomy. Multi-variate time series are generated by instruments that have some intrinsic energy resolution (ability to distinguish photons of different energy or wavelength), or by combining simultaneous data from several observatories that are sensitive in different wavebands (e.g. X-rays, optical and radio).

At bottom, all of time-series analysis is concerned with time-varying signals that contain a non-trivial element of randomness. Depending on the target, the astronomer may be interested in deterministic ‘signals’ that are ‘contaminated’ by random fluctuations. Such signals might be transients—such as bursts of X-rays from the surface of a neutron star, of high-energy radiation from gamma-ray bursts (GRBs), or a supernova—or periodic signals such as from rotating stars, pulsars (radio pulsations because of neutron star rotation), or modulations owing to the periodic orbits of binary stars, or planets. In these cases, the astronomer is usually interested in recovering the deterministic component and testing models or estimating parameters, e.g. burst luminosity and decay time, rotation period, etc. In other cases, the ‘noise’ itself may represent the fundamentally stochastic output of an interesting physical system, as in turbulent accretion flows around black holes. Here, the astronomer is interested in comparing the statistical properties of the observations with those of different physical models, or using the intrinsic luminosity variations to ‘map out’ spatial structure. These projects, and many others, are completely dependent on time-series data and analysis; our only access to the properties of physical interest is through their signature on the time variability of the light we receive.

Many of the different types of variable astronomical objects, and the types of variability they exhibit, are introduced in the nice review by Feigelson [2] given at a conference on Time Series Analysis in Astronomy and Meteorology, which took place in 1993.^{2} In the two decades since that meeting the quality and quantity of astronomical time series has increased dramatically. In the field of X-ray astronomy, NASA’s *Rossi X-ray Timing Explorer* (RXTE) produced thousands of datasets containing high-frequency timing data (time scales as short as millisecond) from the proportional counter array (PCA), and 15 years of daily monitoring of the X-ray brightness of X-ray binaries with the All-Sky Monitor (ASM). The *Swift* mission has captured gamma-ray, X-ray and optical data of hundreds of GRBs and other high-energy transients. Optical astronomy has been transformed by large ground-based telescopes (the 8m class) and the now ubiquitous CCD detectors, and can routinely provide high-precision (in both absolute time and brightness) light curves of variable targets. Powerful optical time domain surveys are now being undertaken by dedicated robotic telescopes in the search for timing signals of, for example, extra-solar planets, supernovae, trans-Neptunian objects, dark matter candidates, producing unprecedented quantities of time domain data. For example, the *CoRoT* and *Kepler* space missions [3,4] and the MACHO, WASP and Pan-STARRS ground-based projects [5]. Radio astronomy is currently enjoying something of a renaissance, and the Dutch-led LOFAR observatory will be sensitive to radio transients on time scales from milliseconds to months [6].

When writing an article on such a broad subject the most difficult question is not what to include, but what to leave out. What follows is rather parochial, being focused on just a few of the active areas and challenges that I have some direct experience with. In particular I will mostly focus on persistent variability, not transients. Although some of the most interesting astrophysical objects—such as supernovae [7], pulsar flares [8], GRBs [9] and stellar capture events [10]—are discovered through their transient emission, such sources are the exclusive subject of another Royal Society Discussion Meeting this year (April 2012). Some of the other challenges in time domain astronomy that I have no room to discuss include high-precision pulsar timing [11], astroseismology [12] and planet detection [13,14].

## 3. Fast variability

Apart from a few rare exceptions (see below) most astronomical observations of particular targets are short, lasting minutes or hours. Observations with ground-based telescopes (e.g. in the radio, infrared (IR) or optical bands) are limited by, for example, the length of a night or the visibility of the target on the sky. Earth-orbiting telescopes may have their observations interrupted as the Earth periodically blocks the telescope’s view of the target, or because of gaps in the telemetry. For both types of observatory the observations are limited by other scheduling and operational constraints. As a result, astronomers tend to treat within-observation and between-observation variability slightly differently. Variability within a single, continuous observation is usually limited to time scales few hours. I shall refer to this within-observation variability as *fast* variability. Time series covering much longer time scales can be built up through repeat observations of the same target, often through dedicated *monitoring campaigns*, and can, in principle, access variations on time scales of ∼few hours–decades but with generally irregular time sampling. Below, I shall refer to such inter-observation variability as *slow* variability.

I will discuss some examples of fast variability from the field of X-ray astronomy. But, first, I should explain what kind of data X-ray astronomy deals with. Data from observations of particular targets by orbiting X-ray detectors are usually packaged as *event lists*, essentially just tables detailing the properties of detector events (most, but not all, of which are triggered by genuine X-ray photons). Each event is associated with several properties that give the arrival time, energy channel and position on the detector, as well as other information. If we are interested in the variability of the target, we form a light curve by selecting a subset of the events^{3} and making a histogram of the event times. See the excellent book by Arnaud *et al*. [15] for more on X-ray detectors, data and processing.

Figure 1 shows example light curves for a bright, variable X-ray source, 4U1608-52, based on observations by the RXTE. The source is an *X-ray binary*, a type of binary star system in our Galaxy comprising a low-mass ‘donor’ star transferring matter onto a neutron star via a disc-shaped accretion flow (figure 2). Figure 1*a*,*b* shows the *slow* variability recorded by the ASM instrument, shown as ≈ daily averages. Figure 1*c* shows the light curve of a single observation lasting less than 1 h with the PCA instrument. The PCA data were recorded with a time resolution of 2^{−16} s but are shown ‘binned’ to a resolution of *Δt*=1 s. The ASM data show accretion *outbursts*, where the source becomes X-ray bright, reaching luminosities of W ( solar luminosities), typically lasting days–weeks (figure 1*a*,*b*). Outbursts are thought to be driven by large-scale instabilities in the accretion flow [16]. During outbursts the source shows both persistent and transient variability. Figure 3*b* illustrates a fast transient phenomenon—an *X-ray burst* lasting about a minute caused by runaway thermonuclear ‘burning’ on the surface of the neutron star (see [17] for a review).

### (a) The power spectrum of fast variations

Figure 1 clearly shows a wide range of variability types. Here, we will concentrate on the persistent, fast variability. This may not be so obvious from figure 1 because of the choice of time binning and the appearance of the X-ray burst. Figure 4 shows a power spectrum taken from a series of observations of the same source during which the source was bright and there were no X-ray bursts. As in many areas of the physical sciences, the power spectrum is often used by astronomers to characterize the variations in more detail. It describes the distribution of *power* (squared variability amplitude) as a function of temporal frequency (∼ timescale^{−1}).

In the case of the 4U1608-52 data, the power spectral estimate shows some weak low-frequency power, a flat spectrum Hz caused by Poissonian fluctuations in the X-ray count rate, and a narrow bump at Hz. This excess of power at high frequencies is a *quasi-periodic oscillation* (QPO), a concentration of variability power in a narrow range of frequencies. This is distinct from ‘strictly periodic’ variability (repeating perfectly after a fixed interval *T*) that would appear as unresolved narrow peaks in the spectrum, and ‘noise’ where power is spread over a wide range.^{4} The peak in the 4U1608-52 power spectrum is an example of a ‘kHz QPO’ (because their frequencies range to Hz in some cases) and have been detected in approximately 20 different neutron star X-ray binaries [18,19]. They are of great interest to X-ray astronomers, being the fastest (persistent) variations of any known astronomical source, and thought to be a result of the relativistic gravity acting on the inner accretion disc close to the neutron star surface [20] since the orbital period close to the neutron star surface is approximately 1 ms.

### (b) Estimating the power spectrum

How do astronomers estimate the power spectrum from these data? The standard procedure is to follow the Bartlett method [21]:

— from the event list form an evenly sampled time series

*x*_{t}with a bin time size*Δt*;— break this into

*M*non-overlapping intervals of length*N*;— compute the periodogram of each interval;

— average the

*M*periodograms.

If the (binned) light curve is *x*_{t} and *X*_{j} is its (complex valued) discrete Fourier transform (DFT) at Fourier frequency *f*_{j}, then the periodogram is where *A* is a normalizing term [22] and *X** is the complex conjugate of *X*. Under simplifying assumptions (e.g. that the data are a realization of a linear, stationary and ergodic process) the Bartlett estimator should converge to the true power spectrum *S*(*f*) in the limit of large *N* and *M*. See van der Klis [23] for a discussion of the standard power spectral methods as applied in X-ray astronomy. This has become the standard method of power spectral estimation in astronomy when the data are evenly sampled (below we mention the case of uneven time sampling). The result is a power spectral estimate at each of the *N*/2 Fourier frequencies *f*_{j}=*j*/*NΔt*, which should have a normal distribution (for large *M*) with standard deviation .

This analysis is completely conventional and has been popular in astronomy for decades. However, researchers working on time series in fields such as biology, geology, medicine and engineering often use different methods, e.g. Tukey–Blackmann, auto-correlation function (ACF), partial ACF, maximum entropy/Burg method, ARMA modelling, and so on. To a large extent, these are all just slightly different ways of calculating and estimating the same underlying function, the power spectrum. But the Bartlett/periodogram approach has a number of features that help maintain its popularity with astronomers. These include: the sensitivity of the periodogram to nearly sinusoidal variations, as might be expected from the rotation of stars or the orbits of binary stars, the ease with which the DFT can be computed (via fast Fourier transform (FFT), algorithms), and its statistical properties. The Bartlett method yields power spectrum estimates at each frequency that are (asymptotically, for large *M* and *N*) independent of each other and normally distributed, which permits the use of simple least squares fitting of parametric models (for maximum likelihood estimation of parameters, and goodness-of-fit tests).

### (c) Tracking the evolution of the power spectrum

The count rate of the source during the observation shown in figure 3 was approximately 680 counts per second, meaning that we have approximately 1 X-ray count per characteristic time scale of the kHz QPO and that the QPO is clearly broad albeit with a high *Q*∼50. This means that we cannot easily examine the kHz QPO in the time domain: there are too few counts per cycle to recover the waveform of the oscillation, and if we average (‘fold’) the light curve on the period corresponding to the peak frequency in the power spectrum—in order to make up for the lack of counts per cycle—the result carries very little meaning because the QPO is not a coherent oscillation with a single, strictly repeating period. Future X-ray missions with larger collecting area may be capable of resolving these approximately ms variations in the time domain [25].

The width (incoherence) of the QPO is partly explained when we examine the dynamic power spectrum, as shown in figure 3*b*. This shows the variability power as a function of frequency and time, and clearly reveals the QPO frequency to be variable on short time scales. But again the data impose limits on what we can observe. It is necessary to average over enough data (e.g. approximately 50×1 s periodograms) to even detect the QPO; this sets a limit on the minimum time scale on which we can observe variations in the QPO properties. Nevertheless, we can go some way to removing the effect of *frequency drift* on the width of the QPO by using the ‘shift-and-add’ method to align the short-time power spectrum estimates to a common frequency [26]. Figure 5 shows the power spectrum estimated by averaging the 1 s periodograms over the entire observation with and without the frequency–shift correction. In this case, the frequencies were estimated by fitting the short-time periodograms with a simple parametric model using maximum likelihood fitting as discussed in Barret & Vaughan [24]. The ‘shift-and-add’ method does indeed recover a narrower QPO, although not without bias.

In Barret & Vaughan [24], we produced a time series of QPO frequencies, and then examined the power spectrum of this by fitting models using the maximum likelihood. This is in essence a *hierarchical* or *multi-level* analysis: we model the power spectrum of frequency variations of a QPO, which in turn are estimates derived from modelling short-time power spectral estimates of the light curve, which in turn is based on an even higher time resolution event list. On the basis of the power spectrum of the QPO frequency variations, we can estimate broadening of the QPO due to frequency drift on any given time scale, and correct this to recover the unobserved ‘intrinsic’ width of the QPO. Using this method, the recovered *Q*-value is at least as high as ∼260. Such high coherence is a challenge for most models of the behaviour of matter close to the neutron star [27]. The power spectrum of frequency variations is limited at high frequencies by the integration time required to determine the QPO frequency, and at low frequencies by the length of the observations. (RXTE is in a low Earth orbit, which means observations are interrupted at least every 96 min.)

This particular example was used to illustrate some of the complexities that arise when analysing even quite strong timing signals. In this case, the source is bright and the variations of interest (the kHz QPO) are relatively strong, persistent and sampled for a large number of cycles (and therefore well resolved in Fourier space). And yet recovering its basic parameters remains a challenge. Issues of high-frequency timing like this have been a concern in X-ray astronomy, and to some extent radio astronomy (e.g. pulsar timing), for many years. High time-resolution instruments are also used in optical and IR astronomy (e.g. [28]), so these kinds of problems are certainly not restricted to X-ray observations.

## 4. Nonlinear time series

Most work on fast timing of persistently variable astronomical sources is based on power spectral methods and related second-order statistics (cross-spectrum, ACF, etc.). These contain information on the statistical moments up to second order only. If the variability processes are linear, Gaussian and stationary, then no information is lost, everything is contained in the second-order statistics. (The phases of the Fourier components are random, and so it is no real loss to throw them away when we square the Fourier transform to estimate the power spectrum.) This remains an implicit (and often untested) assumption, rather than an open area of investigation, which may in part be due to the fact that constructing physical models of the power spectra has proved to be a serious challenge. However, the assumption of linearity and Gaussianity cannot be strictly true since we observe high-amplitude fluctuations (e.g. r.m.s. ∼20%), but the luminosity fluctuations cannot extend below zero. Other statistics, which are sensitive to the non-Gaussian and nonlinear properties of data, may provide more discriminating tests of existing models or point the way towards new ideas.

One potentially powerful diagnostic is the bispectrum or bicoherence. Where the power spectrum is a second-order spectrum, which can be thought of in terms of the product of two Fourier transforms [*X*( *f*_{j}) *X**( *f*_{j})], the bispectrum is a third-order spectrum formed from the product of three Fourier transforms [*X*( *f*_{j}) *X*( *f*_{k}) *X**( *f*_{j}+*f*_{k})]. The computation and interpretation of the bispectrum (and higher order spectra in general) is treated by Brillinger & Rosenblatt [29] and Nikias & Petropulu [30], and has been applied in geology [31], plasma physics [32] and speech processing [33].

The bispectrum is a complex quantity defined for the pair of frequencies ( *f*_{j},*f*_{k}). The phase of the bispectrum is often called the *biphase* and its modulus-squared amplitude is called the *bicoherence* (once it has been suitably normalized). The bispectrum provides a measure of the nonlinear interaction of frequency components in the time series. If the variations at the frequencies *f*_{j},*f*_{k},*f*_{j}+*f*_{k} are independent of each other they will have independent, random phases, and the bispectrum will tend to zero. If the variations at frequencies *f*_{j} and *f*_{k} are multiplicatively coupled, they will contribute a ‘side-band’ component at *f*_{j}+*f*_{k}. The phases at the three frequencies will be correlated, and the bispectrum will not average to zero. These properties make the bispectrum sensitive to the non-Gaussian and nonlinear characteristics of time series.

Figure 6 shows the bicoherence of the black hole X-ray binary system Cygnus X-1, estimated using a long series of X-ray observations taken in 1996 with RXTE [34]. (The source shows strong, persistent and reasonably stationary aperiodic variability, with not strong QPOs.) The non-trivial structure in the bicoherence indicates frequency coupling, but interpreting the bicoherence is not straightforward. The bicoherences at each bifrequency ( *f*_{j},*f*_{k}) are not independently distributed, and the Poisson noise in these data causes spurious bicoherence at high frequencies (the Poisson noise fluctuations are uncorrelated with the source flux but are not independent, since the amplitude of Poisson fluctuations scales with ). At high frequencies, the Poisson noise dominates over the intrinsic variations in the source flux (which drop off rapidly with increasing frequency). Nevertheless, bicoherence is now starting to be used as a powerful test for more advanced models of variability from black holes and neutron stars (e.g. [35,36]).

Figure 7 illustrates a different approach to the same data, discussed by Uttley *et al*. [37] and Vaughan & Uttley [34]. During this particular observation, the shape of the power spectrum of Cygnus X-1, at least at frequencies approximately 10^{−3}−10^{2} Hz, appears to remain quite constant through more than 2 days of observation. However, its amplitude does not. The normalization but not the shape of the power spectrum seems to scale (on average) with the X-ray count rate: brighter intervals are also more variable intervals. This ‘r.m.s.–flux’ relation is now quite well established in accreting systems (e.g. [38,39] and references therein). This can be thought of as an example of non-stationarity, since the variability properties are changing with time, or as a ‘static’ nonlinear transformation of the series,^{5} since normalizing an interval by its mean flux stabilizes the variance (there is connection here to Box–Cox transformations). However, even after accounting for this effect there remain some residual nonlinearities that are unexplained.

## 5. Non-stationary time series

We have, so far, simply assumed that the time series we are interested in are stationary (or, strictly, that they are random realizations of a stationary stochastic process). This is a necessary assumption for some analysis methods; for example, to obtain a consistent power spectral estimate using the Bartlett method, we assume that periodograms calculated for different time intervals are drawn from the same underlying process and so may be averaged to produce a meaningful estimate. But one thing we do know about these systems is that often the variability is not stationary.

The kHz QPO shown in figure 3 is clearly not stationary because its frequency varies on time scales as short as s. However, it may be the case that this higher level of variability is itself stationary, i.e. the frequency variations can be considered as a stationary stochastic process (as with the r.m.s.–flux effect mentioned above). Even this is not true on longer time scales because these sources undergo periods of outbursts and quiescence (figure 1), and move through a variety of ‘states’ that display characteristically different variability (and other) properties. Figure 8 shows some X-ray light curves of a black hole X-ray binary, GRS 1915+105. This particular source seems to move through a few reasonably well-defined ‘states’ (with specific high-frequency timing and energy spectral properties), but the pattern and speed of the transitions are quite remarkable and remain largely unexplained [40].

This non-stationarity imposes another limit on what we can achieve with conventional methods. It may be the case that on short time scales the system is ‘locally’ stationary, and so we can measure, for example, power spectral properties and examine their time evolution on the frequency–time plane. There is now a large volume of work on the evolution of power spectral features in X-ray binaries (see [18] for a review). But, in some cases, it may not be meaningful to estimate the power spectrum if the system is non-stationary on the shortest time scales on which we are able to estimate a power spectrum.

## 6. Slow, persistent variability

There are two facts about the fast timing of sources like X-ray binaries that help keep the field lively. The first is that there is a rich phenomenology: a range of power spectral features (noise, low-frequency QPOs, kHz QPOs, pulsations) that evolve systematically with time, in addition to transient events (rapid bursts, accretion outbursts) that can be detected, identified, classified and compared with other observables in the hope of finding associations that might help explain their origins [41,18,19]. The second reason is more pragmatic: there are a lot of good data. There are thousands of observations from RXTE of X-ray binaries over a wide range of luminosities and states, covering frequencies up to kHz, often with high count rates, and there are now many good fast optical and IR time series, sometimes simultaneous with the X-ray ones. The high photon event rates (meaning a relatively small contribution to the power from Poissonian fluctuations), high time resolution and multiple observations means that, in many cases, we can compute very good power spectra (or related statistical summaries of the variability), in the sense of having small variance and low bias on the estimates. And where there is non-stationarity the data are often sufficient to resolve the varying power spectrum in the time domain. But not all astronomical time series are like this.

At the other end of the frequency spectrum, we see much slower variations from active galactic nuclei (AGNs). AGNs are the luminous cores of galaxies thought to be powered by accretion onto a supermassive black hole (*M*_{BH}∼10^{6}–10^{9} solar masses). Over the past few years, several similarities have been observed in the X-ray variability of nearby AGNs and X-ray binaries, supporting the idea of *black hole unification* [42]. An important part of this project is the reliable estimation of power spectra from AGNs, which are expected to vary like X-ray binaries but slower by a factor of approximately 10^{5−8}, commensurate with their much larger black hole masses [43,44]. The time scales of interest for AGN variability studies therefore range from approximately 100 s to years, which leads to data analysis problems that are usually not present for the fast X-ray binary data.

If we take the power spectrum of Cygnus X-1 (figure 7) and scale the frequencies by a factor to those expected for an AGN, we would expect to see the power spectrum as ‘red’ (i.e. increasing to lower frequencies such as *f*^{−α} with ) down to frequencies as low as approximately 10^{−6} Hz (time scales of weeks), and there are good reasons for expecting the red spectrum to extend much lower. The sources also tend to have much lower count rates (fluxes) and so the Poissonian fluctuations in the count rate (measurement error) are relatively stronger. Clearly, an observation lasting only approximately an hour, such as that shown in figure 3, will not be sufficient to obtain a good power spectrum estimate for a much fainter, slower AGN. On the time scales of most observations, we are looking at ‘long memory’ processes when we observe AGNs, with significant power outside the observed frequency bandpass.

In order to obtain a useful power spectrum estimate, we need low variance, from large *M*, equivalent to many repeated samples at each frequency, and low bias, from long segments, i.e. large *N*, and good time sampling *Δt*, covering the full range of frequencies that show an interesting (non-white) spectrum. But it is almost never practical to observe an AGN continuously for many, many years, at high time resolution (e.g. seconds–minutes), to meet these criteria. AGN time series come from ‘monitoring campaigns’ comprising regular short snapshot observations of a target extending over a long period, or continuous but much shorter ‘long look’ observations, or combinations of the two. The monitoring approach often suffers from the effects of *aliasing* of high-frequency power [45,46], while the shorter, intensive observations often suffer from *leakage* of lower frequency power [47,45]. In both cases, significant power from frequencies outside the observed frequency bandpass biases the results of simple spectral estimators, which can be understood in terms of the effect of the window function (the time-sampling pattern) on the Fourier transform of the time series. Figure 9 shows a rare exception to the rule that long time-scale monitoring data are sparsely sampled: it shows data from the *Kepler* mission that continuously monitors the brightness of thousands of sources in a single patch of sky. Even with such good data there is a severe problem of spectral leakage of power from low frequencies.

It is a fundamental fact of spectral analysis that the window function biases power spectrum estimates [49,50], but this bias should disappear in the limit of large *N* (and small *Δt*) for stationary processes. The expectation value of the periodogram *I*_{j} at Fourier frequency *f*_{j} is
6.1
where *S*( *f*) is the true spectral density function, *F*( *f*) is the Fejer kernel, and *f*_{Nyq}=1/(2*Δt*) is the Nyquist frequency. The above is valid when there is no significant power in the spectrum above *f*_{Nyq}—i.e. the variations are well resolved by sampling at discrete times *Δt*. If there is significant high-frequency power above *f*_{Nyq} then additional terms (power from high frequencies, called *aliases*) also bias the periodogram.

The convolution by the Fejer kernel distorts the periodogram away from the true spectrum, and results from the finite duration of the time series [49], ch. 6. The expected difference between the true spectrum and the periodogram is the bias on the periodogram: *b*_{I}( *f*_{j}) = *E*[*S*( *f*_{j}) − *I*_{j}]. The Fejer kernel depends on the time-series duration, *T*=*NΔt*, and has a main lobe of width *Δf*∼1/*T* plus oscillatory side-lobes that decay as ∼1/*f*^{2} either side. It is these side-lobes that cause the leakage—the transfer of power between distant Fourier frequencies.

The basic problem in AGN timing analysis is that we cannot obtain sufficiently long and well-sampled time series to give negligible bias. This also means that there are not enough data to average over *M* non-overlapping intervals to give normally distributed averages. The time-series literature contains many discussions of power spectral estimation for ‘long memory’ processes (e.g. [51–54]). Popular methods often involve some kind of *tapering* of the time series, a modification of the window function to reduce leakage, as in the multiple taper method ([54] and references therein), or *differencing*, i.e. examination of *Δ*(*t*)=*x*(*t*)−*x*(*t*−1), rather than *x*(*t*).

### (a) The forward problem

The approach being developed by astronomers to solve this problem is based on ‘forward fitting’. Rather than trying to deconvolve the effects of the window function from the power spectrum estimate (the *inverse problem*), we include the window function biases in model power spectra that can be compared with the data (without tapering, etc.). The approach could be said to be *parametric*, being based on comparison of the data with specific well-defined models, rather than the *non-parametric* approach of deconvolution, but this is reasonable as the goal of much of this work is to compare physical models with data. The ‘response’ method (the work of Uttley *et al.* [45] based on the earlier work of Done *et al.* [55] and Green *et al.* [56] relies heavily on simulations. The essence of the method is as follows:

— estimate the power spectrum of the data

*x*(*t*_{i}),*i*=1,2,…,*N*using, for example, the periodogram (with optional smoothing/averaging);— define a power spectrum model

*S*(*f*;*θ*) in terms of parameters*θ*;repeat

*N*_{sim}times:— generate a random time series based on

*S*(*f*;*θ*);— sample the simulated series in the same manner as the real data (i.e. extract a subset from times

*t*_{i});— compute the power spectral estimate of the sampled simulation;

— at each frequency compute the mean and variance–covariance of the

*N*_{sim}estimates;— compare the mean spectral estimates from simulations with those from the real data.

The simulation step needs to be handled carefully. Usually, the algorithm of Timmer & König [57]^{6} is used to form random data in Fourier space that are then inverse Fourier transformed into a random time series. It is usually necessary to produce simulations much longer, and with finer time sampling, than the real data, from which a subset of the points is extracted in order to simulate the process of sampling a limited, discrete time series from a continuous process. Furthermore, simulations produced in this manner are realizations of stationary, linear, Gaussian processes. Nonlinearities such as the ‘r.m.s.–flux’ effect discussed above can be included in the simulation if necessary [37].

There is still some discussion about how best to perform the statistical comparison of data and simulations—e.g. how to define a fit statistic to use when searching for the best fitting model and its parameters. (It is not trivial to apply, for example, the *Whittle* likelihood in the presence of leakage and/or aliasing biases.) But the main advantage of this approach is that it should account correctly for the effect of the finite and discrete time sampling on the distribution of the power spectral estimates at each frequency (aliasing and leakage), by ‘folding’ the model spectra though the same sampling process. Other observational effects such as measurement errors or nonlinear detector response can also be included in the simulation procedure. The main disadvantage is computing time—a large number of (potentially long) time-series simulations needs to be performed for each choice of model and parameter values *θ*.

This approach has proved to be quite powerful for extracting reliable information about power spectral models from time series of (relatively) slowly varying AGNs that are (relatively) poorly sampled. But work remains to be done. Can this simulation-based approach be replaced by a faster analytical equivalent? What statistic(s) should be used for data–model comparison? Can non-stationarity be taken into account? There is also a problem with long-term monitoring campaigns in that the data are usually not very evenly sampled in time, yet standard Fourier methods (and many time domain methods) are best understood for evenly sampled time series. Generalizations of some methods to unevenly sampled data, such as the Lomb–Scargle periodogram [61], have been developed for special cases, such as searching for sinusoidal signals in white noise. To my knowledge, the effect of uneven sampling on the recovery of red power spectra remains an open problem [62,63].

## 7. Multi-variate time series

Bivariate, or more generally multi-variate, time series offer a new dimension to the analysis and potential for model testing. In many cases, the only information we can obtain on astrophysical phenomena is contained in the time and wavelength distribution of the light we receive. We might have observations from two telescopes operating simultaneously in different wavebands (e.g. optical and X-rays), or from one detector that tags events according to their energy (wavelength) so that we can form time series in different ranges of photon energy (wavelength). Depending on the situation the analysis of such data falls under different names: often this might be called *spectral variability*, meaning the time variability of the photon energy/wavelength spectrum (not the power spectrum!). Another situation is when we receive light from multiple images of the same event and form a time series from each image, such as with gravitationally lensed quasars [64].

In some cases having data in both time and wavelength dimensions allows for the recovery of spatial information from objects so distant they appear only as points in our best images. For example, the periodic rotation of single stars or the revolution of binary stars about their centre of mass can produce periodic variations in the light received on Earth when there is anisotropic emission such as star spots or a structured accretion disc. The rotation of these systems causes periodic modulations in wavelengths of atomic features (the relativistic Doppler shift), and by resolving the wavelength (and hence radial velocity) changes as a function of the phase of the period it is possible to extract information about the spatial structure by the method of Doppler tomography [65]. This can be viewed as a problem in multi-variate time-series analysis where we observe simultaneous time series for each of many distinct Doppler velocities.

A related method is that of ‘reverberation mapping’ of AGNs [66]. AGNs often show strong atomic emissions lines, such as the hydrogen Balmer series lines, excited by photoionization and recombination from the centrally generated continuum emission. The lines are therefore expected to respond to variations in the luminosity of the ionizing continuum, and time delays between continuum and line variations relate to the light travel time (and hence distance) between continuum and line-emitting regions. More specifically, the goal is to recover the *impulse response*^{7} of the line (its response as a function of time to an instantaneous impulse) and relate this to the spatial distribution of the emitting gas. If *x*(*t*) is the driving continuum, then the line emission *y*(*t*) is given by
7.1
where *ψ*(*t*) is the response function. As with Doppler tomography, the time variations are used to probe unresolved spatial structure. But AGNs do not vary periodically like rotating or binary stars and so the delays are seen as ‘echoes’ of the variations of a (red) noise process. Figure 10 shows an example of data obtained from a sustained reverberation mapping campaign. The recovery of the time response of the emission line is hampered by the irregular time sampling and the very red (low-frequency dominated) nature of the continuum variations, and the fact that the line response is likely to be nonlinear and non-stationary at some level.

The time delays between correlated time series of different but related origin tell us about the causal connections within the system. In reverberation mapping observations, we see the continuum increase in flux, and the line responds in a correlated manner with some (reproducible) delay. In some cases, we do not understand the physical relationship between different parts of the wavelength spectrum of sources and so correlation, time delays and multi-variate time-series methods in general are used to uncover these causal links.

### (a) Cross-correlation and cross-spectrum methods

As with most time-series problems one can consider the problem in the time domain or the frequency (Fourier) domain. The primary analysis tool for comparing two time series in the time domain is the cross-correlation or cross-covariance function [49,50]. If we have two series *x*_{t} and *y*_{t}, the cross-covariance is defined as
7.2
and the cross-correlation is . The *sample* cross-covariance is usually simple to calculate for long, evenly sampled, simultaneous time series. (*E* is the expectation.)

But when the two time series are not evenly sampled, perhaps not even sampled simultaneously, and are contaminated by measurement errors, it is less obvious how to perform the calculation. The solutions used by astronomers typically involve some kind of interpolation of one or both time series to evenly spaced and/or simultaneous times, or averaging together pairs of points with similar time differences [68]. Even when the issue of computation is resolved there remain questions for interpreting the cross-correlation: when is a peak an indication of a meaningful correlation between *x*_{t} and *y*_{t} if these are non-white noise processes? how should we estimate the time delay and its confidence intervals? or, more generally, how do we recover the response function? Here, as with power spectrum estimation, simulation-based Monte Carlo methods are proving to be useful [69].

The cross-spectrum is a related tool used for comparing the properties of two simultaneous time series in the frequency domain [49,50]. It is the Fourier transform of the cross-covariance,^{8} or the product of the Fourier transforms of the two light curves: *C*( *f*)=*X**( *f*)*Y* ( *f*). In principle, at least, the cross-spectrum contains the same information as the cross-covariance. But the cross-spectrum can often be computed quickly (using a FFT), and can be a more direct route to the response function relating the two time series. The cross-spectrum is complex valued, represented in the complex plane by an amplitude and a phase. If we write the complex Fourier transforms of the two light curves as *X*( *f*)=|*X*( *f*)|e^{iϕx( f)} and *Y* ( *f*)=|*Y* ( *f*)|e^{iϕy( f)}, then the cross-spectrum can be written as
7.3

In practice, the estimation of the cross-spectrum from a finite sample of real data usually requires a more involved procedure (e.g. Vaughan & Nowak [70]). Usually, the cross-spectrum is divided into its modulus and phase. The phase component, represents the phase difference between the two light curves (which corresponds to a time delay *τ*( *f*)=*Δϕ*( *f*)/2*πf* at that frequency). The squared magnitude of the cross-spectrum |*C*( *f*)|^{2} is used to define the *coherence* of the two light curves (different from the *coherence* of a peak in a power spectrum). The cross-spectral phase and coherence have been used to great effect on the evenly sampled X-ray data from sources like Cygnus X-1 [71,72], revealing complex, frequency-dependent time delays between X-rays of different energies. Recently, such methods have provided evidence for the causal connections between different parts of the accretion flow in X-ray binaries and in AGNs [73,74].

In order to illustrate the cross-spectrum, consider two time series *x*(*t*) and *y*(*t*), evenly and simultaneously sampled, which are related by a response function (equation (??eq7.1)). From the convolution theorem, we have
7.4

where *Ψ*( *f*) is the Fourier transform of *ψ*(*t*). Therefore,
7.5

where *S*_{x}( *f*)=|*X*_{f}|^{2} is the power spectrum of *x*(*t*) (neglecting a normalizing factor). Estimates of the cross-spectrum *C*_{xy} and power spectrum *S*_{x} therefore lead to a simple estimate of the response *Ψ*( *f*). If the two light curves are related by a single, linear response function (e.g. a delay or a smoothing), then they will have unity coherence (in expectation and in the absence of bias). Nonlinear relations and partially independent variations result in lower coherence.

Given *L* evenly sampled time series (corresponding to channels *l*=1,2,…,*L*) one could compute *L*(*L*−1)/2 cross-spectra, but in most cases there will be a great deal of redundancy between these. The phase delay (at a given frequency) between channels 1 and 3 cannot be independent of the phase delays between channels 1 and 2 and between channels 2 and 3. There will usually be fewer response functions relating the *L* bands than pairs of time series to form cross-spectra [75]. Another approach being explored by, for example, Uttley *et al.* [73] is to form only *L* cross-spectra by comparing the time series for each channel *l* against the summed time series of all the other channels. This has some utility when the channels are adjacent in wavelength and the number of underlying varying components to the emission is small compared with *L*.

Currently, the frequency domain approach (cross-spectrum) is applied to data that are evenly sampled, such as those from short, intensive X-ray observations, whereas the time domain approach (cross-covariance) is often applied to unevenly sampled data such as those from long-term AGN monitoring. But these tools may not be optimal for those tasks. Some outstanding problems in this area are: how best to compute reliable time delays (and, more generally, response functions) from unevenly and non-simultaneously sampled light curves? how best to recover transfer functions from very red time series (where again spectral leakage can bias the results)? how best to model multi-channel time series in terms of their simultaneous power spectra and cross-spectra (or time domain equivalents)?

### (b) Time domain methods

It seems to be the case that a great deal of astronomical time-series analysis is based around spectral methods, rather than time domain methods. This choice is largely motivated by scientific and practical considerations. There may also be historical and sociological factors involved. But I certainly do not wish to give the impression that only spectral methods are used in astronomy. For example, Scargle [76,77] discussed in great detail the time domain modelling of random time series in astronomy. Kelly *et al.* [78] and Miller *et al.* [79] developed tools for recovering the power spectra of AGNs by fitting to the data in the time domain. Press & Rybicki [80] and Zu *et al.* [81] developed time domain methods for reconstructing missing data from irregularly sampled AGN time series in order to recover time delays from pairs of unevenly sampled time series (the method is similar to the Kalman filter). Another method often used on unevenly sampled data is the *structure function*, although Emmanoulopoulos *et al*. [82] give strong arguments against its value as an analysis tool.

## 8. Concluding remarks

All the data presented in this paper have been obtained from public archives. The RXTE archive contains many thousands of observations of rapidly varying X-ray binaries and slowly varying AGNs; the *Swift* archive contains data on gamma-ray, X-ray and ultraviolet variability of GRBs but also a host of other transient and persistent sources; and time-series data from other wavebands for thousands of sources are being made available right now (e.g. the public *Kepler* archive). I encourage readers to download some data and explore the huge range of data and physics available to astronomical time-series analysis. It may be that we are just beginning to explore the science possible with astronomical time-series data [83].

## Acknowledgements

I would like to thank the meeting organizers (Nick Jones and Tom Maccarone) for allowing me to submit the manuscript rather late. I also extend thinks to Phil Uttley and the two referees for their valuable comments on earlier versions of this manuscript. This research has made use of NASA’s Astrophysics Data System.

## Footnotes

One contribution of 17 to a Discussion Meeting Issue ‘Signal processing and inference for the physical sciences’.

↵1 This is sometimes called the ‘Wolfer sunspot number’ in time-series textbooks. See Izenman [1].

↵2 Data files for all the examples in Feigelson’s review are available from http://xweb.nrl.navy.mil/timeseries/timeseries.html.

↵3 For example, events that fall within particular energy channels and within a particular part of the detector image. Not all events recorded by the detector will be caused by X-rays from the target source. Some event selection is required to reduce or remove background events, non-X-ray events and other contaminants.

↵4 Astronomers typically define QPOs in terms of the width of the peak,

*δf*, relative to its centroid frequency,*f*_{0}:*Q*=*f*_{0}/*δf*the*quality factor*(also called*coherence*) of the oscillation. Periodic, quasi-periodic and noise processes are loosely defined in terms of their*Q*values: unresolved peaks with very high*Q*are strictly periodic, peaks with finite*Q*down to*Q*=2 are defined as QPOs (quasi-periodic) and even broader features (*Q*<2) are usually called noise (aperiodic).↵5 Strictly, we should consider non-stationary and nonlinearity to be properties of models or processes, not data, which are but realizations of the process. In the present case, the model may be described as either nonlinear or non-stationary.

↵6 This algorithm seems to have been discovered several times. Very similar methods are discussed by Ripley [58], § 4.5, Davies & Harte [59] and Davis

*et al*. [60].↵7 In some astronomy contexts, this is also called as the

*transfer function*.↵8 The cross-spectrum and cross-covariance are Fourier counterparts, just as the power spectrum and auto-covariance are Fourier counterparts.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.