### 6.5 The two-body problem in general relativity

The largest effort in gravitational radiation theory in recent years has been to study the two-body problem using various approximations. The reason is that gravitationally bound binary systems are likely to be important gravitational wave sources, and until the evolution of such a system is thoroughly understood, it will not be possible to extract the maximum possible information from the observations.

From Figure 2, we see that ground-based detectors will be sensitive to compact binaries with mass in the range of while LISA will be sensitive to the mass range . As we have seen in Section 3, most classes of binary sources will follow orbits that evolve strongly due to gravitational radiation reaction. In the case of ground-based detectors, they will all merge within a year of entering the observation band. In the case of LISA, we might observe sources (both stellar mass binaries as well as SMBH binaries), whose frequency hardly changes.

In contrast to Newtonian gravity, modeling a bound binary in general relativity is complicated by the existence of gravitational radiation and the nonlinearity of Einstein’s equations. It must therefore be done approximately. The three most important approximation methods for solving gravitational wave problems are:

• The post-Newtonian scheme. This is a combination of a low-velocity expansion ( small) and a weak-field expansion ( small), in which the two small parameters are linked because a gravitationally-bound binary satisfies the virial relation , even in relativity. The zero-order solution is the Newtonian binary system. The post-Newtonian (PN) approximation has now been developed to a very high order in because the velocities in late-stage binaries, just before coalescence, are very high.
• Perturbation theory. This is an expansion in which the small parameter is the mass-ratio of the binary components. The zero-order solution is the field of the more massive component, and linear field corrections due to the second component determine the binary’s orbital motion and the emitted radiation. This approximation is fully relativistic at all orders. It is being used to study the signals emitted by compact stars and stellar-mass black holes as they fall into SMBHs, an important source for LISA.
• Numerical approaches. With numerical relativity one can in principle simulate any desired relativistic system, no matter how strong the fields or high the velocities. It is being used to study the final stage of the evolution of binaries, including their coalescence, after the PN approximation breaks down. Although it deals with fully relativistic and nonlinear general relativity, the method needs to be regarded as an approximate one, since spacetime is not resolved to infinite precision. The accuracy of a numerical simulation is normally judged by performing convergence tests, that is by doing the simulation at a variety of resolutions and showing that there are no unexpected differences between them.

We will review the physics that can be learned from models using each of these approximation schemes. But first we treat a subject that is common to all binaries that evolve due to radiation reaction, which is that one can estimate their distance from a gravitational wave observation.

#### 6.5.1 Binaries as standard candles: distance estimation

Astronomers refer to systems as standard candles if their intrinsic luminosity is known, so that when the apparent luminosity of a particular system is measured, then its distance can be deduced. As mentioned in Section 3.4.2, radiating binaries have this property, if one can measure the effects of radiation reaction on their orbits [332]. Because of the one-dimensional nature of gravitational wave data, some scientists have begun calling these standard sirens [197]. Over cosmological distances, the distance measured from the observation is the luminosity distance. We discuss in Section 8 below how this can be used to determine the Hubble constant and even the acceleration of the universe in methods independent of any cosmic distance ladder.

#### 6.5.2 Numerical approaches to the two-body problem

From the point of view of relativity, the simplest two-body problem is that of two black holes. There are no matter fields and no point particles, just pure gravity. Therefore, the physics is entirely governed by Einstein’s equations, which are highly nonlinear and rather difficult to solve. A number of teams have worked for over three decades towards developing accurate numerical solutions for the coalescence of two black holes, using fully three-dimensional numerical simulations.

A breakthrough came in early 2005 with Pretorius [300] announcing the results from the first stable simulation ever, followed by further breakthroughs by two other groups [10563] with successful simulations. The main results from numerical simulations of nonspinning black holes are rather simple. Indeed, just as the EOB had predicted, and probably contrary to what many people had expected, the final merger is just a continuation of the adiabatic inspiral, leading on smoothly to merger and ringdown. In Figure 10 we show the results from one of the numerical simulations (right panel) and that of the EOB (left panel), both for the same initial conditions. There is also good agreement in the prediction of the total energy emitted by the system, being 5.0% ( 0.4%) (for a review see [301]) and 3.1% [99], by numerical simulations and EOB, respectively, as well as the spin of the final black hole (respectively, 0.69 and 0.8) that results from the merger.

The total energy emitted and the spin angular momentum of the black hole both depend on the spin angular momenta of the parent black holes and how they are aligned with respect to the orbital angular momentum. In the test-mass limit, it is well known that the last stable orbit of a test particle in prograde orbit will be closer to, and that of a retrograde orbit will be farther from, the black hole as compared to the Schwarzschild case. Thus, prograde orbits last longer and radiate more compared to retrograde orbits. The same is true even in the case of spinning black holes of comparable masses; the emitted energy will be greater when the spins are aligned with the orbital angular momentum and least when they are anti-aligned. For instance, for two equal mass black holes, each with its spin angular momentum equal to 0.76, the total energy radiated in the aligned (anti-aligned) case is 6.7% (2.2%) and the spin of the final black hole is 0.89 (0.44) [106296]. Heuristically, in the aligned case the black holes experience a repulsive force, deferring the merger of the two bodies to a much later time than in the anti-aligned case, where they experience an attractive force, accelerating the merger.

Detailed comparisons [13828590] show that we should be able to deploy the analytical templates from EOB [102136139140] (and other approximants [25]) that better fit the numerical data in our searches. With the availability of merger waveforms from numerical simulations and analytical templates, it will now be possible to search for compact binary coalescences with a greater sensitivity. The visibility of the signal improves significantly for binaries with their component masses in the range . Currently, an effort is underway to evaluate how to make use of numerical relativity simulations in gravitational wave searches [276], which should help to increase the distance reach of interferometric detectors by a factor of two and correspondingly nearly an order-of-magnitude increase in event rate.

Numerical relativity simulations have now greatly matured, allowing a variety of different studies. Some are studying the effect of the spin orientations of the component black holes on the linear momentum carried away by the final black hole, fancifully called kicks [1909464106181296]; some have focused on the dependence of the emitted waveform phase and energy on the mass ratio; and yet others have strived to evolve the system with high accuracy and for a greater number of cycles so as to push the techniques of numerical relativity to the limit [9190].

Of particular interest are the numerical values of black hole kicks that have been obtained for certain special configurations of the component spins. Velocities as large as 4000 km s–1 have been reported by several groups, but such velocities are only achieved when both black holes have large spins. Such velocities are in excess of escape velocities typical of normal galaxies and are, therefore, of great astronomical significance. These high velocities, however, are not seen for generic geometries of the initial spin orientations; therefore, their astronomical significance is not yet clear.

What is the physics behind kicks? Beamed emission of radiation from a binary could result in imparting a net linear momentum to the final black hole. The radiation could be beamed either because the masses of the two black holes are not the same (resulting in asymmetric emission in the orbital plane) or because of the precession of the orbital plane arising from spin-orbit and spin-spin interactions, or both. In the case of black holes with unequal masses, the largest kick one can get is around 170 km s–1, corresponding to a mass ratio of about 3:1. It was really with the advent of numerical simulations that superkicks begin to be realized, but only when black holes had large spins. The spin-orbit configurations that produce large kicks are rather unusual and at first sight unexpected. When the component black holes are both of the same mass and have equal but opposite spin angular momenta that lie in the orbital plane, frame dragging can lead to tilting and oscillation of the orbital plane, which, in the final phases of the evolution, could result in a rather large kick [301]. SMBHs are suspected to have large spins and, therefore, the effect of spin on the evolution of a binary and the final spin and kick velocity could be of astrophysical interest too.

Curiously, a recent optical observation of a distant quasar, SDSS J0927 12.65+294344.0, could well be the first identification of a superkick, causing the SMBH to escape from the parent galaxy [223]. From a fundamental physics point of view, kicks offer a new way of testing frame dragging in the vicinity of black holes, but much work is needed in this direction.

More recently, there has been an effort to understand and predict [100311137] the spin of the final black hole, which should help in further exploring interesting regions of the spin parameter space. In the relatively simple case of two black holes with equal and aligned spins of magnitude , but unequal masses, with the symmetric mass ratio being , Rezzolla et al. [311] have obtained an excellent fit for the final spin of the black hole by enforcing basic constraints from the test-mass limit:

where , , , , and . The top and middle panels of Figure 11 compare as functions of black hole spin and the symmetric mass ratio the goodness of their fit (blue short-dashed line, top panels) with the predictions of numerical simulations (circles and stars) from different groups (AEI [313], FAU–Jena [253], Jena [77] and Goddard [102]). Their residuals (red dotted lines, bottom panels) are less than a percent over the entire parameter space observed. These figures also show the fits obtained for the equal-mass but variable-spin case (green long-dashed line, left panel) [100] and for the nonspinning but unequal-mass case (green long-dashed line, middle panel) [137].

For the simple case of two equal mass black holes with aligned spins, the above analytical formula predicts that minimal and maximal final spin values of and , respectively [311]. More interestingly, one can now ask what initial configurations of the mass ratios and spins would lead to the formation of a Schwarzschild black hole (i.e., [201], which defines the boundary of the region on one side of which lie systems for which the spin of the final black hole flips relative to the initial total angular momentum (bottom panel in Figure 11).

Finally, the evolution of binaries composed of nonspinning bodies is characterized by a single parameter, namely the ratio of the masses of the two black holes. The study of systems with different mass ratios has allowed relativists to fit numerical waveforms with phenomenological waveforms [25]. The advantage of the latter waveforms is that one is able to more readily carry out data analysis in any part of the parameter space without needing the numerical data over the entire signal manifold.

Numerical relativity is still in its infancy and the parameter space is quite large. In the coming years more accurate simulations should become available, allowing the computation of waveforms with more cycles and less systematic errors. However, the challenge remains to systematically explore the effect of different spin orientations, mass ratios and eccentricity. One area that has not been explored using perturbative methods or post-Newtonian theory is that of intermediate–mass-ratio inspirals. These are systems with moderate mass ratios of order 100:1, where neither black-hole perturbation theory nor post-Newtonian approximation might be adequate. Yet, the prospect for detecting such systems in ground and space-based detectors is rather high. Numerical relativity simulations might be the only way to set up effectual search templates for such systems.

#### 6.5.3 Post-Newtonian approximation to the two-body problem

For the interpretation of observations of neutron-star–binary coalescences, which might be detected within five years by upgraded detectors that are now taking data, it is necessary to understand their orbital evolution to a high order in the PN expansion. The first effects of radiation reaction are seen at 2.5 PN order (i.e., at order beyond Newtonian gravity), but we probably have to have control in the equations of motion over the expansion at least to 3.5 PN order beyond the first radiation reaction (i.e., to order beyond Newtonian dynamics). There are many approaches to this, and we can not do justice here to the enormous effort that has gone into this field in recent years and refer the reader to the Living Reviews by Blanchet [81] and by Futamase & Itoh [169].

Most work on this problem so far has treated a binary system as if it were composed of two point masses. This is, strictly speaking, inconsistent in general relativity, since the masses should form black holes of finite size. Blanchet, Damour, Iyer, and collaborators [80] have avoided this problem by a method that involves generalized functions. They first expand in the nonlinearity parameter, and, when they have reached sufficiently high order, they obtain the velocity expansion of each order. By ordering terms in the post-Newtonian manner they have developed step-by-step the approximations up to 3.5 PN order.

A different team, led by Will, works with a different method of regularizing the point-particle singularity and compares its results with those of Blanchet et al. at each order [84]. There is no guarantee that either method can be continued successfully to any particular order, but so far they have worked well and are in agreement. Their results form the basis of the templates that are being designed to search for binary coalescences.

An interesting way of extending the validity of the expansion that is known to any order is to use Padé approximants [134135] (rational polynomials) of the fundamental quantities in the theory, namely the orbital energy and the gravitational wave luminosity. This has worked rather well in improving the convergence of PN theory. Buonanno and Damour [9899] have proposed an EOB approach to two-body dynamics, which makes it possible to compute the orbit of the binary and hence the phasing of the gravitational waves emitted beyond the last stable orbit into the merger and ringdown phases in the evolution of the black hole binary. This analytical approach has been remarkably successful and gained a lot of ground after the recent success in numerical relativity (see Section 6.5.2).

Other methods have been applied to this problem. Futamase [168] introduced a limit that combines the nonlinearity and velocity expansions in different ways in different regions of space, so that the orbiting bodies themselves have a regular (finite relativistic self-gravity) limit, while their orbital motion is treated in a Newtonian limit. This should not fail at any order [169], and has demonstrated its robustness by arriving at the same results as the other approaches, at least through 3 PN order. But it has a degree of arbitrariness in choosing initial data (see [330]) that could cause problems for gravitational wave search templates that integrate orbits for a long period of time.

Linear calculations of point particles around black holes are of interest in themselves and also for checking results of the full two-body calculations. These are well-developed for certain situations, e.g., [356261]. But the general equation of motion for such a body, taking into account all nongeodesic effects, has not yet been cast into a form suitable for practical calculations [107303]. This field is reviewed by two separate Living Reviews [294322].

Matched filtering, discussed in Section 5.1, is a plausible method of testing the validity of different approaches to computing the inspiral and merger waveforms from binary systems. Though a single observation is not likely to settle the question as to which methods are correct, a catalogue of events will help to evaluate the accuracy of different approaches by studying the statistics of the SNRs they measure.

##### 6.5.3.1 Post-Newtonian expansions of energy and luminosity.
Post-Newtonian calculations yield the expansion of the gravitational binding energy and the gravitational wave luminosity as a function of the post-Newtonian expansion parameter . This is related to the frequency of the dominant component of gravitational waves emitted by the binary system by

where is the total mass of the system. The expansions for a circular binary are [828381]

where is Euler’s constant.

##### 6.5.3.2 Evolution equation for the orbital phase.
Starting from these expressions, one requires that gravitational radiation comes at the expense of the binding energy of the system (see, e.g., [135]):

the energy balance equation. This can then be used to compute the (adiabatic) evolution of the various quantities as a function of time. For instance, the rate of change of the orbital velocity ( being the total mass) can be computed using:
where . Supplemented with a differential equation for ,
one can solve for the evolution of the system’s orbital velocity. Similarly, the evolution of the orbital phase can be computed using

##### 6.5.3.3 Phasing formulas.
The foregoing evolution equations for the orbital phase can be solved in several equivalent ways [135], each correct to the required post-Newtonian order, but numerically different from one another. For instance, one can retain the rational polynomial in Equation (118) and solve the two differential equations numerically, thereby obtaining the time evolution of . Alternatively, one might re-expand the rational function as a polynomial in , truncate it to order (where is the order to which the luminosity is given), thereby obtaining a parametric representation of the phasing formula in terms of polynomial expressions in :

where and are a reference phase and time, respectively. The standard post-Newtonian phasing formula goes one step further and inverts the second of the relations above to express as a polynomial in (again truncated to appropriate order), which is then substituted in the first of the expressions above to obtain a phasing formula as an explicit function of time:
In the above formulas and , being the time at which the two stars merge together and the gravitational wave frequency formally diverges.

##### 6.5.3.4 Waveform polarizations.
The post-Newtonian formalism also gives the two polarizations and as multipole expansions in powers of the parameter . To lowest order, the two polarizations of the radiation from a binary with a circular orbit, located at a distance , with total mass and symmetric mass ratio , are given by

where is the inclination of the orbital plane with the line of sight and is the velocity parameter introduced earlier. An interferometer will record a certain combination of the two polarizations given by
where the beam pattern functions and are those discussed in Section 4.2.1. In the case of ground-based instruments, the signal duration is pretty small, at most 15 min for neutron star binaries and smaller for heavier systems. Consequently, one can assume the source direction to be unchanging during the course of observation and the above combination produces essentially the same functional form of the waveforms as in Equation (122). Indeed, it is quite straightforward to show that
where
The factor is a function of the various angles and lies in the range [0, 1] with an RMS value of 2/5 (see Section 4.2.1, especially the discussion following Equation (62)).

These waveforms form the basis for evaluating the science that can be extracted from future observations of neutron star and black hole binaries. We will discuss the astrophysical and cosmological measurements that are made possible with such high precision waveforms in several sections that follow (6.5.5 and 8.3). It is clear from the expressions for the waveform polarizations that, at the lowest order, the radiation from a binary is predominantly emitted at twice the orbital frequency. However, even in the case of quasi-circular orbits the waves come off at other harmonics of the orbital frequency. As we shall see below, these harmonics are very important for estimating the parameters of a binary, although they do not seem to contribute much to the SNR of the system.

#### 6.5.4 Measuring the parameters of an inspiraling binary

The issue of parameter estimation in the context of black hole binaries has received a lot of attention [115129165662956746]. Most authors have used the covariance matrix for this purpose, although Markov Chain Monte Carlo (MCMC) techniques have also been used occasionally [117319318121], especially in the context of LISA [367122123124126]. Covariance matrix is often the preferred method, as one can explore a large parameter space without having to do expensive Monte Carlo simulations. However, when the parameter space is large, covariance matrix is not a reliable method for estimating parameter accuracies, especially at low SNRs [6566372]; but at high SNRs, as in the case of SMBH binaries in LISA, the problem might be that our waveforms are not accurate enough to facilitate a reliable extraction of the source parameters [130]. Although MCMC methods can give more reliable estimates, they suffer from being computationally extremely expensive. However, they are important in ascertaining the validity of results based on the covariance matrix, at least in a small subset of the parameter space, and should probably be employed in assessing parameter accuracies of candidate gravitational wave events.

In what follows we shall summarize the most recent work on parameter estimation in ground and space-based detectors for binaries with and without spin and the improvements brought about by including higher harmonics.

##### 6.5.4.1 Ground-based detectors – nonspinning components.
In Figure 12 we have plotted the one-sigma uncertainty in the measurement of the time of coalescence, chirpmass and symmetric mass ratio for initial and advanced LIGO and VIRGO [46]. The plots show errors for sources all producing a fixed SNR of 10 (left panels) or all at a fixed distance of 300 Mpc (right panels). The fractional error in chirpmass, even at a modest SNR of 10, can be as low as a few parts in ten thousand for stellar mass binaries, but the error stays around 1%, even for heavier systems that have only a few cycles in a detector’s sensitivity band. Error in the mass ratio is not as small, increasing to 100% at the higher end of the mass range explored. Thus, although the chirpmass can be measured to a good accuracy, poor estimation of the mass ratio means that the individual masses of the binary cannot be measured very well. Note also that the time of coalescence of the signal is determined pretty well, which means that we would be able to measure the location of the system in the sky quite well. At a given SNR the accuracy is better in the case of low-mass binaries, since they spend a longer duration and a greater number of cycles in the detector band and the chirpmass can be determined better than the mass ratio, since to first order the frequency evolution of a binary is determined only by the chirpmass.

##### 6.5.4.2 Measuring the parameters of supermassive black hole binaries in LISA.
In the case of LISA, the merger of SMBHs produces events with extremely large SNRs, even at a redshift of (100s to several thousands depending on the chirpmass of the source). Therefore, one expects to measure the parameters of a merger event in LISA to a phenomenal accuracy. Figure 13 depicts the distribution of the errors for a binary consisting of two SMBHs of masses at a redshift of  [232]. The distribution was obtained for ten thousand samples of the system corresponding to random orientations of the binary at random sky locations with the starting frequency greater than 3 × 10–5 Hz and the ending frequency corresponding to the last stable orbit. Each plot in Figure 13 shows the results of computations for binaries consisting of black holes with and without spins. Even in the absence of spin-induced modulations in the waveform, the parameter accuracies are pretty good. Note that spin-induced modulations in the waveform enable a far better estimation of parameters, chirpmass accuracy improving by more than an order of magnitude and reduced mass accuracy by two orders of magnitude. It is because of such accurate measurements that it will be possible to use SMBH mergers to test general relativity in the strong field regime of the theory (see below).

Although Figure 13 corresponds to a binary with specific masses, the trends shown are found to be true more generically for other systems too, the actual parameter accuracies and improvements due to spin both depending on the specific system studied.

#### 6.5.5 Improvement from higher harmonics

The results discussed so far use the restricted post-Newtonian approximation in which the waveform polarizations contain only twice the orbital frequency, neglecting all higher-order corrections (including those to the second harmonic). The full waveform is a post-Newtonian expansion of the two polarizations as a power-series in and consists of terms that have not only the dominant harmonic at twice the orbital frequency, but also other harmonics of the waveform. Schematically, the full waveform can be written as [81374]

where is the symmetric mass ratio, the first sum (index ) is over the different harmonics of the waveform and the second sum (index ) is over the different post-Newtonian orders. Note that post-Newtonian order weighs down the importance of higher-order amplitude corrections by an appropriate factor of the small parameter . In the restricted post-Newtonian approximation one keeps only the lowest-order term. Since happens to be zero, the dominant term corresponds to and , containing twice the orbital frequency.

The various signal harmonics, and the associated additional structure in the waveform, can potentially enhance our ability to measure the parameters of a binary to a greater accuracy. The reason we can expect to do so can be seen by looking at the spectra of gravitational waves with and without these harmonics. For a binary that is oriented face on with respect to a detector only the second harmonic is seen, while for any other orientation the radiation is emitted at all other harmonics, the influence of the harmonics becoming more pronounced as the inclination angle changes from 0 to . Figure 14 compares, in the frequency band of ground-based detectors, the spectrum of a source using the restricted post-Newtonian approximation (left panel) to the full waveform. In both cases the source is inclined to the line of sight at 45 degrees.

Following is a list of improvements brought about by higher harmonics. In the case of ground-based detectors Van Den Broeck and Sengupta [374375] found that, when harmonics are included, the SNR hardly changes, but is always smaller, relative to a restricted waveform. However, the presence of frequencies higher than twice the orbital frequency means that it will be possible to observe heavier systems, increasing the mass reach of ground-based detectors by a factor of 2 to 3 in advanced LIGO and third generation detectors [374375]. The same effect was found in the case of LISA too, allowing LISA to observe SMBH masses up to a  [49]. More than the increased mass reach, the harmonics reduce the error in the estimation of the chirpmass, symmetric mass ratio and the time of arrival by more than an order of magnitude for stellar-mass black hole binaries. The same is true to a greater extent in the case of SMBH binaries, allowing as well a far greater accuracy in the measurement of the luminosity distance and sky resolution in LISA’s observation of these sources [50364]. For instance, Figure 15 [364] shows the gain in LISA’s angular resolution for two massive black-hole–binary mergers as a consequence of using higher harmonics for a specific orientation of the binary. Improvements of order 10 to 100 can be seen over large regions of the sky. This improved performance of LISA makes it a good probe of dark energy [50] (see Section 8.3).

A word of caution is in order with regard to the improvements brought about by higher harmonics. If the sensitivity of a detector has an abrupt lower frequency cutoff, or falls off rapidly below a certain frequency, then the harmonics bring about a more dramatic improvement than when the sensitivity falls off gently. Higher harmonics, nevertheless, always help in reducing the random errors associated with the measurement of parameters of a coalescing black-hole binary.