Introduction

1 Introduction

The theory of gravitational radiation from isolated sources, in the context of general relativity, is a fascinating science that can be explored by means of what was referred to in the XVIIIth century France as l’analyse sublime: The analytical (i.e., mathematical) method, and more specifically the resolution of partial differential equations. Indeed, the field equations of general relativity, when use is made of the harmonic-coordinate conditions, take the form of a quasi-linear hyperbolic differential system of equations, involving the famous wave operator or d’Alembertian [140 ]. The resolution of that system of equations constitutes a problème bien posé in the sense of Hadamard [236 *, 104 ], and which is amenable to an analytic solution using approximation methods.

Nowadays, the importance of the field lies in the exciting comparison of the theory with contemporary astrophysical observations, of binary pulsars like the historical Hulse–Taylor pulsar PSR 1913+16 [250 ], and, in a forthcoming future, of gravitational waves produced by massive and rapidly evolving systems such as inspiralling compact binaries. These should be routinely detected on Earth by the network of large-scale laser interferometers, including the advanced versions of the ground-based interferometers LIGO and VIRGO, with also GEO and the future cryogenic detector KAGRA. The first direct detection of a coalescence of two black holes has been achieved on September 14, 2015 by the advanced LIGO detector [1 ]. Further ahead, the space-based laser interferometer LISA (actually, the evolved version eLISA) should be able to detect supermassive black-hole binaries at cosmological distances.

To prepare these experiments, the required theoretical work consists of carrying out a sufficiently general solution of the Einstein field equations, valid for a large class of matter systems, and describing the physical processes of the emission and propagation of the gravitational waves from the source to the distant detector, as well as their back-reaction onto the source. The solution should then be applied to specific sources like inspiralling compact binaries.

For general sources it is hopeless to solve the problem via a rigorous deduction within the exact theory of general relativity, and we have to resort to approximation methods. Of course the ultimate aim of approximation methods is to extract from the theory some firm predictions to be compared with the outcome of experiments. However, we have to keep in mind that such methods are often missing a clear theoretical framework and are sometimes not related in a very precise mathematical way to the first principles of the theory.

The flagship of approximation methods is the post-Newtonian approximation, which has been developed from the early days of general relativity [303 *]. This approximation is at the origin of many of the great successes of general relativity, and it gives wonderful answers to the problems of motion and gravitational radiation of systems of compact objects. Three crucial applications are:

The motion of N point-like objects at the first post-Newtonian approximation level [184 *], is taken into account to describe the solar system dynamics (motion of the centers of mass of planets);
The gravitational radiation-reaction force, which appears in the equations of motion at the second-and-a-half post-Newtonian (2.5PN) order [148 *, 147 *, 143 *, 142 *], has been experimentally verified by the observation of the secular acceleration of the orbital motion of the Hulse–Taylor binary pulsar PSR 1913+16 [399 *, 400 *, 398 *];
The analysis of gravitational waves emitted by inspiralling compact binaries – two neutron stars or black holes driven into coalescence by emission of gravitational radiation – necessitates the prior knowledge of the equations of motion and radiation field up to very high post-Newtonian order.

This article reviews the current status of the post-Newtonian approach to the problems of the motion of inspiralling compact binaries and their emission of gravitational waves. When the two compact objects approach each other toward merger, the post-Newtonian expansion will lose accuracy and should be taken over by numerical-relativity computations [359 , 116 , 21 ]. We shall refer to other review articles like Refs. [233 , 187 *] for discussions of numerical-relativity methods. Despite very intensive developments of numerical relativity, the post-Newtonian approximation remains indispensable for describing the inspiral phase of compact binaries to high accuracy, and for providing a powerful benchmark against which the numerical computations are tested.

Part A of the article deals with general post-Newtonian matter sources. The exterior field of the source is investigated by means of a combination of analytic post-Minkowskian and multipolar approximations. The physical observables in the far-zone of the source are described by a specific set of radiative multipole moments. By matching the exterior solution to the metric of the post-Newtonian source in the near-zone the explicit expressions of the source multipole moments are obtained. The relationships between the radiative and source moments involve many non-linear multipole interactions, among them those associated with the tails (and tails-of-tails, etc.) of gravitational waves.

Part B is devoted to the application to compact binary systems, with particular emphasis on black hole binaries with spins. We present the equations of binary motion, and the associated Lagrangian and Hamiltonian, at the third post-Newtonian (3PN) order beyond the Newtonian acceleration. The gravitational-wave energy flux, taking consistently into account the relativistic corrections in the binary’s moments as well as the various tail effects, is derived through 3.5PN order with respect to the quadrupole formalism. The binary’s orbital phase, whose prior knowledge is crucial for searching and analyzing the signals from inspiralling compact binaries, is deduced from an energy balance argument (in the simple case of circular orbits).

All over the article we try to state the main results in a form that is simple enough to be understood without the full details; however, we also outline some of the proofs when they present some interest on their own. To emphasize the importance of some key results, we present them in the form of mathematical theorems. In applications we generally show the most up-to-date results up to the highest known post-Newtonian order.¹

1.1 Analytic approximations and wave generation formalism

The basic problem we face is to relate the asymptotic gravitational-wave form $hij$ generated by some isolated source, at the location of a detector in the wave zone of the source, to the material content of the source, i.e., its stress-energy tensor $Tαβ$ , using approximation methods in general relativity.² Therefore, a general wave-generation formalism must solve the field equations, and the non-linearity therein, by imposing some suitable approximation series in one or several small physical parameters. Some important approximations that we shall use in this article are the post-Newtonian method (or non-linear $1∕c$ -expansion), the post-Minkowskian method or non-linear iteration ( $G$ -expansion), the multipole decomposition in irreducible representations of the rotation group (or equivalently $a$ -expansion in the source radius), the far-zone expansion ( $1∕R$ -expansion in the distance to the source), and the perturbation in the small mass limit ( $ν$ -expansion in the mass ratio of a binary system). In particular, the post-Newtonian expansion has provided us in the past with our best insights into the problems of motion and radiation. The most successful wave-generation formalisms make a gourmet cocktail of these approximation methods. For reviews on analytic approximations and applications to the motion and the gravitational wave-generation see Refs. [404 , 142 *, 144 , 145 *, 405 , 421 , 46 , 52 , 378 *]. For reviews on black-hole pertubations and the self-force approach see Refs. [348 *, 373 , 177 *, 23 *].

The post-Newtonian approximation is valid under the assumptions of a weak gravitational field inside the source (we shall see later how to model neutron stars and black holes), and of slow internal motions.³ The main problem with this approximation, is its domain of validity, which is limited to the near zone of the source – the region surrounding the source that is of small extent with respect to the wavelength of the gravitational waves. A serious consequence is the a priori inability of the post-Newtonian expansion to incorporate the boundary conditions at infinity, which determine the radiation reaction force in the source’s local equations of motion.

The post-Minkowskian expansion, by contrast, is uniformly valid, as soon as the source is weakly self-gravitating, over all space-time. In a sense, the post-Minkowskian method is more fundamental than the post-Newtonian one; it can be regarded as an “upstream” approximation with respect to the post-Newtonian expansion, because each coefficient of the post-Minkowskian series can in turn be re-expanded in a post-Newtonian fashion. Therefore, a way to take into account the boundary conditions at infinity in the post-Newtonian series is to control first the post-Minkowskian expansion. Notice that the post-Minkowskian method is also upstream (in the previous sense) with respect to the multipole expansion, when considered outside the source, and with respect to the far-zone expansion, when considered far from the source.

The most “downstream” approximation that we shall use in this article is the post-Newtonian one; therefore this is the approximation that dictates the allowed physical properties of our matter source. We assume mainly that the source is at once slowly moving and weakly stressed, and we abbreviate this by saying that the source is post-Newtonian. For post-Newtonian sources, the parameter defined from the components of the matter stress-energy tensor $αβ T$ and the source’s Newtonian potential $U$ by

is much less than one. This parameter represents essentially a slow motion estimate $𝜖 ∼ v∕c$ , where $v$ denotes a typical internal velocity. By a slight abuse of notation, following Chandrasekhar et al. [122 *, 124 *, 123 *], we shall henceforth write formally $𝜖 ≡ 1∕c$ , even though $𝜖$ is dimensionless whereas $c$ has the dimension of a velocity. Thus, $1 ∕c ≪ 1$ in the case of post-Newtonian sources. The small post-Newtonian remainders will be denoted $n 𝒪 (1∕c )$ . Furthermore, still following Refs. [122 *, 124 *, 123 *], we shall refer to a small post-Newtonian term with formal order $𝒪 (1∕cn)$ relative to the Newtonian acceleration in the equations of motion, as $n2$ PN.

We have $|U ∕c2|1∕2 ≪ 1∕c$ for sources with negligible self-gravity, and whose dynamics are therefore driven by non-gravitational forces. However, we shall generally assume that the source is self-gravitating; in that case we see that it is necessarily weakly (but not negligibly) self-gravitating, i.e., $21∕2 |U ∕c | = 𝒪 (1 ∕c)$ .⁴ Note that the adjective “slow-motion” is a bit clumsy because we shall in fact consider very relativistic sources such as inspiralling compact binaries, for which $v∕c$ can be as large as 50% in the last rotations, and whose description necessitates the control of high post-Newtonian approximations.

At the lowest-order in the Newtonian limit $1∕c → 0$ , the gravitational waveform of a post-Newtonian matter source is generated by the time variations of the quadrupole moment of the source. We shall review in Section 1.2 the utterly important “Newtonian” quadrupole moment formalism [183 *, 285 *]. Taking into account higher post-Newtonian corrections in a wave generation formalism will mean including into the waveform the contributions of higher multipole moments, beyond the quadrupole. Post-Newtonian corrections of order $𝒪 (1∕cn)$ beyond the quadrupole formalism will still be denoted as $n 2$ PN. Building a post-Newtonian wave generation formalism must be concomitant to understanding the multipole expansion in general relativity.

The multipole expansion is one of the most useful tools of physics, but its use in general relativity is difficult because of the non-linearity of the theory and the tensorial character of the gravitational interaction. In the stationary case, the multipole moments are determined by the expansion of the metric at spatial infinity [219 , 238 , 384 ], while, in the case of non-stationary fields, the moments, starting with the quadrupole, are defined at future null infinity. The multipole moments have been extensively studied in the linearized theory, which ignores the gravitational forces inside the source. Early studies have extended the Einstein quadrupole formula [given by Eq. (4 *) below] to include the current-quadrupole and mass-octupole moments [332 *, 333 *], and obtained the corresponding formulas for linear momentum [332 , 333 , 30 *, 358 ] and angular momentum [339 *, 134 ]. The general structure of the infinite multipole series in the linearized theory was investigated by several works [369 , 367 , 343 , 403 *], from which it emerged that the expansion is characterized by two and only two sets of moments: Mass-type and current-type moments. Below we shall use a particular multipole decomposition of the linearized (vacuum) metric, parametrized by symmetric and trace-free (STF) mass and current moments, as given by Thorne [403 *]. The expressions of the multipole moments as integrals over the source, valid in the linearized theory but irrespective of a slow motion hypothesis, have been worked out in [309 , 119 , 118 , 154 *]. In particular, Damour & Iyer [154 *] obtained the complete closed-form expressions for the time-dependent mass and spin multipole moments (in STF guise) of linearized gravity.

In the full non-linear theory, the (radiative) multipole moments can be read off the coefficient of $1∕R$ in the expansion of the metric when $R → +∞$ , with a null coordinate $T − R ∕c = const$ . The solutions of the field equations in the form of a far-field expansion (power series in $1∕R$ ) have been constructed, and their properties elucidated, by Bondi et al. [93 *] and Sachs [368 *]. The precise way under which such radiative space-times fall off asymptotically has been formulated geometrically by Penrose [337 *, 338 *] in the concept of an asymptotically simple space-time (see also Ref. [220 *]). The resulting Bondi–Sachs–Penrose approach is very powerful, but it can answer a priori only a part of the problem, because it gives information on the field only in the limit where $R → +∞$ , which cannot be connected in a direct way to the actual matter content and dynamics of the source. In particular the multipole moments that one considers in this approach are those measured at infinity – we call them the radiative multipole moments. These moments are distinct, because of non-linearities, from some more natural source multipole moments, which are defined operationally by means of explicit integrals extending over the matter and gravitational fields.

An alternative way of defining the multipole expansion within the complete non-linear theory is that of Blanchet & Damour [57 *, 41 *], following pioneering works by Bonnor and collaborators [94 *, 95 *, 96 *, 251 *] and Thorne [403 *]. In this approach the basic multipole moments are the source moments, rather than the radiative ones. In a first stage, the moments are left unspecified, as being some arbitrary functions of time, supposed to describe an actual physical source. They are iterated by means of a post-Minkowskian expansion of the vacuum field equations (valid in the source’s exterior). Technically, the post-Minkowskian approximation scheme is greatly simplified by the assumption of a multipolar expansion, as one can consider separately the iteration of the different multipole pieces composing the exterior field.⁵ In this “multipolar-post-Minkowskian” (MPM) formalism, which is physically valid over the entire weak-field region outside the source, and in particular in the wave zone (up to future null infinity), the radiative multipole moments are obtained in the form of some non-linear functionals of the more basic source moments. A priori, the method is not limited to post-Newtonian sources; however, we shall see that, in the current situation, the closed-form expressions of the source multipole moments can be established only in the case where the source is post-Newtonian [44 *, 49 *]. The reason is that in this case the domain of validity of the post-Newtonian iteration (viz. the near zone) overlaps the exterior weak-field region, so that there exists an intermediate zone in which the post-Newtonian and multipolar expansions can be matched together. This is a standard application of the method of matched asymptotic expansions in general relativity [114 *, 113 *, 7 *, 357 *].

To be more precise, we shall show how a systematic multipolar and post-Minkowskian iteration scheme for the vacuum Einstein field equations yields the most general physically admissible solution of these equations [57 *]. The solution is specified once we give two and only two sets of time-varying (source) multipole moments. Some general theorems about the near-zone and far-zone expansions of that general solution will be stated. Notably, we show [41 *] that the asymptotic behaviour of the solution at future null infinity is in agreement with the findings of the Bondi–Sachs–Penrose [93 *, 368 *, 337 *, 338 *, 220 *] approach to gravitational radiation. However, checking that the asymptotic structure of the radiative field is correct is not sufficient by itself, because the ultimate aim, as we said, is to relate the far field to the properties of the source, and we are now obliged to ask: What are the multipole moments corresponding to a given stress-energy tensor $αβ T$ describing the source? The general expression of the moments was obtained at the level of the second post-Newtonian (2PN) order in Ref. [44 *], and was subsequently proved to be in fact valid up to any post-Newtonian order in Ref. [49 *]. The source moments are given by some integrals extending over the post-Newtonian expansion of the total (pseudo) stress-energy tensor $αβ τ$ , which is made of a matter part described by $αβ T$ and a crucial non-linear gravitational source term $αβ Λ$ . These moments carry in front a particular operation of taking the finite part ( $ℱ 𝒫$ as we call it below), which makes them mathematically well-defined despite the fact that the gravitational part $Λ αβ$ has a spatially infinite support, which would have made the bound of the integral at spatial infinity singular (of course the finite part is not added a posteriori to restore the well-definiteness of the integral, but is proved to be actually present in this formalism). The expressions of the moments had been obtained earlier at the 1PN level, albeit in different forms, in Ref. [59 *] for the mass-type moments [see Eq. (157a) below], and in Ref. [155 *] for the current-type ones.

The wave-generation formalism resulting from matching the exterior multipolar and post-Minkowskian field [57 *, 41 *] to the post-Newtonian source [44 *, 49 *] is able to take into account, in principle, any post-Newtonian correction to both the source and radiative multipole moments (for any multipolarity of the moments). The relationships between the radiative and source moments include many non-linear multipole interactions, because the source moments mix with each other as they “propagate” from the source to the detector. Such multipole interactions include the famous effects of wave tails, corresponding to the coupling between the non-static moments with the total mass $M$ of the source. The non-linear multipole interactions have been computed within the present wave-generation formalism up to the 3.5PN order in Refs. [60 *, 50 *, 48 *, 74 *, 197 *]. Furthermore, the back-reaction of the gravitational-wave emission onto the source, up to the 1.5PN order relative to the leading order of radiation reaction, has also been studied within this formalism [58 *, 43 *, 47 *]. Now, recall that the leading radiation reaction force, which is quadrupolar, occurs already at the 2.5PN order in the source’s equations of motion. Therefore the 1.5PN “relative” order in the radiation reaction corresponds in fact to the 4PN order in the equations of motion, beyond the Newtonian acceleration. It has been shown that the gravitational-wave tails enter the radiation reaction at precisely the 1.5PN relative order, i.e., 4PN absolute order [58 *]. A systematic post-Newtonian iteration scheme for the near-zone field, formally taking into account all radiation reaction effects, has been obtained, fully consistent with the present formalism [357 *, 75 *].

A different wave-generation formalism has been devised by Will & Wiseman [424 *] (see also Refs. [422 *, 335 *, 336 *]), after earlier attempts by Epstein & Wagoner [185 *] and Thorne [403 *]. This formalism has exactly the same scope as the one of Refs. [57 *, 41 *, 44 *, 49 *], i.e., it applies to any isolated post-Newtonian sources, but it differs in the definition of the source multipole moments and in many technical details when properly implemented [424 *]. In both formalisms, the moments are generated by the post-Newtonian expansion of the pseudo-tensor $αβ τ$ , but in the Will–Wiseman formalism they are defined by some compact-support integrals terminating at some finite radius enclosing the source, e.g., the radius $ℛ$ of the near zone. By contrast, in Refs. [44 *, 49 *], the moments are given by some integrals covering the whole space ( $ℝ3$ ) and regularized by means of the finite part $ℱ 𝒫$ . Nevertheless, in both formalisms the source multipole moments, which involve a whole series of relativistic corrections, must be coupled together in a complicated way in the true non-linear solution; such non-linear couplings form an integral part of the radiative moments at infinity and thereby of the observed signal. We shall prove in Section 4.3 the complete equivalence, at the most general level, between the two formalisms.

1.2 The quadrupole moment formalism

The lowest-order wave generation formalism is the famous quadrupole formalism of Einstein [183 ] and Landau & Lifshitz [285 *]. This formalism applies to a general isolated matter source which is post-Newtonian in the sense of existence of the small post-Newtonian parameter $𝜖$ defined by Eq. (1 *). However, the quadrupole formalism is valid in the Newtonian limit $𝜖 → 0$ ; it can rightly be qualified as “Newtonian” because the quadrupole moment of the matter source is Newtonian and its evolution obeys Newton’s laws of gravity. In this formalism the gravitational field $hTTij$ is expressed in a transverse and traceless (TT) coordinate system covering the far zone of the source at retarded times,⁶ as

where $R = |X |$ is the distance to the source, $T − R ∕c$ is the retarded time, $N = X ∕R$ is the unit direction from the source to the far away observer, and $1 𝒫ijab = 𝒫ia𝒫jb − 2𝒫ij𝒫ab$ is the TT projection operator, with $𝒫ij = δij − NiNj$ being the projector onto the plane orthogonal to $N$ . The source’s quadrupole moment takes the familiar Newtonian form

where $ρ$ is the Newtonian mass density. The total gravitational power emitted by the source in all directions around the source is given by the Einstein quadrupole formula

Our notation $ℱ$ stands for the total gravitational energy flux or gravitational “luminosity” of the source. Similarly, the total angular momentum flux is given by

where $𝜖abc$ denotes the standard Levi-Civita symbol with $𝜖123 = 1$ .

Associated with the latter energy and angular momentum fluxes, there is also a quadrupole formula for the radiation reaction force, which reacts on the source’s dynamics in consequence of the emission of waves. This force will inflect the time evolution of the orbital phase of the binary pulsar and inspiralling compact binaries. At the position $(x,t)$ in a particular coordinate system covering the source, the reaction force density can be written as [114 *, 113 *, 319 *]

This is the gravitational analogue of the damping force of electromagnetism. However, notice that gravitational radiation reaction is inherently gauge-dependent, so the expression of the force depends on the coordinate system which is used. Consider now the energy and angular momentum of a matter system made of some perfect fluid, say

The specific internal energy of the fluid is denoted $Π$ , and obeys the usual thermodynamic relation $dΠ = − P d(1∕ρ)$ where $P$ is the pressure; the gravitational potential obeys the Poisson equation $ΔU = − 4πG ρ$ . We compute the mechanical losses of energy and angular momentum from the time derivatives of $E$ and $Ji$ . We employ the usual Eulerian equation of motion $ρdvi∕dt = − ∂ P + ρ∂ U + F reac i i i$ and continuity equation $∂ ρ + ∂(ρvi) = 0 t i$ . Note that we add the small dissipative radiation-reaction contribution $reac Fi$ in the equation of motion but neglect all conservative post-Newtonian corrections. The result is

where one recognizes the fluxes at infinity given by Eqs. (4 *) and (5 *), and where the second terms denote some total time derivatives made of quadratic products of derivatives of the quadrupole moment. Looking only for secular effects, we apply an average over time on a typical period of variation of the system; the latter time derivatives will be in average numerically small in the case of quasi-periodic motion (see e.g., [103 ] for a discussion). Hence we obtain

where the brackets denote the time averaging over an orbit. These balance equations encode the secular decreases of energy and angular momentum by gravitational radiation emission.

The cardinal virtues of the Einstein–Landau–Lifshitz quadrupole formalism are: Its generality – the only restrictions are that the source be Newtonian and bounded; its simplicity, as it necessitates only the computation of the time derivatives of the Newtonian quadrupole moment (using the Newtonian laws of motion); and, most importantly, its agreement with the observation of the dynamics of the binary pulsar PSR 1913+16 [399 *, 400 *, 398 *]. Indeed let us apply the balance equations (9) to a system of two point masses moving on an eccentric orbit modelling the binary pulsar PSR 1913+16 – the classic references are [340 *, 339 *]; see also [186 , 415 ]. We use the binary’s Newtonian energy and angular momentum,

where $a$ and $e$ are the semi-major axis and eccentricity of the orbit and $m 1$ and $m2$ are the two masses. From the energy balance equation (9a) we obtain first the secular evolution of $a$ ; next changing from $a$ to the orbital period $P$ using Kepler’s third law,⁷ we get the secular evolution of the orbital period $P$ as

The last factor, depending on the eccentricity, comes out from the orbital average and is known as the Peters & Mathews [340 *] “enhancement” factor, so designated because in the case of the binary pulsar, which has a rather large eccentricity $e ≃ 0.617$ , it enhances the effect by a factor $∼ 12$ . Numerically, one finds $⟨dP ∕dt⟩ = − 2.4 × 10−12$ , a dimensionless number in excellent agreement with the observations of the binary pulsar [399 , 400 , 398 ]. On the other hand the secular evolution of the eccentricity $e$ is deduced from the angular momentum balance equation (9b) [together with the previous result (11 *)], as

Interestingly, the system of equations (11 *) – (12 *) can be thoroughly integrated in closed analytic form. This yields the evolution of the eccentricity [339 *]:

where $c0$ denotes an integration constant to be determined by the initial conditions at the start of the binary evolution. When $e ≪ 1$ the latter relation gives approximately $2 19∕9 e ≃ c0P$ .

For a long while, it was thought that the various quadrupole formulas would be sufficient for sources of gravitational radiation to be observed directly on Earth – as they had proved to be amply sufficient in the case of the binary pulsar. However, further works [139 *]⁸ and [87 *, 138 *] showed that this is not true, as one has to include post-Newtonian corrections to the quadrupole formalism in order to prepare for the data analysis of future detectors, in the case of inspiralling compact binaries. From the beautiful test of the orbital decay (11 *) of the binary pulsar, we can say that the post-Newtonian corrections to the “Newtonian” quadrupole formalism – which we shall compute in this article – have already received a strong observational support.

1.3 Problem posed by compact binary systems

Inspiralling compact binaries, containing neutron stars and/or black holes, are likely to become the bread-and-butter sources of gravitational waves for the detectors LIGO, VIRGO, GEO and KAGRA on ground, and also eLISA in space. The two compact objects steadily lose their orbital binding energy by emission of gravitational radiation; as a result, the orbital separation between them decreases, and the orbital frequency increases. Thus, the frequency of the gravitational-wave signal, which equals twice the orbital frequency for the dominant harmonics, “chirps” in time (i.e., the signal becomes higher and higher pitched) until the two objects collide and merge.

The orbit of most inspiralling compact binaries can be considered to be circular, apart from the gradual inspiral, because the gravitational radiation reaction forces tend to circularize the motion rapidly. This effect is due to the emission of angular momentum by gravitational waves, resulting in a secular decrease of the eccentricity of the orbit, which has been computed within the quadrupole formalism in Eq. (12 *). For instance, suppose that the inspiralling compact binary was long ago (a few hundred million years ago) a system similar to the binary pulsar system, with an orbital frequency $−4 Ω0 ≡ 2π ∕P0 ∼ 10 rad∕s$ and a rather large orbital eccentricity $e0 ∼ 0.6$ . When it becomes visible by the detectors on ground, i.e., when the gravitational wave signal frequency reaches about $f ≡ Ω ∕π ∼ 10 Hz$ , the eccentricity of the orbit should be $e ∼ 10− 6$ according to the formula (13 *). This is a very small eccentricity, even when compared to high-order relativistic corrections. Only non-isolated binary systems could have a non negligible eccentricity. For instance, the Kozai mechanism [283 *, 300 *] is one important scenario that produces eccentric binaries and involves the interaction between a pair of binaries in the dense cores of globular clusters [315 ]. If the mutual inclination angle of the inner binary is strongly tilted with respect to the outer compact star, then a resonance occurs and can increase the eccentricity of the inner binary to large values. This is one motivation for looking at the waves emitted by inspiralling binaries in non-circular, quasi-elliptical orbits (see Section 10).

The main point about modelling the inspiralling compact binary is that a model made of two structureless point particles, characterized solely by two mass parameters $ma$ and possibly two spins $Sa$ (with $a = 1, 2$ labelling the particles), is sufficient in first approximation. Indeed, most of the non-gravitational effects usually plaguing the dynamics of binary star systems, such as the effects of a magnetic field, of an interstellar medium, of the internal structure of extended bodies, are dominated by gravitational effects. The main justification for a model of point particles is that the effects due to the finite size of the compact bodies are small. Consider for instance the influence of the Newtonian quadrupole moments $Qa$ induced by tidal interaction between two neutron stars. Let $aa$ be the radius of the stars, and $r12$ be the distance between the two centers of mass. We have, for tidal moments,

where $k a$ are the star’s dimensionless (second) Love numbers [321 ], which depend on their internal structure, and are, typically, of the order unity. On the other hand, for compact objects, we can introduce their “compactness” parameters, defined by the dimensionless ratios

and equal $∼ 0.2$ for neutron stars (depending on their equation of state). The quadrupoles $Qa$ will affect the Newtonian binding energy $E$ of the two bodies, and also the emitted total gravitational flux $ℱ$ as computed using the Newtonian quadrupole formula (4 *). It is known that for inspiralling compact binaries the neutron stars are not co-rotating because the tidal synchronization time is much larger than the time left till the coalescence. As shown by Kochanek [276 ] the best models for the fluid motion inside the two neutron stars are the so-called Roche–Riemann ellipsoids, which have tidally locked figures (the quadrupole moments face each other at any instant during the inspiral), but for which the fluid motion has zero circulation in the inertial frame. In the Newtonian approximation, using the energy balance equation (9a), we find that within such a model (in the case of two identical neutron stars with same mass $m$ , compactness $K$ and Love number $k$ ), the orbital phase reads

where “const” denotes a numerical coefficient of order one, $ϕ0$ is some initial phase, and $x ≡ (Gm Ω ∕c3)2∕3$ is a standard dimensionless post-Newtonian parameter of the order of $∼ 1∕c2$ (with $Ω = 2π∕P$ the orbital frequency). The first term in the right-hand side of Eq. (16 *) corresponds to the gravitational-wave damping of two point masses without internal structure; the second term is the finite-size effect, which appears as a relative correction, proportional to $(x ∕K )5$ , to the latter radiation damping effect. Because the finite-size effect is purely Newtonian, its relative correction $∼ (x∕K )5$ should not depend on the speed of light $c$ ; and indeed the factors $1∕c2$ cancel out in the ratio $x∕K$ . However, the compactness $K$ of neutron stars is of the order of 0.2 say, and by definition of compact objects we can consider that $K$ is formally of the order of unity or one half; therefore the factor $1∕c2$ it contains in (15 *) should not be taken into account when estimating numerically the effect. So the real order of magnitude of the relative contribution of the finite-size effect in Eq. (16 *) is given by the factor $x5$ alone. This means that for compact objects the finite-size effect should roughly be comparable, numerically, to a post-Newtonian correction of magnitude $5 10 x ∼ 1 ∕c$ namely 5PN order. This is a much higher post-Newtonian order than the one at which we shall investigate the gravitational effects on the phasing formula. Using $k ∼ 1$ , $K ∼ 0.2$ and the bandwidth of detectors between 10 Hz and 1000 Hz, we find that the cumulative phase error due to the finite-size effect amounts to less that one orbital rotation over a total of $∼$ 16 000 produced by the gravitational-wave damping of two neutron stars. The conclusion is that the finite-size effects can in general be neglected in comparison with purely gravitational-wave damping effects. The internal structure of the two compact bodies is “effaced” and their dynamics and radiation depend only, in first approximation, on the masses (and possibly spins); hence this property has been called the “effacement” principle of general relativity [142 *]. But note that for non-compact or moderately compact objects (such as white dwarfs for instance) the Newtonian tidal interaction dominates over the radiation damping. The constraints on the nuclear equation of state and the tidal deformability of neutron stars which can be inferred from gravitational wave observations of neutron star binary inspirals have been investigated in Refs. [320 *, 200 , 414 ]. For numerical computations of the merger of two neutron stars see Refs. [187 , 249 ].

Inspiralling compact binaries are ideally suited for application of a high-order post-Newtonian wave generation formalism. These systems are very relativistic, with orbital velocities as high as $0.5c$ in the last rotations (as compared to $−3 ∼ 10 c$ for the binary pulsar), so it is not surprising that the quadrupole-moment formalism (2 *) – (6 *) constitutes a poor description of the emitted gravitational waves, since many post-Newtonian corrections are expected to play a substantial role. This expectation has been confirmed by measurement-analyses [139 *, 137 *, 198 *, 138 *, 393 *, 346 *, 350 *, 284 *, 157 *], which have demonstrated that the post-Newtonian precision needed to implement successfully the optimal filtering technique for the LIGO/VIRGO detectors corresponds grossly, in the case of neutron-star binaries, to the 3PN approximation, or $1∕c6$ beyond the quadrupole moment approximation. Such a high precision is necessary because of the large number of orbital rotations that will be monitored in the detector’s frequency bandwidth, giving the possibility of measuring very accurately the orbital phase of the binary. Thus, the 3PN order is required mostly to compute the time evolution of the orbital phase, which depends, via Eq. (9a), on the center-of-mass binding energy $E$ and the total gravitational-wave energy flux $ℱ$ .

In summary, the theoretical problem is two-fold: On the one hand $E$ , and on the other hand $ℱ$ , are to be computed with 3PN precision or better. To obtain $E$ we must control the 3PN equations of motion of the binary in the case of general, not necessarily circular, orbits; as for $ℱ$ it necessitates the application of a 3PN wave generation formalism. It is remarkable that such high PN approximation is needed in preparation for the LIGO and VIRGO data analyses. As we shall see, the signal from compact binaries contains the signature of several non-linear effects which are specific to general relativity. We thus have the possibility of probing, experimentally, some aspects of the non-linear structure of Einstein’s theory [84 *, 85 *, 15 *, 14 *].

1.4 Post-Newtonian equations of motion

By equations of motion we mean the explicit expression of the accelerations of the bodies in terms of the positions and velocities. In Newtonian gravity, writing the equations of motion for a system of N particles is trivial; in general relativity, even writing the equations in the case N = 2 is difficult. The first relativistic terms, at the 1PN order, were derived by Lorentz & Droste [303 ]. Subsequently, Einstein, Infeld & Hoffmann [184 *] obtained the 1PN corrections for N particles by means of their famous “surface-integral” method, in which the equations of motion are deduced from the vacuum field equations, and are therefore applicable to any compact objects (be they neutron stars, black holes, or, perhaps, naked singularities). The 1PN-accurate equations were also obtained, for the motion of the centers of mass of compact bodies, by Fock [201 ] (see also Refs. [341 , 330 ]).

The 2PN approximation was tackled by Ohta et al. [324 *, 327 *, 326 *, 325 *], who considered the post-Newtonian iteration of the Hamiltonian of N point-particles. We refer here to the Hamiltonian as a “Fokker-type” Hamiltonian, which is obtained from the matter-plus-field Arnowitt–Deser–Misner (ADM) Hamiltonian by eliminating the field degrees of freedom. The 2.5PN equations of motion were obtained in harmonic coordinates by Damour & Deruelle [148 *, 147 *, 175 *, 141 *, 142 *], building on a non-linear (post-Minkowskian) iteration of the metric of two particles initiated in Ref. [31 ]. The corresponding result for the ADM-Hamiltonian of two particles at the 2PN order was given in Ref. [169 *] (see also Refs. [375 , 376 ]). The 2.5PN equations of motion have also been derived in the case of two extended compact objects [280 *, 234 *]. The 2.5PN equations of two point masses as well as the near zone gravitational field in harmonic-coordinate were computed in Ref. [76 *].⁹

Up to the 2PN level the equations of motion are conservative. Only at the 2.5PN order does the first non-conservative effect appear, associated with the gravitational radiation emission. The equations of motion up to that level [148 , 147 *, 175 , 141 , 142 *], have been used for the study of the radiation damping of the binary pulsar – its orbital $P˙$ [142 *, 143 *, 173 ]. The result was in agreement with the prediction of the quadrupole formalism given by (11 *). An important point is that the 2.5PN equations of motion have been proved to hold in the case of binary systems of strongly self-gravitating bodies [142 *]. This is via the effacing principle for the internal structure of the compact bodies. As a result, the equations depend only on the “Schwarzschild” masses, $m1$ and $m2$ , of the neutron stars. Notably their compactness parameters $K1$ and $K2$ , defined by Eq. (15 *), do not enter the equations of motion. This has also been explicitly verified up to the 2.5PN order by Kopeikin et al. [280 *, 234 *], who made a physical computation à la Fock, taking into account the internal structure of two self-gravitating extended compact bodies. The 2.5PN equations of motion have also been obtained by Itoh, Futamase & Asada [256 *, 257 *] in harmonic coordinates, using a variant (but, much simpler and more developed) of the surface-integral approach of Einstein et al. [184 *], that is valid for compact bodies, independently of the strength of the internal gravity.

At the 3PN order the equations of motion have been worked out independently by several groups, by means of different methods, and with equivalent results:

Jaranowski & Schäfer [261 *, 262 *, 263 *], and then with Damour [162 *, 164 *], employ the ADM-Hamiltonian canonical formalism of general relativity, following the line of research initiated in Refs. [324 *, 327 *, 326 *, 325 *, 169 *];
Blanchet & Faye [69 *, 71 *, 70 *, 72 *], and with de Andrade [174 *] and Iyer [79 *], founding their approach on the post-Newtonian iteration initiated in Ref. [76 *], compute directly the equations of motion (instead of a Hamiltonian) in harmonic coordinates;
Itoh & Futamase [255 *, 253 *] (see [213 ] for a review), continuing the surface-integral method of Refs. [256 *, 257 *], obtain the complete 3PN equations of motion in harmonic coordinates, without need of a self-field regularization;
Foffa & Sturani [203 *] derive the 3PN Lagrangian in harmonic coordinates within the effective field theory approach pioneered by Goldberger & Rothstein [223 *].

It has been shown [164 *, 174 *] that there exists a unique “contact” transformation of the dynamical variables that changes the harmonic-coordinates Lagrangian of Ref. [174 *] (identical to the ones issued from Refs. [255 *, 253 *] and [203 *]) into a new Lagrangian, whose Legendre transform coincides with the ADM-Hamiltonian of Ref. [162 *]. The equations of motion are therefore physically equivalent. For a while, however, they depended on one unspecified numerical coefficient, which is due to some incompleteness of the Hadamard self-field regularization method. This coefficient has been fixed by means of a better regularization, dimensional regularization, both within the ADM-Hamiltonian formalism [163 *], and the harmonic-coordinates equations of motion [61 *]. These works have demonstrated the power of dimensional regularization and its adequateness to the classical problem of interacting point masses in general relativity. By contrast, notice that, interestingly, the surface-integral method [256 *, 257 *, 255 *, 253 *] by-passes the need of a regularization. We devote Section 6 to questions related to the use of self-field regularizations.

The effective field theory (EFT) approach to the problems of motion and radiation of compact binaries, has been extensively developed since the initial proposal [223 *] (see [206 ] for a review). It borrows techniques from quantum field theory and consists of treating the gravitational interaction between point particles as a classical limit of a quantum field theory, i.e., in the “tree level” approximation. The theory is based on the effective action, defined from a Feynman path integral that integrates over the degrees of freedom that mediate the gravitational interaction.¹⁰ The phase factor in the path integral is built from the standard Einstein–Hilbert action for gravity, augmented by a harmonic gauge fixing term and by the action of particles. The Feynman diagrams naturally show up as a perturbative technique for solving iteratively the Green’s functions. Like traditional approaches [163 *, 61 *] the EFT uses the dimensional regularization.

Computing the equations of motion and radiation field using Feynman diagrams in classical general relativity is not a new idea by itself: Bertotti & Plebanski [35 ] defined the diagrammatic tree-level perturbative expansion of the Green’s functions in classical general relativity; Hari Dass & Soni [240 ]¹¹ showed how to derive the classical energy-loss formula at Newtonian approximation using tree-level propagating gravitons; Feynman diagrams have been used for the equations of motion up to 2PN order in general relativity [324 *, 327 , 326 , 325 ] and in scalar-tensor theories [151 *]. Nevertheless, the systematic EFT treatment has proved to be powerful and innovative for the field, e.g., with the introduction of a decomposition of the metric into “Kaluza–Klein type” potentials [277 ], the interesting link with the renormalization group equation [222 *], and the systematization of the computation of diagrams [203 *].

The 3.5PN terms in the equations of motion correspond to the 1PN relative corrections in the radiation reaction force. They were derived by Iyer & Will [258 *, 259 *] for point-particle binaries in a general gauge, relying on energy and angular momentum balance equations and the known expressions of the 1PN fluxes. The latter works have been extended to 2PN order [226 ] and to include the leading spin-orbit effects [428 ]. The result has been then established from first principles (i.e., not relying on balance equations) in various works at 1PN order [260 *, 336 *, 278 *, 322 *, 254 *]. The 1PN radiation reaction force has also been obtained for general extended fluid systems in a particular gauge [43 *, 47 *]. Known also is the contribution of gravitational-wave tails in the equations of motion, which arises at the 4PN order and represents a 1.5PN modification of the radiation damping force [58 *]. This 1.5PN tail-induced correction to the radiation reaction force was also derived within the EFT approach [205 , 215 ].

The state of the art on equations of motion is the 4PN approximation. Partial results on the equations of motion at the 4PN order have been obtained in [264 *, 265 *, 266 *] using the ADM Hamiltonian formalism, and in [204 *] using the EFT. The first derivation of the complete 4PN dynamics was accomplished in [166 *] by combining the local contributions [264 *, 265 *, 266 ] with a non-local contribution related to gravitational wave tails [58 *, 43 *], with the help of the result of an auxiliary analytical self-force calculation [36 *]. The non-local dynamics of [166 *] has been transformed in Ref. [167 *] into a local Hamiltonian containing an infinite series of even powers of the radial momentum. A second computation of the complete 4PN dynamics (including the same non-local interaction as in [166 *], but disagreeing on the local interaction) was accomplished in [33 *] using a Fokker Lagrangian in harmonic coordinates. Further works [39 , 248 ] have given independent confirmations of the results of Refs. [166 *, 167 ]. More work is needed to understand the difference between the results of [166 ] and [33 ].

An important body of works concerns the effects of spins on the equations of motion of compact binaries. In this case we have in mind black holes rather than neutron stars, since astrophysical stellar-size black holes as well as super-massive galactic black holes have spins which can be close to maximal. The dominant effects are the spin-orbit (SO) coupling which is linear in spin, and the spin-spin (SS) coupling which is quadratic. For maximally spinning objects, and adopting a particular convention in which the spin is regarded as a 0.5PN quantity (see Section 11), the leading SO effect arises at the 1.5PN order while the leading SS effect appears at 2PN order. The leading SO and SS effects in the equations of motion have been determined by Barker & O’Connell [27 *, 28 *] and Kidder, Will & Wiseman [275 *, 271 *]. The next-to-leading SO effect, i.e., 1PN relative order corresponding to 2.5PN order, was obtained by Tagoshi, Ohashi & Owen [394 *], then confirmed and completed by Faye, Blanchet & Buonanno [194 *]. The results were also retrieved by two subsequent calculations, using the ADM Hamiltonian [165 *] and using EFT methods [292 *, 352 *]. The ADM calculation was later generalized to the N-body problem [241 *] and extended to the next-to-leading spin-spin effects (including both the coupling between different spins and spin square terms) in Refs. [387 *, 389 *, 388 *, 247 *, 243 *], and the next-to-next-to-leading SS interactions between different spins at the 4PN order [243 *]. In the meantime EFT methods progressed concurrently by computing the next-to-leading 3PN SS and spin-squared contributions [354 *, 356 *, 355 *, 293 *, 299 *], and the next-to-next-to-leading 4PN SS interactions for different spins [294 *] and for spin-squared [298 *]. Finally, the next-to-next-to-leading order SO effects, corresponding to 3.5PN order equivalent to 2PN relative order, were obtained in the ADM-coordinates Hamiltonian [242 *, 244 *] and in the harmonic-coordinates equations of motion [307 *, 90 *], with complete equivalence between the two approaches. Comparisons between the EFT and ADM Hamiltonian schemes for high-order SO and SS couplings can be found in Refs. [295 , 299 *, 297 *]. We shall devote Section 11 to spin effects (focusing mainly on spin-orbit effects) in black hole binaries.

So far the status of post-Newtonian equations of motion is very satisfying. There is mutual agreement between all the results obtained by means of many different approaches and techniques, whenever they can be compared: point particles described by Dirac delta-functions or extended post-Newtonian fluids; surface-integrals methods; mixed post-Minkowskian and post-Newtonian expansions; direct post-Newtonian iteration and matching; EFT techniques versus traditional expansions; harmonic coordinates versus ADM-type coordinates; different processes or variants of the self-field regularization for point particles; different ways to including spins within the post-Newtonian approximation. In Part B of this article, we present complete results for the 3.5PN equations of motion (including the 1PN radiation reaction), and discuss the conservative part of the equations in the case of quasi-circular orbits. Notably, the conservative part of the dynamics is compared with numerical results for the gravitational self-force in Section 8.4.

1.5 Post-Newtonian gravitational radiation

The second problem, that of the computation of the gravitational waveform and the energy flux $ℱ$ , has to be solved by application of a wave generation formalism (see Section 1.1). The earliest computations at the 1PN level beyond the quadrupole moment formalism were done by Wagoner & Will [416 *], but based on some ill-defined expressions of the multipole moments [185 *, 403 *]. The computations were redone and confirmed by Blanchet & Schäfer [86 *] applying the rigorous wave generation formalism of Refs. [57 *, 60 *]. Remember that at that time the post-Newtonian corrections to the emission of gravitational waves had only a purely academic interest.

The energy flux of inspiralling compact binaries was then completed to the 2PN order by Blanchet, Damour & Iyer [64 *, 224 *], and, independently, by Will & Wiseman [424 *, 422 *], using their own formalism; see Refs. [66 , 82 *] for joint reports of these calculations. The energy flux has been computed using the EFT approach in Ref. [221 ] with results agreeing with traditional methods.

At the 1.5PN order in the radiation field, appears the first contribution of “hereditary” terms, which are a priori sensitive to the entire past history of the source, i.e., which depend on all previous times up to $t → − ∞$ in the past [60 *]. This 1.5PN hereditary term represents the dominant contribution of tails in the wave zone. It has been evaluated for compact binaries in Refs. [426 , 87 *] by application of the formula for tail integrals given in Ref. [60 *]. Higher-order tail effects at the 2.5PN and 3.5PN orders, as well as a crucial contribution of tails generated by the tails themselves (the so-called “tails of tails”) at the 3PN order, were obtained in Refs [45 *, 48 *].

The 3PN approximation also involves, besides the tails of tails, many non-tail contributions coming from the relativistic corrections in the (source) multipole moments of the compact binary. Those have been almost completed in Refs. [81 *, 73 *, 80 *], in the sense that the result still involved one unknown numerical coefficient, due to the use of the Hadamard regularization. We shall review in Section 6 the computation of this parameter by means of dimensional regularization [62 *, 63 *], and shall present in Section 9 the most up-to-date results for the 3.5PN energy flux and orbital phase, deduced from the energy balance equation. In recent years all the results have been generalized to non-circular orbits, including both the fluxes of energy and angular momentum, and the associated balance equations [10 *, 9 *, 12 *]. The problem of eccentric orbits will be the subject of Section 10.

Besides the problem of the energy flux there is the problem of the gravitational waveform itself, which includes higher-order amplitude corrections and correlatively higher-order harmonics of the orbital frequency, consistent with the post-Newtonian order. Such full post-Newtonian waveform is to be contrasted with the so-called “restricted” post-Newtonian waveform which retains only the leading-order harmonic of the signal at twice the orbital frequency, and is often used in practical data analysis when searching the signal. However, for parameter estimation the full waveform is to be taken into account. For instance it has been shown that using the full waveform in the data analysis of future space-based detectors like eLISA will yield substantial improvements (with respect to the restricted waveform) of the angular resolution and the estimation of the luminosity distance of super-massive black hole binaries [16 , 17 , 410 ].

The full waveform has been obtained up to 2PN order in Ref. [82 *] by means of two independent wave generations (respectively those of Refs. [57 *, 44 *] and [424 *]), and it was subsequently extended up to the 3PN order in Refs. [11 *, 273 , 272 *, 74 *]. At that order the signal contains the contributions of harmonics of the orbital frequency up to the eighth mode. The motivation is not only to build accurate templates for the data analysis of gravitational wave detectors, but also to facilitate the comparison and match of the high post-Newtonian prediction for the inspiral waveform with the numerically-generated waveforms for the merger and ringdown. For the latter application it is important to provide the post-Newtonian results in terms of a spin-weighted spherical harmonic decomposition suitable for a direct comparison with the results of numerical relativity. Recently the dominant quadrupole mode $(ℓ,m ) = (2, 2)$ in the spin-weighted spherical harmonic decomposition has been obtained at the 3.5PN order [197 *]. Available results will be provided in Sections 9.4 and 9.5.

At the 2.5PN order in the waveform appears the dominant contribution of another hereditary effect called the “non-linear memory” effect (or sometimes Christodoulou effect) [128 *, 427 *, 406 *, 60 *, 50 *]. This effect was actually discovered using approximation methods in Ref. [42 *] (see [60 *] for a discussion). It implies a permanent change in the wave amplitude from before to after a burst of gravitational waves, which can be interpreted as the contribution of gravitons in the known formulas for the linear memory for massless particles [99 ]. Note that the non-linear memory takes the form of a simple anti-derivative of an “instantaneous” term, and therefore becomes instantaneous (i.e., non-hereditary) in the energy flux which is composed of the time-derivative of the waveform. In principle the memory contribution must be computed using some model for the evolution of the binary system in the past. Because of the cumulative effect of integration over the whole past, the memory term, though originating from 2.5PN order, finally contributes in the waveform at the Newtonian level [427 *, 11 *]. It represents a part of the waveform whose amplitude steadily grows with time, but which is nearly constant over one orbital period. It is therefore essentially a zero-frequency effect (or DC effect), which has rather poor observational consequences in the case of the LIGO-VIRGO detectors, whose frequency bandwidth is always limited from below by some cut-off frequency $fseismic > 0$ . Non-linear memory contributions in the waveform of inspiralling compact binaries have been thoroughly computed by Favata [189 *, 192 *].

The post-Newtonian results for the waveform and energy flux are in complete agreement (up to the 3.5PN order) with the results given by the very different technique of linear black-hole perturbations, valid when the mass of one of the bodies is small compared to the other. This is the test-mass limit $ν → 0$ , in which we define the symmetric mass ratio to be the reduced mass divided by the total mass, $ν ≡ μ∕m$ such that $ν = 1∕4$ for equal masses. Linear black-hole perturbations, triggered by the geodesic motion of a small particle around the black hole, have been applied to this problem by Poisson [345 *] at the 1.5PN order (following the pioneering work [216 ]), by Tagoshi & Nakamura [393 *], using a numerical code up to the 4PN order, and by Sasaki, Tagoshi & Tanaka [372 , 395 *, 397 *] (see also Ref. [316 ]), analytically up to the 5.5PN order. More recently the method has been improved and extended up to extremely high post-Newtonian orders: 14PN [209 ] and even 22PN [210 ] orders – but still for linear black-hole perturbations.

To successfully detect the gravitational waves emitted by spinning black hole binaries and to estimate the binary parameters, it is crucial to include spins effects in the templates, most importantly the spin-orbit effect which is linear in spins. The spins will affect the gravitational waves through a modulation of their amplitude, phase and frequency. Notably the orbital plane will precess in the case where the spins are not aligned or anti-aligned with the orbital angular momentum, see e.g., Ref. [8 *]. The leading SO and SS contributions in the waveform and flux of compact binaries are known from Refs. [275 *, 271 *, 314 *]; the next-to-leading SO terms at order 2.5PN were obtained in Ref. [53 *] after a previous attempt in [328 ]; the 3PN SO contribution is due to tails and was computed in Ref. [54 *], after intermediate results at the same order (but including SS terms) given in [353 ]. Finally, the next-to-next-to-leading SO contributions in the multipole moments and the energy flux, corresponding to 3.5PN order, and the next-to-leading SO tail corresponding to 4PN order, have been obtained in Refs. [89 *, 306 *]. The next-to-leading 3PN SS and spin-squared contributions in the radiation field were derived in Ref. [88 *]. In Section 11 we shall give full results for the contributions of spins (at SO linear level) in the energy flux and phase evolution up to 4PN order.

A related topic is the loss of linear momentum by gravitational radiation and the resulting gravitational recoil (or “kick”) of black-hole binary systems. This phenomenon has potentially important astrophysical consequences [313 ]. In models of formation of massive black holes involving successive mergers of smaller “seed” black holes, a recoil with sufficient velocity could eject the system from the host galaxy and effectively terminate the process. Recoils could eject coalescing black holes from dwarf galaxies or globular clusters. Even in galaxies whose potential wells are deep enough to confine the recoiling system, displacement of the system from the center could have important dynamical consequences for the galactic core.

Post-Newtonian methods are not ideally suited to compute the recoil of binary black holes because most of the recoil is generated in the strong field regime close to the coalescence [199 *]. Nevertheless, after earlier computations of the dominant Newtonian effect [30 , 199 ]¹² and the 1PN relative corrections [425 ], the recoil velocity has been obtained up to 2PN order for point particle binaries without spin [83 *], and is also known for the dominant spin effects [271 *]. Various estimations of the magnitude of the kick include a PN calculation for the inspiraling phase together with a treatment of the plunge phase [83 ], an application of the effective-one-body formalism [152 ], a close-limit calculation with Bowen–York type initial conditions [385 ], and a close-limit calculation with initial PN conditions for the ringdown phase [288 , 290 ].

In parallel the problem of gravitational recoil of coalescing binaries has attracted considerable attention from the numerical relativity community. These computations led to increasingly accurate estimates of the kick velocity from the merger along quasicircular orbits of binary black holes without spins [115 , 20 ] and with spins [117 ]. In particular these numerical simulations revealed the interesting result that very large kick velocities can be obtained in the case of spinning black holes for particular spin configurations.

	Abstract
1	Introduction
	1.1	Analytic approximations and wave generation formalism
	1.2	The quadrupole moment formalism
	1.3	Problem posed by compact binary systems
	1.4	Post-Newtonian equations of motion
	1.5	Post-Newtonian gravitational radiation
A	Post-Newtonian Sources
2	Non-linear Iteration of the Vacuum Field Equations
	2.1	Einstein’s field equations
	2.2	Linearized vacuum equations
	2.3	The multipolar post-Minkowskian solution
	2.4	Generality of the MPM solution
	2.5	Near-zone and far-zone structures
3	Asymptotic Gravitational Waveform
	3.1	The radiative multipole moments
	3.2	Gravitational-wave tails and tails-of-tails
	3.3	Radiative versus source moments
4	Matching to a Post-Newtonian Source
	4.1	The matching equation
	4.2	General expression of the multipole expansion
	4.3	Equivalence with the Will–Wiseman formalism
	4.4	The source multipole moments
5	Interior Field of a Post-Newtonian Source
	5.1	Post-Newtonian iteration in the near zone
	5.2	Post-Newtonian metric and radiation reaction effects
	5.3	The 3.5PN metric for general matter systems
	5.4	Radiation reaction potentials to 4PN order
B	Compact Binary Systems
6	Regularization of the Field of Point Particles
	6.1	Hadamard self-field regularization
	6.2	Hadamard regularization ambiguities
	6.3	Dimensional regularization of the equations of motion
	6.4	Dimensional regularization of the radiation field
7	Newtonian-like Equations of Motion
	7.1	The 3PN acceleration and energy for particles
	7.2	Lagrangian and Hamiltonian formulations
	7.3	Equations of motion in the center-of-mass frame
	7.4	Equations of motion and energy for quasi-circular orbits
	7.5	The 2.5PN metric in the near zone
8	Conservative Dynamics of Compact Binaries
	8.1	Concept of innermost circular orbit
	8.2	Dynamical stability of circular orbits
	8.3	The first law of binary point-particle mechanics
	8.4	Post-Newtonian approximation versus gravitational self-force
9	Gravitational Waves from Compact Binaries
	9.1	The binary’s multipole moments
	9.2	Gravitational wave energy flux
	9.3	Orbital phase evolution
	9.4	Polarization waveforms for data analysis
	9.5	Spherical harmonic modes for numerical relativity
10	Eccentric Compact Binaries
	10.1	Doubly periodic structure of the motion of eccentric binaries
	10.2	Quasi-Keplerian representation of the motion
	10.3	Averaged energy and angular momentum fluxes
11	Spinning Compact Binaries
	11.1	Lagrangian formalism for spinning point particles
	11.2	Equations of motion and precession for spin-orbit effects
	11.3	Spin-orbit effects in the gravitational wave flux and orbital phase
	Acknowledgments
	References
	Footnotes
	Updates
	Figures
	Tables