Prospects for observing extreme-mass-ratio inspirals with LISA
One of the key astrophysical sources for the Laser Interferometer Space Antenna (LISA) are the inspirals of stellar-origin compact objects into massive black holes in the centres of galaxies. These extreme-mass-ratio inspirals (EMRIs) have great potential for astrophysics, cosmology and fundamental physics. In this paper we describe the likely numbers and properties of EMRI events that LISA will observe. We present the first results computed for the interferometer that was the new baseline mission submitted in January 2017 in response to the ESA L3 mission call. In addition, we attempt to quantify the astrophysical uncertainties in EMRI event rate estimates by considering a range of different models for the astrophysical population. We present both likely event rates and estimates for the precision with which the parameters of the observed sources could be measured. We finish by discussing the implications of these results for science using EMRIs.
The first gravitational wave (GW) detections by LIGO [1, 2] in 2015 began the new field of GW astronomy. These events indicated that the ground-based interferometers will see many more systems over the coming decades . However, seismic noise limits the sensitivity of ground-based detectors to frequencies above and hence to systems with total mass no more than a few hundred solar masses. Observational evidence indicates that the centres of most galaxies contain massive black holes (MBHs), with masses between a few tens of thousands and a few billion solar masses . Mergers involving MBHs are powerful sources of GWs, but these GWs can only be observed from space. ESA selected GW detection from space as the theme to be addressed by the L3 mission in the Cosmic Vision programme, scheduled for launch in 2034 . The technology for a space-based interferometer was successfully demonstrated in 2016 by the LISA Pathfinder satellite , which paved the way for an ESA call for L3 mission proposals that closed in early 2017. Given the success of LISA Pathfinder, it is highly likely that the successful proposal will be the Laser Interferometer Space Antenna (LISA) .
LISA will have sensitivity to GWs in the millihertz band, which are generated by merging systems having total mass in the range –. LISA is expected to observe a variety of sources, including hour-period compact binaries (primarily white dwarf–white dwarf binaries) in the Milky Way, the early inspiral phase of some of the heavier stellar-origin black hole binaries that LIGO has shown to exist , and possibly cosmological sources such as GWs from cosmic string cusps or a stochastic background generated by phase transitions in the early Universe. However, the primary source for LISA are the MBHs in the centres of galaxies. These MBHs will generate GWs either when they merge with other MBHs, which is expected to occur following mergers between their host galaxies, or when they merge with much smaller compact objects formed from stars. These latter systems are called extreme-mass-ratio inspirals (EMRIs), because of the large difference in mass between the two objects involved. Galactic MBHs are typically surrounded by a stellar cluster, and EMRIs occur when compact objects formed as the end point of the evolution of stars in that cluster are captured and gradually fall into the central MBH. EMRIs occur over long timescales, meaning that – cycles of gravitational radiation are generated while the compact object is in the strong gravitational field region close to the central MBH. In addition, the process by which EMRIs begin, through a scattering capture of the compact object by the MBH, tends to lead to high initial eccentricities, with the result that EMRI orbits are expected to retain significant eccentricity and non-equatorial inclinations when they are radiating in the LISA band. These properties lead to a great richness in EMRI gravitational waveforms, created by a superposition of the three fundamental frequencies — the orbital frequency, the perihelion precession frequency and the frequency of precession of the orbital plane. This complexity contains a great deal of information that can be used for science. EMRI GWs can be used both to measure the parameters of the system to high precision and to detect small deviations in the waveforms from the predictions of general relativity that are indicative of a breakdown of the theory. EMRIs thus have tremendous potential for astrophysics, cosmology and fundamental physics; this will be discussed again at the end of this article, but we refer the reader to [9, 10] for a more comprehensive discussion.
In this article we will compute the likely numbers of EMRI events and the precision with which LISA will be able to measure their parameters. There have been previous studies computing expected numbers of EMRI events . However, the astrophysical model employed in those calculations was a combination of simple power laws and no attempt was made to quantify the uncertainties in that model. EMRI parameter-estimation studies have also been carried out [12, 13], but only for a small sample of representative cases and not for a full astrophysical population. We address both of these shortcomings in the current article. We compute event rates for several different astrophysical models that cover the range that is plausible given current observational constraints and compute estimates of the parameter-estimation accuracies for all the events in each population. In addition, we compute these results for the first time for a LISA detector with six laser links, which was proposed as the new mission baseline in the response to the ESA call in January 2017 .
This article is a summary of a paper that has recently been submitted  and we refer the interested reader to that paper for more details of the calculations and final results. This paper is organised as follows. In Section 2, we describe the models we employ for the detector sensitivity, the EMRI waveforms and the astrophysical population of EMRI sources. In Section 3 we present a comparison of the likely numbers of EMRI events that will be observed under each model and the parameter-estimation precision. In Section 4 we describe the implications of these results for science with EMRIs, before finishing in Section 5 with a short summary of conclusions.
To estimate the detectability of EMRIs by LISA, we must model three different things — the intrinsic astrophysical population of EMRI events; the GWs generated by EMRI events; and the sensitivity of our detector to those GWs. We describe the models we use for each of these in this section.
2.1 Intrinsic astrophysical population
Our astrophysical population model includes models for the MBH mass and spin distribution, the intrinsic rate of EMRIs occurring in systems containing black holes with particular properties and the properties of the compact object in the EMRI.
MBH mass distribution Previous estimates  of EMRI rates have used a simple power law MBH mass function of the form
which is known to be a good fit to the observed MBH population in the relevant range . Values of and are usually assumed. A more sophisticated model that self-consistently evolves the MBH population is now available , and has been used to predict the expected number of MBH binary mergers that will be observed by LISA . We will denote this self-consistent MBH model as “B12”. This model generates mass functions that are roughly consistent with the above power-law form, with and . To obtain a lower bound on the EMRI rate we will also consider a strict power-law MBH mass function with and . We will denote this model as “G10”.
During a major merger, we expect the stellar cusp around the MBH to be disrupted. The cusp will regrow, but this takes a quarter of the relaxation time for a Milky Way-like galaxy , during which the standard EMRI picture does not apply. We include this effect in our model by assuming an MBH is unavailable as a host for EMRIs for a time
following a merger with total binary mass and mass ratio . The pre-factor depends on the assumed – relation, which is discussed below. We cannot apply this correction for the “G10” MBH population, as we do not track mergers in that case.
MBH spin distribution The “B12” model also tracks MBH spins. These are driven partially by mergers but primarily by accretion. We denote the MBH spin distribution computed in “B12” by “a98” as it predicts  that most MBHs will be spinning close to the maximal imposed limit, i.e., with . To understand the significance of this assumption we also consider a “flat” model with spins uniformly distributed in the range , and a conservative “a0” model in which all spins are set to .
EMRI rate per black hole We base our EMRI rate per MBH on 
If this rate was maintained over the age of the Universe, it would overgrow lower mass black holes, particularly when we account for the fact that for every EMRI there will be – “failed EMRIs” that directly plunge . Therefore, we add two corrections. First, we compute , where is the relaxation time in the MBH stellar cusp and is the time it would take to deplete the stellar cusp of compact objects through EMRIs and direct plunges. Second, we compute , where is the total mass that would be accreted by the MBH over its evolution history according to (3). The final EMRI rate is taken to be . The factor corrects for the time taken for the stellar loss cone to be repopulated by diffusion, while the factor imposes a limit (of one e-fold) on the total mass that an MBH can gain by compact object accretion.
Black hole – relation The factor in (2) and the stellar relaxation time depend on the velocity dispersion in the stellar cluster. This can be determined from the black hole mass via the – relation. We use the model of  as our reference model “Gultekin09”, which predicts , plus an optimistic case “GrahamScott13” , which predicts , and a pessimistic case “KormendyHo13” , which predicts .
Compact object properties EMRI rates are expected to be dominated by inspirals of black holes, as these are visible to much greater distances and are intrinsically enhanced by mass segregation [9, 11]. Typically the mass, , of the inspiralling object has been taken to be fixed at . Given recent LIGO results [1, 3] we will consider both this standard case and one in which all EMRIs are assumed to have . We take the eccentricity of the EMRI orbit at plunge to be uniformly distributed in the range . EMRIs on prograde orbits into more rapidly spinning MBHs can be more easily captured as inspirals rather than direct plunges. We account for this by enhancing the EMRI rate by a factor that is a function of the MBH spin , and by using the inclination distribution .
|ID||MF||MBH spin||Erosion||– relation||CO mass||EMRI rate |
Overall, we use different models that are constructed using different combinations of the above ingredients. The models are summarised in Table 1. M1 is our reference model. M10 and M11 will be used only to estimate event rates, as degeneracies present for non-spinning EMRIs make parameter-estimation calculations based on the Fisher matrix more challenging. Full parameter-estimation results will appear in the longer journal article based on this work.
2.2 EMRI gravitational waveforms
The extreme-mass-ratio in an EMRI system means that the emitted GWs can be computed accurately using black hole perturbation theory. This involves calculating the self-force acting on the inspiralling object. We refer the reader to  for details on this approach. Self-force calculations have not yet been completed at the necessary order for objects on arbitrary orbits about spinning black holes, and those calculations that do exist for simpler cases are computationally intensive and therefore not well suited to studies of populations such as this one. As an alternative, there are two different families of kludge waveforms available. These are approximate models that have been constructed to include the most important qualitative features of EMRI waveforms, while not necessarily being quantitatively accurate. The analytic kludge (AK) model  constructs an EMRI waveform starting from the simple GW emission from a Keplerian orbit, and then imposing precession and evolution of the orbit onto the system. The precessional and evolution effects are governed by low-order post-Newtonian equations, which are accurate in the weak-field, low-velocity regime. The AK model is cheap to generate, but it rapidly goes out of phase with more accurate waveforms. The numerical kludge (NK) model [27, 28] uses Kerr geodesics as the basis for constructing EMRI orbits, adding inspiral determined by evolution equations that are based on post-Newtonian expansions, but augmented by fits to perturbation theory results. The NK waveform is generated from the trajectory using a flat-space emission formula . NK waveforms are much more faithful to accurate inspirals than AK waveforms, staying in phase down to the last stages of inspiral ; however, they are more computationally expensive.
For the purposes of this study we need to be able to generate large numbers of gravitational waveforms to determine the detectability and parameter-estimation precision for astrophysical populations of EMRI events. We therefore use the computationally cheapest of the three models, the AK. While these waveforms are not faithful to the true EMRI signals produced in nature, it is believed that they capture the richness of real EMRI gravitational waveforms and therefore will provide an accurate guide to signal to noise ratios and parameter-estimation precisions. To attempt to characterise the uncertainty from the waveform choice, we consider two different versions of the AK. In the original AK model , the waveforms are terminated when the object reaches the Schwazrschild innermost stable orbit (ISO). This is inaccurate for prograde inspirals into spinning black holes, which can get much closer to the central black hole. As the spins of the black holes in our population tend to be quite high, this is likely to significantly underestimate the detectability of many events. A simple fix is to continue the inspiral until the EMRI reaches the Kerr ISO. As the post-Newtonian expressions on which the AK model is built break down in the strong-field regime at the end of the inspiral, this continuation is not going to be accurate and will most likely lead to an overestimate of detectability. Therefore, we expect the two sets of results to bracket the true answer.
2.3 LISA sensitivity to EMRIs
An EMRI will be detectable by LISA if it is sufficiently loud. This can be characterised by a requirement on the matched filtering signal-to-noise ratio (SNR) , which the source has in the detector data stream
This expression depends on the final detector configuration through , the power spectral density of noise in the detector. We use the LISA noise model described in 
where , , and . We set the detector arm length to Gm. The complexities of EMRI data analysis mean that a relatively high SNR will be required for confident detection. The value was historically assumed , but results from the Mock LISA Data Challenge suggest that EMRIs with SNR as low as could be identified , albeit under somewhat idealised conditions. We will present results for both thresholds.
We assume conservatively that the duration of the LISA mission is two years and that any EMRIs that are detected must plunge during the mission lifetime. EMRIs accumulate SNR gradually over the whole inspiral, so LISA will also observe systems that are in the final stages of inspiral but do not reach plunge over the mission duration. The event rates reported here are therefore conservative in this regard. We use two years as this was the default used in the suite of recent studies connected to possible LISA configurations, such as . The recent mission proposal  uses a longer nominal mission duration of four years. To extrapolate to longer mission durations we can multiply by the ratio of the mission duration to two years, so a four year mission will have approximately twice as many EMRI events as reported here. Again, this will be conservative, as some EMRI signals that do not accumulate enough SNR over two years will pass the SNR threshold when integrated over a longer observation.
To evaluate parameter-estimation precision we use the Fisher matrix, defined by
where denotes the waveform model parameters. An estimate of the precision of measurement is given by (no sum over ).
3.1 Event rates
Table 2 indicates the number of EMRI events that would be detected over the two year mission lifetime for each model, using SNRs computed for AK waveforms with the standard Schwarzschild plunge condition. For the reference model we would expect to observe several hundred events with SNR above , of which approximately one third would have SNR above . The number of observed EMRIs can be no more than a factor of larger (M12), if we assume that MBHs are always available as EMRI hosts and MBHs gain no mass from direct plunges of compact objects. The EMRI rate could be as low as , if we assume a mass function that falls steeply toward lower MBH masses and that MBHs have small spins (M11). That model is overly pessimistic, and M5 and M8 probably give more realistic lower bounds, which are a factor of smaller. The majority of the models predict several hundred EMRI events, so this figure appears to be fairly robust to astrophysical uncertainties.
If the inspiralling compact objects in EMRIs tend to be more massive (M4), we see a comparable number of events. While such EMRIs can be seen to larger redshift, this is compensated by a decrease in the intrinsic rate to prevent MBH overgrowth. Table 3 shows the corresponding results computed using the Kerr plunge criterion to terminate the EMRI. These rates are somewhat higher, as expected, but only by a factor –.
In our model we have prevented MBH overgrowth through compact object accretion by imposing a constraint on the EMRI rate. An alternative solution is that the number of light MBHs is significantly reduced, because these MBHs grow rapidly from EMRI consumption. Results in the table are divided into mass bins and we see that between and of the events have MBH mass greater than , so even if no lighter MBHs exist we should see approximately one hundred events in the reference model. The rate drops significantly, by a factor of or more, if we impose a cut off at , when using the Schwarzschild plunge condition. However, a higher proportion of the events are heavier mass when we use the Kerr plunge condition. This is because in that model the EMRI can get closer to the MBH, shifting the GW emission to higher frequencies and hence providing sensitivity to heavier MBHs. The prospects for significant numbers of EMRI detections are therefore good even if the number of lower mass MBHs is significantly depleted.
|Number of events in mass range|
|M1||20 (10)||240 (60)||110 (50)||10 (0)||380 (130)|
|M2||30 (10)||190 (50)||70 (30)||0 (0)||290 (90)|
|M3||20 (0)||310 (90)||510 (220)||40 (20)||880 (340)|
|M4||70 (20)||280 (130)||80 (50)||0 (0)||440 (200)|
|M5||0 (0)||10 (0)||20 (10)||0 (0)||30 (10)|
|M6||20 (0)||270 (70)||210 (90)||20 (10)||520 (180)|
|M7||230 (50)||2190 (600)||1040 (480)||60 (40)||3530 (1170)|
|M8||0 (0)||30 (10)||10 (10)||0 (0)||50 (10)|
|M9||20 (10)||210 (60)||110 (50)||10 (10)||350 (130)|
|M10||30 (10)||240 (70)||100 (40)||10 (10)||370 (130)|
|M11||0 (0)||0 (0)||1 (0)||0 (0)||1 (0)|
|M12||230 (50)||2420 (670)||1730 (730)||180 (110)||4560 (1560)|
|Number of events in mass range|
|M1||20 (0)||260 (60)||230 (100)||80 (60)||590 (230)|
|M2||20 (0)||210 (50)||160 (70)||50 (40)||440 (160)|
|M3||10 (0)||360 (90)||1000 (470)||240 (180)||1620 (750)|
|M4||50 (10)||300 (150)||140 (100)||30 (30)||520 (280)|
|M5||0 (0)||10 (0)||40 (20)||40 (30)||90 (50)|
|M6||20 (0)||300 (80)||430 (200)||200 (150)||960 (440)|
|M7||190 (40)||2390 (600)||2110 (930)||730 (510)||5420 (2090)|
|M8||0 (0)||30 (10)||30 (10)||10 (10)||70 (30)|
|M9||20 (0)||230 (60)||160 (70)||30 (20)||430 (160)|
|M10||30 (10)||240 (70)||100 (40)||10 (10)||370 (130)|
|M11||0 (0)||0 (0)||1 (0)||0 (0)||1 (0)|
|M12||190 (40)||2700 (680)||3710 (1690)||1830 (1380)||8440 (3790)|
Table 4 shows how the number of events detected depends on the configuration of the detector for the reference model M1. If the final LISA configuration is similar to the , four-link NGO model used in , the number of events would be about a factor of smaller. If a more sensitive configuration was launched, then the event rate could be increased by a factor of about three. The change in the number of detected events is similar for the other astrophysical models, except for the high mass compact object model (M4), for which the decrease in number of events going to the configuration is only a factor of , and the increase going to the configuration is only a factor of . This difference arises because the EMRIs in M4 are intrinsically louder and so all detector configurations can detect most of the EMRI events out to moderate redshift.
|4 links||6 links||4 links||6 links||4 links||6 links|
3.2 Parameter-estimation accuracy
Figure 1 shows the distribution of parameter-estimation precisions that we would expect for the set of EMRIs observed in each population model. We present results for the redshifted masses of the MBH and compact object, for the spin of the MBH and for the eccentricity at plunge. We have normalised results to a fixed SNR of (accuracies scale as ) and we see that after doing that the precision is essentially independent of the population model. The distribution of SNRs is determined by the spatial distribution of sources, which is uniform in comoving volume and hence also similar in each model. However, the models that predict larger numbers of events will have some sources at high SNRs which will have correspondingly more precise parameter determinations. We see that we would expect to measure the redshifted masses of the components and the orbital characteristics such as eccentricity to precisions of –, while also measuring MBH spins to –. Figure 2 shows corresponding results for EMRI measurements of source distances and sky positions. Luminosity distance should typically be determined to a few to a few tens of percent, and sky location should be measured to about a thousandth of a steradian, which is a few square degrees. Sky position will be slightly less well measured for EMRIs with larger mass inspiraling objects, since such sources are observable for a smaller fraction of a LISA orbit. The luminosity distance accuracy is also approximately the accuracy with which the intrinsic masses of the sources will be determined, since the redshift will be inferred from the luminosity distance. We would therefore expect to determine the intrinsic MBH mass to better than in the majority of cases.
These results were computed using the Schwarzschild plunge condition. Results computed with the Kerr plunge condition are typically a factor of better for the intrinsic parameters (and comparable for the extrinsic parameters), so these results are most likely conservative.
4 Implications for EMRI science
In  it was shown that the classic LISA configuration would be able to measure the slope of the mass function of MBHs in the LISA range to a precision of if EMRI events were observed. This is the precision with which the mass function is currently known. This conclusion is robust to assumptions about the configuration of the detector. In the astrophysical models presented here, we expect more than EMRI events in two years for all but the most pessimistic model. In most cases, the number of EMRIs expected is several hundreds. The precision scales as the square root of the number of events, so with a factor of more EMRIs we would measure the slope of the mass function to a precision of . This would provide constraints on the low mass end of the MBH mass function that can not be determined in any other way. We note, however, that  assumed the dependence of the EMRI rate on MBH mass was known, which is not the case due to the various complications discussed here. LISA will actually determine the product of the MBH mass function with the rate per MBH to the aforementioned precision. Combining EMRI observations with those of MBH mergers could break this degeneracy, but more work is required to properly understand this.
EMRI observations can be used to probe cosmological parameters. A GW observation cannot determine the redshift of an EMRI source and electromagnetic counterparts are unlikely to be observed. However, constraints can also be determined statistically. McLeod and Hogan  showed that by using the electromagnetically determined redshifts of all galaxies within any given LISA EMRI error box, a statistical constraint on the Hubble constant could be determined. A determination of to precision would be possible provided LISA observes EMRIs at redshift . This calculation was done for the classic LISA configuration and assumed that LISA could determine luminosity distance and sky position to precisions and . The new baseline will not be able to measure parameters to the same precision. However, these precisions will be reached for some events. In Table 5 we show how many EMRI events will be observed at redshift in each model and how many will be observed with the assumed parameter-estimation accuracy. We see that most models predict the requisite number of EMRIs, but we would only expect events with the assumed precision. We would expect a factor of worse precision with a factor of fewer events, and so it is clear that a measurement of will be possible in most cases. The LISA mission is likely to be longer than the two years assumed here, which will increase the number of accurately determined EMRIs. Moreover, all of the EMRIs will contribute something to the measurement, so assuming only the well localised EMRIs are useful is clearly conservative. Thus it is likely that EMRI observations with LISA will provide interesting cosmological constraints, unless the total number of observed EMRIs is toward the low end.
|Schwarzschild plunge condition||Kerr plunge condition|
|Model||small error||small error|
4.3 Fundamental physics
EMRIs are excellent probes of fundamental physics, as they generate many cycles in the strong field regime. In particular, EMRIs can be used to map out the space-time structure outside the central MBH, and identify any deviations from the Kerr metric structure predicted by general relativity. There has been a lot of research in this area and a full discussion of tests of fundamental physics with LISA can be found in . An example of the application of EMRIs to fundamental physics is to ask with what precision an EMRI could measure a deviation in the quadrupole moment of a black hole away from the Kerr metric value. Such a study was performed in , using a gravitational waveform model constructed by adding the leading order effect of a quadruple deviation into the AK model . We have repeated the same calculation for the current EMRI populations. We assume that EMRIs are occurring in a Kerr MBH spacetime, i.e., the true , but compute a Fisher matrix including this quadrupole deviation parameter. The corresponding element of the inverse Fisher Matrix is then a measure of how large a quadrupole deviation could be present in an observed EMRI before we would be able to identify that it was there. These results are shown in Figure 3, again normalised to a fixed SNR . We see that every EMRI should provide a constraint on deviations from the Kerr metric at a level of –, irrespective of our assumptions about the underlying astrophysical model.
We have described a comprehensive study of the prospects for detection of EMRIs with LISA. Our study has attempted to quantify, for the first time, the uncertainties in EMRI rate predictions arising from astrophysical uncertainties, as well as updating predictions for the new LISA baseline configuration used in . We find that LISA should observe several hundred EMRI events over two years, with uncertainties of about one order of magnitude in each direction. These predictions are robust to the distribution of MBH spins and the possible depletion of MBHs at low masses. For all of these events LISA will determine the intrinsic parameters to high precision (sub-percent accuracy), determine sky location to a few square degrees and determine luminosity distance to . The parameter-estimation results are largely independent of the astrophysical model assumptions.
These EMRI detections have tremendous scientific potential. EMRIs will provide new constraints on the mass function and spin distribution of MBHs in the LISA range. EMRIs will provide constraints on cosmological parameters that are competitive to existing constraints, but with completely independent systematics. Finally, EMRIs will provide strong tests of aspects of the theory of general relativity. We have shown that all of these scientific objectives will be realised irrespective of the true nature of the astrophysical population. The science will be enhanced if the number of observed EMRIs is at the high end of the plausible range, but these goals will still be met even if the rates are nearer to the low end.
AS is supported by the Royal Society. E. Barausse and E. Berti acknowledge support from the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This work has made use of the Horizon Cluster, hosted by the Institut d’Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us. E. Berti was supported by NSF Grant No. PHY-1607130 and by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Programme. PAS acknowledges support from the Ramón y Cajal Programme of the Ministry of Economy, Industry and Competitiveness of Spain. PAS’s work has been partially supported by the CAS President’s International Fellowship Initiative.
-  Abbott, B P et al. 2016 Phys. Rev. Lett. 116 061102
-  Abbott, B P et al. 2016 Phys. Rev. Lett. 116 241103
-  Abbott, B P et al. 2016 Phys. Rev. X 6 041015
-  Croton D J et al. 2006 Mon. Not. Roy. Astron. Soc. 365 11
-  The eLISA Consortium 2013 The Gravitational Universe arxiv:1375.5720
-  Armano M et al. 2016 Phys. Rev. Lett. 116 231101
-  Amaro-Seoane P et al. 2017 Laser Interferometer Space Antenna arxiv:1702.00786
-  Sesana A 2016 Phys. Rev. Lett. 116 2311021
-  Amaro-Seoane P, Gair J R, Freitag M, Miller M C, Mandel I, Cutler C and Babak S 2007 Class. Quantum Grav. 24 R113
-  Gair J R, Vallisneri M, Larson S L and Baker J G 2013 Liv. Rev. Rel. 16 7 (cited on February 28th 2017)
-  Gair J R 2009 Class. Quantum Grav. 26 094034
-  Barack L and Cutler C J 2004 Phys. Rev. D 69 082005
-  Huerta E A and Gair J R 2009 Phys. Rev. D 79 084021
-  Babak S, Gair J R, Sesana A, Barausse E, Sopuerta C F, Berry C P L, Berti E, Amaro-Seoane P, Petiteau A and Klein A 2017 Phys. Rev. D submitted, arxiv:1703.09722
-  Gair J R, Tang C and Volonteri M 2010 Phys. Rev. D 81 104014
-  Barausse E 2012 Mon. Not. Roy. Astron. Soc. 423 2533
-  Klein A et al. 2016 Phys. Rev. D 93 024003
-  Preto M and Amaro-Seoane 2010 Astrophys. J. 708 L42
-  Sesana A, Barausse E, Dotti M and Rossi E M 2014 Astrophys. J. 794 104
-  Amaro-Seoane P and Preto M 2011 Class. Quantum Grav. 28 094017
-  Merritt D 2015 Astrophys. J. 814 57
-  Gültekin K et al. 2009 Astrophys. J. 698 198
-  Graham A W and Scott N 2013 Astrophys. J. 764 151
-  Kormendy J and Ho L C 2013 Ann. Rev. Astron. & Astrophys. 51 511
-  Amaro-Seoane P, Sopuerta C F and Freitag M D 2013 Mon. Not. Roy. Astron. Soc. 429 3155
-  Poisson E, Pound A and Vega I 2011 Liv. Rev. Rel. 14 7 (cited on February 28th 2017)
-  Gair J R and Glampedakis K 2006 Phys. Rev. D 73 064037
-  Babak S, Fang H, Gair J R, Glampedakis K and Hughes S A 2007 Phys. Rev. D 75 024005
-  Gair J R et al. 2004 Class. Quantum Grav. 21 S1595
-  Babak S et al. 2010 Class. Quantum Grav. 27 084009
-  MacLeod C L and Hogan C J 2008 Phys. Rev. D 77 3512
-  Barack L and Cutler C 2007 Phys. Rev. D 75 042003