Improved timefrequency analysis of extrememassratio inspiral signals in mock LISA data
Abstract
The planned Laser Interferometer Space Antenna (LISA) is expected to detect gravitational wave signals from extrememassratio inspirals (EMRIs) of stellarmass compact objects into massive black holes. The long duration and large parameter space of EMRI signals makes data analysis for these signals a challenging problem. One approach to EMRI data analysis is to use timefrequency methods. This consists of two steps: (i) searching for tracks from EMRI sources in a timefrequency spectrogram, and (ii) extracting parameter estimates from the tracks. In this paper we discuss the results of applying these techniques to the latest round of the Mock LISA Data Challenge, Round 1B. This analysis included three new techniques not used in previous analyses: (i) a new Chirpbased Algorithm for Track Search for track detection; (ii) estimation of the inclination of the source to the line of sight; (iii) a MetropolisHastings Monte Carlo over the parameter space in order to find the best fit to the tracks.
I Introduction
The extrememassratio inspirals (EMRIs) of compact objects (white dwarfs, neutron stars or stellarmass black holes) with mass into massive black holes with mass can serve as excellent probes of strongfield general relativity (see Pau () and references therein for details). The detection of gravitational waves from EMRIs is one of the main targets for the planned gravitationalwave detector LISA (Laser Interferometer Space Antenna) LISA (). Astrophysical uncertainties mean that rate predictions are somewhat uncertain, but tens to thousands of EMRIs may be observed over the duration of the LISA mission Pau ().
Detection of EMRIs is not an easy task, however. The expected gravitational wave (GW) signals from EMRIs will be buried in instrumental noise and a foreground of galactic whitedwarf binaries. The instantaneous signal amplitude will typically be an order of magnitude smaller than the amplitude of noise and foreground fluctuations in the detector at the same frequency. Matched filtering allows the identification of weak sources buried in noise, but this relies on the generation of templates of possible signals present in the data. The EMRI parameter space is fourteen dimensional. Following the notation of Barack and Cutler AK (), these parameters are: the central black hole mass, , the mass of the inspiraling object, , the dimensionless spin parameter of the central black hole, , the orbital frequency, , and orbital eccentricity, , specified at some moment of time, , the inclination of the orbit to the equatorial plane of the central black hole, , the ecliptic latitude, , and longitude, , of the source on the sky, the orientation of the spin axis of the central black hole, and , the distance to the source, , and three phase angles, , and , which specify, respectively, the azimuthal orbital phase, periapsis precession phase, and phase of the orbitalplane precession at . The large dimensionality of the EMRI parameter space, combined with the long observation timescale ( year), makes fully coherent matched filtering using a template bank impractical as a detection method gairetal ().
Timefrequency searches are one possible alternative to coherent matched filtering. These have a much lower computational cost, albeit at the price of lower sensitivity. Timefrequency methods consist of building a spectrogram of the signal by dividing it into shorter segments and performing a Fourier transform on each of these segments, identifying possible tracks in the spectrogram, and then finally estimating source parameters from the tracks. We previously used timefrequency searches to analyse the Mock LISA Data Challenge (MLDC) Arnaud1 () Round 2 data set proc13 (), and found that we were able to get good parameter recovery. That previous analysis used the Hierarchical Algorithm for Clusters and Ridges (HACR) hacr07 () to identify tracks by clustering bright pixels in binned spectrograms, followed by a simple parameter estimation routine that used the track shape to infer the six intrinsic EMRI parameters and used the power variation along the tracks to estimate the two sky position angles proc13 ().
The HACR algorithm looks for clusters of any shape, and uses a binning technique that is based on rectangular pixel bins hacr07 (); wengair05 (); gairwen05 (). This is a rather inefficient way to search for EMRI tracks, which have a characteristic chirping shape. The most recent release of the MLDC, Round 1B, was a repeat of the type of EMRI data previously released in Round 2 Arnaud2 (). For this analysis we developed and used an improved track detection algorithm, which is tuned to look for EMRI type signals. We call this algorithm the Chirpbased Algorithm for Track Search (CATS; see CATS () for additional details). For the Round 1B analysis, we also applied parameterestimation methodology that had been improved relative to the techniques used in Round 2 proc13 (). The improvements included an estimation of the inclination of the source to the line of sight and a MetropolisHastings Monte Carlo search to find the parameters that best reproduced the detected tracks. The MLDC Round 1B data set consisted of five separate EMRI challenge data sets, 1B.3.11B.3.5; each challenge set comprised a single EMRI signal embedded in noise, with signaltonoise ratio (SNR) between 40 and 110 Arnaud2 (). We were able to detect EMRI signals in all of these data sets, and estimated 9 parameters for each challenge data set: the six intrinsic parameters — and , the two sky direction angles, and , and the angle between the SMBH spin vector and the line of sight to the source. We also independently estimated the plunge time, but this is determined by the other intrinsic parameters and was used to improve the determination of those parameters. Overall, the parameter recovery from the timefrequency techniques was very good and better than achieved in the Round 2 analysis proc13 ().
The paper is organised as follows. In Section II we describe the new track detection techniques and in Section III we describe the improvements to our parameter estimation algorithms. In Section IV we report the results of our analysis of the Round 1B data and compare these to the true parameters of the challenge data sets. Section V includes a summary and a discussion of possible future improvements to the search pipeline.
Ii Track detection
Our track search for the MLDC Round 2 analysis used HACR hacr07 (); proc13 (). We made a first cut at the Round 1B analysis using the same technique. We were able to detect several bright tracks for sources 1B.3.1, 1B.3.2 and 1B.3.3, and a cluster for 1B.3.4. The 1B.3.4 detection was very noisy, which was disappointing since this was a high SNR source. This failure prompted us to look at improved techniques more tuned to EMRI detection.
ii.1 Chirp Track Search
HACR does not use all of the available information about the signal. Specifically, we expect the tracks in the timefrequency spectrogram to be chirping curves characterized by three parameters: the frequency and its first two time derivatives and . (The third time derivatives of frequency along the tracks are vanishingly small for loweccentricity signals like those in challenges 1B.3.11B.3.3, and can be neglected for challenges 1B.3.4 and 1B.3.5 if sufficiently short sections of the tracks are considered.) Moreover, the first and second time derivative of frequency must be positive. This additional information can help to detect tracks. In fact, when we scanned the timefrequency spectrogram by eye, we found some tracks that HACR didn’t report, because we intuitively apply this additional information about the expected track shape when performing a visual search for tracks.
To make use of this information, we built a Chirpbased Algorithm for Track Search (CATS; see CATS () for additional details) to aid in track detection. This algorithm works as follows:

Construct and Time Delay Interferometry (TDI) channels for the challenge data sets.

For each channel, construct a spectrogram by dividing the data into time segments of equal duration, then Fourier transforming the data in each time segment after multiplying it by a Hanning window to reduce edge effects. The time segment durations were seconds for challenge 1B.3.1, seconds for challenges 1B.3.2 and 1B.3.3, and seconds for challenges 1B.3.4 and 1B.3.5.

Whiten the two spectrograms by dividing the signal power in each pixel by the expected LISA noise power spectral density, then construct a joint spectrogram by adding the two normalized spectrograms.

Define the starting and ending time of a chirp search. Construct a threedimensional grid of parameters in the space. For each point in the parameter space, build a potential track. Then add up the power in all pixels along this track between the starting and ending time of the search to obtain the total track power. (This is, effectively, a variation on a templatebased matchedfiltering search in the space.)

Find the brightest track and claim it as a detection. Set the power in all pixels along this track in the spectrogram to large negative values in order to keep future detected tracks from intersecting the current track. Repeat this step until all tracks are detected.
For this round, we did not attempt to search the whole spectrogram using this technique, nor did we search over the entire parameter space of . Instead we used the HACR algorithm to identify clusters in the spectrogram arising from the presence of signals in the data stream, and the approximate shapes of these tracks. We then targeted those areas for followup with the CATS algorithm. At the time of the Round 1B analysis, we had not robustly determined the correct thresholds for track detection. In the future, the track power threshold will be set based on a desired value of the false alarm probability. For this round, we set thresholds by hand, but performed a number of “sanity checks” to decide whether a given track detection was to be believed. These checks included:

Compare the detected tracks with a visual observation of the spectrogram.

Check if the tracks have sidetracks due to FFT bleeding, i.e., tracks that are parallel to them but one pixel off on each side.

Check if there is a clear harmonic structure typical to EMRIs, with sidebands caused by orbitalplane precession, i.e., several tracks that are nearly parallel (same ) with a fixed offset from each other.

Check if the track parameters are reasonable: e.g., if presumed sections of a track are found at early and late times, the frequency and frequency derivative at the later time should be greater than those at the earlier time.
In all the Round 1B challenge data sets, CATS successfully found all tracks that were detected by HACR or could be made out by eye. For some of the challenges, CATS found additional sidebands and CATS was also able to find several tracks in the Round 1B.3.5 data set, in which HACR was unable to identify any interesting clusters. This was particularly pleasing, since previous work wengair05 (); gairwen05 (); hacr07 () had indicated that sources of low SNR, like 1B.3.5, would be out of reach of timefrequency methods.
We estimated the systematic uncertainties of CATS by comparing the results of track searches using different time steps and different gridding of the chirp parameter space. These comparisons suggest uncertainties on the order of Hz, Hz/s, and Hz/s. Of course, the inherent errors in chirp parameter estimation of lowSNR signals, particularly for short tracks, could be much greater than these systematic uncertainties.
ii.2 Radon Transform Search
We have also been developing and testing another timefrequency technique that utilizes the modified Radon transform. For this search, we sum power in pixels along expected tracks in the timefrequency spectrogram of a given source. The expected tracks are calculated by integrating the known evolution equations for EMRIs for a given set of initial parameters. Tracks of the harmonics and sidebands of orbital frequencies are included in the summation, and we search through the whole parameter space for the maximum signaltonoise ratio. This technique was not fully developed at the time of submission, so although the results were compared to the results obtained from the CATS algorithm, only the latter results were submitted for this round. In the future, further improvement of the Radon transform algorithm will be to needed in order to overcome problems of getting stuck in local maxima of the parameter space. We are also investigating similar techniques using the Hough transform which sums over number counts (instead of powers) of bright pixels along the track. A search using the Radon or Hough transform normally takes minutes to complete on a laptop when a good initial guess has been made, or hours starting from a poor guess.
Iii Parameter Extraction
Parameter extraction based on the tracks identified by CATS consisted of two steps. In the first step we made rough estimates of the parameters. This was the approach used for parameter extraction in our analysis of the Round 2 data sets proc13 (); expowsymp (). An EMRI waveform is characterized by three fundamental frequencies — the orbital frequency, , the perihelion precession frequency, and the frequency of precession of the orbital plane, . Each track identified in the spectrogram is a harmonic of these three frequencies characterized by three integers, , with . The analytic kludge (AK) waveforms of Barack and Cutler AK (), which is the model used to generate the MLDC data sets, are quadrupole in nature and precession effects are included in an ad hoc way by precessing the observer about the source. This restricts and, in practice, harmonics with are suppressed relative to the harmonics and are not detected by a timefrequency analysis.
If three tracks are detected at a given time — two harmonics plus a sideband for one of those harmonics, then the separation of the harmonics is a measure of , the separation of the sideband is a measure of and, if the harmonic numbers of one of the harmonics can be guessed, the frequency of that harmonic then gives . For a specification of and , determines , determines and is a parameter that has been measured directly. The rate of change of frequency of one harmonic, , then determines . If the plunge time and the frequency of at least one harmonic at some other time can also be estimated, it is possible to iterate on and to determine all six intrinsic parameters for the source.
The motion of the LISA detector around the sun leads to a modulation in the detector response that depends on sky position. From the power variation along a track it is therefore possible to estimate the sky position of the source, by using the lowfrequency approximation to the LISA response cutler98 (); proc13 (). For the Round 1B analysis, we also tried to use the relative power in the sidebands to estimate the orientation of the source relative to our line of sight. Physically, the relative sideband power can only depend on the inclination of the source to the line of sight, , which is given in terms of the standard waveform parameters by . The relative harmonic power depends only on the absolute value of , and so we can choose . It is not possible to determine the two source orientation angles, , from the single number and so for the MLDC entry we returned as an additional parameter.
It turns out that, for a fixed value of , the relative power in different harmonics depends only on and the inclination of the orbit, (more details will be given in CATS ()). The power ratio in sidebands can thus be used to estimate and also to improve the determination of . This is illustrated in Figure 1. The curves in this figure show the power in the various sidebands relative to the power in the harmonic as a function of , for a fixed value of . The horizontal dashed line shows the power ratio estimated from the spectrogram, from which we infer an estimate . In the analysis of the Round 1B data, we found for all the data sets (specifically for Challenges 1B.3.1–1B.3.5 respectively). This suggested the determination of was weak. The true values were respectively, which implies the measurement was quite accurate, except for 1B.3.1 in which the sideband detected was very weak, and so the harmonic power ratio was noise dominated. The value seems to be typical for points drawn from the MLDC priors, and the tf spectrogram has limited ability to pin down more accurately. However, the power ratio did help to improve the estimation of for these data sets.
After determining rough parameter estimates via the above methods, the second stage of parameter estimation, which was used for the first time in this round of the MLDC, was to refine the estimates using an implementation of the MetropolisHastings Monte Carlo (MHMC) algorithm. The MHMC algorithm is described in many standard textbooks (e.g., gregory ()) and is being used extensively for the detection of various sources in the MLDC Arnaud2 (). An MHMC search constructs a sequence of points in parameter space, . From a point , a value for the next point in the chain is proposed by drawing from a proposal distribution . A random parameter is then drawn from a uniform distribution and if the move is accepted and the chain moves to , otherwise the move is rejected and . Here is the MetropolisHastings ratio
(1) 
is the prior on the parameters and is the likelihood of the data given the parameters . For the timefrequency MHMC search, we used priors that were uniform within the parameter ranges specified for the MLDC, and based the loglikelihood on the overlap between the observed whitened timefrequency spectrogram and a theoretical spectrogram built from the AK model (see CATS () for additional details). Our proposal distribution was constructed from the Fisher Information Matrix for our waveform model.
For the Round 1B analysis, the MHMC chains were seeded at the rough parameters estimated in the first stage and were allowed to run for points. The search used plunge eccentricity and plunge time as parameters rather than initial frequency and eccentricity, with a prior on the plunge time estimated from the detected timefrequency tracks. The MHMC chains refined the initial parameter estimates, but did not move far from the initial parameters, as is clear by comparing the “Rough” and “Final” parameters in Table 1. The chains moved somewhat further for 1B.3.5, for which a misidentification of the harmonics led to a poor initial parameter estimation. We also ran chains that started at random points in the parameter space, and these converged to the same parameter estimates, which gave us confidence that the code was working as expected.
Iv Results
In Table 1 we tabulate the true and recovered parameters for our search of the MLDC Round 1B data. In the table “Rough” denotes the first, rough parameter estimates, “Final” denotes the refined estimates obtained from the MHMC search and “True” denotes the true parameters. We see that, overall, we obtained good estimates of all the parameters, except for the lowest SNR challenge, 1B.3.5. The “Rough” parameter estimation was already quite good in most cases, but the parameter estimation was further improved by the MHMC, except in the case of challenge 1B.3.1 where the final parameters were worse than the rough parameters. Our sky position estimator has a degeneracy between antipodeal sky positions and , which reflects the real degeneracy between these points in the low frequency approximation cutler98 (). We used the convention that we returned values with , which meant that we returned the antipodeal position for sources 1B.3.4 and 1B.3.5. The sky position estimates were otherwise very good. We make some source specific comments below.
Challenge  S  (mHz)  

Rough  0.491  4.9  0.698  10.5  9.56  0.1920  0.2102  0.424  
1B.3.1  Final  0.4941  4.939  0.6667  10.404  9.787  0.1921  0.1886  0.1938 
True  0.5526  4.9104  0.6982  10.2961  9.5180  0.1920  0.2144  0.4395  
Rough  0.393  4.59  0.641  9.593  5.118  0.3422  0.2251  1.48  
1B.3.2  Final  0.4028  4.656  0.6371  9.817  5.250  0.3423  0.2017  1.423 
True  0.3597  4.6826  0.6380  9.7711  5.2156  0.3423  0.2079  1.4358  
Rough  1.100  0.628  0.525  9.713  5.2285  0.3426  0.1958  0.921  
1B.3.3  Final  1.013  0.7348  0.5318  9.756  5.250  0.3426  0.1936  0.9091 
True  0.9817  0.7097  0.5333  9.6973  5.2197  0.3426  0.1993  0.9282  
Rough  1.508  1.19  0.62  10.04  9.51  0.8473  0.4598  1.663  
1B.3.4  Final  0.7805  4.168  0.6133  9.989  9.518  0.8512  0.4564  1.658 
True  0.9802  0.9787  0.6251  10.1047  9.558  0.8514  0.4506  1.6707  
Rough  1.493  0.565  0.610  9.925  9.504  0.8823  0.4435  1.39  
1B.3.5  Final  0.726  4.45  0.618  10.49  9.500  0.887  0.417  1.47 
True  1.1092  1.0876  0.6583  9.7905  1.0334  0.8322  0.4269  2.3196 
iv.1 Challenge 1B.3.1
The high SMBH mass meant that there were a number of neardegeneracies in the parameter space for this challenge, which made individual parameters difficult to resolve despite the high SNR. This problem was exacerbated by the fact that the inclination was such that only one bright sideband was seen for each harmonic. Although other, weak, sidebands were detected, the MHMC moved towards lower eccentricity and inclination since this put more of the signal power into the brightest harmonic and hence “improved” the fit. In a future implementation of the search, we plan to include a penalty associated with parameter space points that do not reproduce all the tracks that have been seen. The final parameters of the MHMC search were such that the strongest detected track for the harmonic should be the sideband, and the secondary track should be the sideband. However, the values of and were such that the sideband should have been brighter than the sideband. Since we could not detect another sideband in the presumed location of , we manually changed the inclination angle , which decouples from the other parameters, to get a fit to the detected tracks. The resulting estimates for and were quite poor as a result. This was the only challenge in which the rough parameter estimates were better than the final parameter estimates. In future searches, we may use only the rough parameter estimation for high mass SMBH sources as it is then easier to control degeneracies.
iv.2 Challenge 1B.3.2
The MHMC search pushed the SMBH mass all the way to the maximum value allowed by the prior, which suggested a lack of accuracy. The initial rough parameter estimate matched the data well, although not quite as well as the final MHMC result. We submitted both results, although this was artificial in a sense, because we felt the final point was a better fit and only thought it was wrong because we had artificially tight prior constraints on . When the prior was relaxed, the best fit was found to be for a slightly higher central black hole mass, around .
iv.3 Challenge 1B.3.3
Again, the MHMC search pushed the SMBH mass to the maximum value allowed by the prior. This was, however, very close to the value detected by the initial rough estimate. Once again, if a less restricted prior was used, the best fit would have been for a slightly higher central black hole mass. For both 1B.3.2 and 1B.3.3, the true black hole mass was quite close to the edge of the prior which is why this issue arose.
iv.4 Challenge 1B.3.4
Although sections of multiple harmonics and sidebands were detected, the duration of each detected segment was quite short. This suggested that parameter estimation might not be as accurate as usual. However, the final parameters reproduced the detected tracks well, which gave us confidence in the results, and it is clear from Table 1 that parameter recovery was very good.
iv.5 Challenge 1B.3.5
The tracks were difficult to find because of the low SNR. Nevertheless, we succeeded in finding a track corresponding to the harmonic with a possible sideband near , and tracks corresponding to the and harmonics, each with one sideband, near s. The short duration of the tracks made accurate parameter estimation difficult, and the best fit parameters did not yield tracks that passed through the tracks detected near (although they were close). We expected these parameters to be the worst of the five Challenges, which is confirmed by Table 1. However, we were unable to detect anything in the corresponding data set in MLDC Round 2, so the performance of our approach in Round 1B is a significant improvement.
V Summary and future plans
The timefrequency analysis of the Round 1B data was very successful and our entry had the most complete parameter recovery out of the three submissions for this challenge. The other challenge entries used MHMC matched filtering. Ultimately that will be more sensitive and have much better parameter recovery than timefrequency analysis, but the timefrequency approach clearly has some parameter estimation power and will potentially be a useful first step in an analysis pipeline. There are a few ways in which the search could be improved. The CATS algorithm, as applied to Round 1B, was not automatic, but interesting regions of parameter space were selected either by eye or by using the output of a HACR search. In the future, we will assess the performance of the CATS algorithm without using prior information. We will also compare the Radon transform algorithm outlined earlier against the CATS algorithm. The parameter extraction cannot be significantly improved. However, the template spectrogram model used in Round 1B to evaluate the likelihood in the MHMC search is only an approximation. It needs to be properly compared against mock data and might be improved. We also want to investigate imposing penalties in the likelihood evaluation to ensure the parameters predict all the tracks that are detected in the spectrogram. It should also be possible to use the total power in the tracks to estimate the total SNR or equivalently the distance to a source, which we have not done so far. The parameter estimation that we have done to date has been very specific to the AK model used for the MLDC data. The general techniques will generalise to more accurate waveform models but the mapping between harmonic frequencies and physical parameters will change. We plan to derive analogous methods for parameter extraction from “numerical kludge” waveforms GG06 (); NK () as a first step towards this nontrivial generalisation.
So far, the MLDC EMRI data sets have each contained a single EMRI signal in instrumental noise. The prior that there is only one signal in the data set is a very powerful one, but also very unrealistic. The true LISA data stream will contain a few tens to as many as a thousand EMRI sources gairetal () plus signals from white dwarf binaries and supermassive black hole mergers. Timefrequency detection algorithms, plus the simple parameter extraction techniques that we have used to date, will have limited effectiveness in such situations. However, it is a useful exercise to see at what point these techniques no longer work, and it is likely that the methods will still be able to detect the brightest sources even in the presence of confusion, although track detection thresholds will have to be modified to account for the confusion “noise”. The MLDC Round 3 will begin to address these issues, as there will be only one EMRI data set containing five overlapping sources. We will explore confusion while attempting to analyse Round 3. However, the sources will also be low SNR, so it is not yet clear whether timefrequency techniques will be able to detect anything. If necessary, we will create similar data containing brighter sources. We also plan to explore the efficiency of these methods when applied to data containing a single EMRI with a reduced galaxy to mimic the effect of confusion from whitedwarf binaries.
Acknowledgements.
JG thanks the Royal Society for support and the Albert Einstein Institute for hospitality and support while part of this work was being completed. IM was partially supported by NASA ATP Grant NNX07AH22G to Northwestern University. LW’s work is supported by the Alexander von Humboldt Foundation’s Sofja Kovalevskaja Programme funded by the German Federal Ministry of Education and Research.References
 (1) AmaroSeoane P, Gair J R, Freitag M, Miller M C, Mandel I, Cutler C J and Babak S, Class. Quantum Grav. 24 R113 (2007).
 (2) Danzman K et al. LISA PrePhase A Report, Report MPQ 233 (1998).
 (3) Barack L and Cutler C, Phys. Rev. D 69 082005 (2004).
 (4) Gair J R, Barack L, Creighton T, Cutler C, Larson S L, Phinney E S & Vallisneri M, Class. Quant. Grav. 21, S1595 (2004).
 (5) Arnaud K A et al., AIP Conf. Proc., 873 619 (2006).
 (6) Gair J R, Mandel I and Wen L. Accepted to a special issue of the Journal of Physics Conference Proceedings, arXiv:0710.5250 (2007).
 (7) Gair J R and Jones G, Class. Quantum Grav. 27 1145 (2006).
 (8) Wen L & Gair J R, Class. Quantum Grav. 22 S445 (2005).
 (9) Gair J R & Wen L, Class. Quantum Grav. 22 S1359 (2005).
 (10) Arnaud K A et al., Class. Quantum Grav. 24 S551 (2007).
 (11) Gair J R & Mandel I, in preparation (2008).
 (12) Cutler C, Phys. Rev. D57 7089 (1998).
 (13) Wen L, Chen Y & Gair J R, AIP Conf. Proc., 873, 595 (2006).
 (14) Gregory P, Bayesian Logical Data Analysis for the Physical Sciences (Cambridge: University Press) (2005).
 (15) Gair J R & Glampedakis K, Phys. Rev. D73 064037 (2006).
 (16) Babak S, Fang H, Gair J R, Glampedakis K and Hughes S A, Phys. Rev. D75 024005 (2007).