{'fraction_non_alphanumeric': 0.037737517687487364, 'fraction_numerical': 0.011837982615726703, 'mean_word_length': 4.39924965893588, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Post-Newtonian theory and the two-body problem

Reliable predictions of general relativity theory are extracted using approximation methods. Among these, the powerful post-Newtonian approximation provides us with our best insights into the problems of motion and gravitational radiation of systems of compact objects. This approximation has reached an impressive mature status, because of important progress regarding its theoretical foundations, and the successful construction of templates of gravitational waves emitted by inspiralling compact binaries. The post-Newtonian predictions are routinely used for searching and analyzing the very weak signals of gravitational waves in current generations of detectors. High-accuracy comparisons with the results of numerical simulations for the merger and ring-down of binary black holes are going on. In this article we give an overview on the general formulation of the post-Newtonian approximation and present up-to-date results for the templates of compact binary inspiral.


Introduction
Although relativists admire the mathematical coherence -and therefore beauty -of Einstein's general relativity, this theory is not easy to manage when drawing firm predictions for the outcome of laboratory experiments and astronomical observations. Indeed only few exact solutions of the Einstein field equations are known, and one is obliged in most cases to resort to approximation methods. It is beyond question that approximation methods in general relativity do work, and in some cases with some incredible precision. Many of the great successes of general relativity were in fact obtained using approximation methods. However because of the complexity of the field equations, such methods become awfully intricate at high approximation orders. On the other hand it is difficult to set up a formalism in which the approximation method would be perfectly well-defined and based on clear premises. Sometimes it is impossible to relate the approximation method to the exact framework of the theory. In this case the only thing one can do is to rely on the approximation method as the only representation of real phenomena, and to discover "empirically" that the approximation works well.
The most important approximation scheme in general relativity is the post-Newtonian expansion, which can be viewed as an expansion when the speed of light c tends to infinity, and is physically valid under the assumptions of weak gravitational field inside the source and of slow internal motion. The post-Newtonian approximation makes sense only in the near-zone of an isolated matter source, defined as r ≪ λ , where λ ≡ c T is the wavelength of the emitted gravitational radiation, with T being a characteristic time scale of variation of the source. The approximation has been formalized in the early days of general relativity by Einstein [1], De Sitter [2,3], Lorentz & Droste [4]. It was subsequently developed notably by Einstein, Infeld & Hoffmann [5], Fock [6,7], Plebanski & Bazanski [8], Chandrasekhar and collaborators [9,10,11], Ehlers and his school [12,13,14], Papapetrou and coworkers [15,16].
Several long-standing problems with the post-Newtonian approximation for general isolated slowly-moving systems have hindered progress untill recently. At high post-Newtonian ordners some divergent Poisson-type integrals appear, casting a doubt on the physical soundness of the approximation. Linked to that, the domain of validity of the post-Newtonian approximation is limited to the near-zone of the source, making it a priori difficult to incorporate into the scheme the condition of no-incoming radiation, to be imposed at past null infinity from an isolated source. In addition, from a mathematical point of view, we do not know what the "reliability" of the post-Newtonian series is, i.e. if it comes from the Taylor expansion of a family of exact solutions.
The post-Newtonian approximation gives wonderful answers to the problems of motion and gravitational radiation, two of general relativity's corner stones. Three crucial applications are: 1. The motion of N point-like objects at the first post-Newtonian level [5] i.e. 1PN, 1 is taken into account to describe the Solar System dynamics (motion of the centers of mass of planets); 2. The gravitational radiation-reaction force, which appears in the equations of motion at the 2.5PN order [17,18,19,20], has been experimentally verified by the observation of the secular acceleration of the orbital motion of the Hulse-Taylor binary pulsar PSR 1913+16 [21,22,23]; 3. The analysis of gravitational waves emitted by inspiralling compact binariestwo neutron stars or black holes driven into coalescence by emission of gravi-tational radiation -necessitate the prior knowledge of the equations of motion and radiation field up to high post-Newtonian order.
Strategies to detect and analyze the very weak signals from inspiralling compact binaries involve matched filtering of a set of accurate theoretical template waveforms against the output of the detectors. Measurement-accuracy analyses have shown that, in order to get sufficiently accurate theoretical templates, one must at least include conservative post-Newtonian effects up to the 3PN level, and radiationreaction effects up to 5.5PN order, i.e. 3PN beyond the leading order of radiation reaction which is 2.5PN [24,25,26,27,28,29]. The appropriate description of inspiralling compact binaries is by two structureless point-particles, characterized solely by their masses m 1 or m 2 (and possibly their spins). 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, the influence of the internal structure of the compact bodies, are dominated by purely gravitational effects. Inspiralling compact binaries are very clean systems which can essentially be described in pure general relativity.
Although point-particles are ill-defined in the exact theory, they are admissible in post-Newtonian approximations. Furthermore the model of point particles can be pushed to high post-Newtonian order, where an a priori more realistic model involving the internal structure of compact bodies would fail through becoming intractable. However there is an important worry: a process of regularization to deal with the infinite self-field of point particles is crucially needed. The regularization should be carefully defined to be implemented at high orders. It should hopefully be followed by a renormalization.
The orbit of inspiralling compact binaries can be considered to be circular, apart from the gradual inspiral, with an excellent approximation. Indeed, gravitational radiation reaction forces tend to circularize rapidly the orbital motion. At each instant during the gradual inspiral, the eccentricity e of the orbit is related to the instantaneous frequency ω ≡ 2π/P by [30] e 2 ≃ const ω − 19/9 (for e ≪ 1) .
For instance one can check that the eccentricity of a system like the binary pulsar PSR 1913+16 will be e ≃ 5 10 −6 when the gravitational waves become visible by the detectors, i.e. when the signal frequency after its long chirp reaches f ≡ ω/π ≃ 10 Hz. Only systems formed very late, near their final coalescence, could have a non-negligible residual eccentricity. The intrinsic rotations or spins of the compact bodies could play an important role, yielding some relativistic spin-orbit and spin-spin couplings, both in the binary's equations of motion and radiation field. The spin of a rotating body is of the order of S ∼ m a v spin , where m and a denote the mass and typical size of the body, and where v spin represents the velocity of the body's surface. In the case of compact bodies we have a ∼ G m/c 2 , and for maximal rotation v spin ∼ c; for such objects the magnitude of the spin is roughly S ∼ G m 2 /c. It is thus customary to introduce a dimensionless spin parameter, generally denoted by χ, defined by We have χ ≤ 1 for black-holes, and χ 0.63 − 0.74 for neutron stars (depending on the equation of state of nuclear matter inside the neutron star). For binary pulsars such as PSR 1913+16, we have χ 5 10 −3 . Considering models of evolution of observed binary pulsar systems when they become close to the coalescence we expect that the spins will make a negligible contribution to the accumulated phase in this case. However, astrophysical observations suggest that black holes can have nonnegligible spins, due to spin up driven by accretion from a companion during some earlier phase of the binary evolution. For a few black holes surrounded by matter, observations indicate a significant intrinsic angular momentum and the spin may even be close to its maximal value. However, very little is known about the blackhole spin magnitudes in binary systems. For black holes rotating near the maximal value the templates of gravitational waves need to take into account the effects of spins, both for a successful detection and an accurate parameter estimation. We devote Part A of this article, including Sections 2 to 7, to a general overview of the formulation of the post-Newtonian approximation for isolated sources. We emphasize that the approximation can be carried out up to any post-Newtonian order, without the aforesaid problem of divergences. The main technique used is the matching of asymptotic expansions which permits to obtain a complete post-Newtonian solution incorporing the correct boundary conditions at infinity. Then Part B, i.e. Sections 8 to 13, will deal with the application to systems composed of compact objects. The subtle issues linked with the self-field regularization of pointparticles are discussed in some details. The results for the two-body equations of motion and radiation field at the state-of-the-art 3PN level are presented in the case of circular orbits appropriate to inspiralling compact binaries. We put the accent on the description of spinning particles and particularly on the spin-orbit coupling effect on the binary's internal energy and gravitational-wave flux.

A Post-Newtonian formalism 2 Einstein field equations
General relativity is based on two independent tenets, the first one concerned with the dynamics of the gravitational field, the second one dealing with the coupling of all the matter fields with the gravitational field. Accordingly the action of general relativity is made of two terms, The first term represents the kinetic Einstein-Hilbert action for gravity and tells that the gravitational field g αβ propagates like a pure spin-2 field. Here R is the Ricci scalar and g = det(g αβ ) is the determinant of the metric. The second term expresses the fact that all matter fields (collectively denoted by Ψ ) are minimally coupled to the metric g αβ which defines the physical lengths and times as measured in local laboratory experiments. The field equations are obtained by varying the action with respect to the metric (such that δ g αβ = 0 when |x µ | → ∞) and form a system of ten differential equations of second order, The Einstein tensor G αβ ≡ R αβ − 1 2 R g αβ is generated by the matter stress-energy tensor T αβ ≡ (2/ √ −g)(δ S matter /δ g αβ ). Four equations give, via the contracted Bianchi identity, the conservation equation for the matter system as which must be solved conjointly with the Einstein field equations for the gravitational field. The matter equation (5) reads also Let us introduce an asymptotically Minkowskian coordinate system such that the gravitational-wave amplitude, defined by h αβ ≡ √ −g g αβ − η αβ , 2 is divergenceless, i.e. satisfies the de Donder or harmonic gauge condition With this coordinate choice the Einstein field equations can be recast into the d'Alembertian equation Here = η µν ∂ µ ∂ ν is the usual (flat-spacetime) d'Alembertian operator. The source term τ αβ can rightly be interpreted as the "effective" stress-energy distribution of the matter and gravitational fields in harmonic coordinates; note that τ αβ is not a tensor and we shall often call it a pseudo-tensor. It is conserved in the sense that ∂ µ τ α µ = 0 .
This is equivalent to the condition of harmonic coordinates (7) and to the covariant conservation (5) of the matter tensor. The pseudo-tensor is made of the contribution of matter fields described by T αβ , and of the gravitational contribution Λ αβ which is a complicated functional of the gravitational field variable h µν and its first and second derivatives. Thus, The point is that Λ αβ is at least quadratic in the field strength h µν , so the field equations (8) are naturally amenable to a perturbative non-linear treatment. The general expression is To select a physically sensible solution of the field equations in the case of a bounded matter system, we must impose a boundary condition at infinity, namely the famous no-incoming radiation condition, which ensures that the system is truly isolated from other bodies in the Universe. In principle the no-incoming radiation condition is to be formulated at past null infinity I − . Here, we can simplify the formulation by taking advantage of the presence of the Minkowski background η αβ to define the no-incoming radiation condition with respect to Minkowskian past null infinity say I − η . Within approximate methods this is legitimate as we can view the gravitational field as propagating on the flat background η αβ ; indeed η αβ does exist at any finite order of approximation.
The no-incoming radiation condition should be such that it suppresses at I − η any homogeneous (regular in R 4 ) solution of the d'Alembertian equation h hom = 0 in a neighborhood of I − η . We have at our disposal the Kirchhoff formula which expresses h hom at some field point (x ′ ,t ′ ) in terms of its values and its derivatives on a sphere centered on x ′ with radius ρ ≡ |x ′ − x| and at retarded time t ≡ t ′ − ρ/c, where dΩ is the solid angle spanned by the unit direction N ≡ (x − x ′ )/ρ. From this formula we deduce the no-incoming radiation condition as the following limit at I − η , i.e. when r ≡ |x| → +∞ with t + r/c = const, 3 lim Now if h αβ satisfies (13), so does the pseudo-tensor τ αβ built on it, and then it is clear that the retarded integral of τ αβ does satisfy the same condition. Thus we infer that the unique solution of the Einstein equation (8) reads where the retarded integral takes the standard form Notice that since τ αβ depends on h µν and its derivatives, the equation (14) is to be viewed as an integro-differential equation equivalent to the Einstein equation (8) with no-incoming radiation condition.

Post-Newtonian iteration in the near zone
In this Section we proceed with the post-Newtonian iteration of the field equations in harmonic coordinates in the near zone of an isolated matter distribution. We have in mind a general hydrodynamical fluid, whose stress-energy tensor is smooth, i.e. T αβ ∈ C ∞ (R 4 ). Thus the scheme a priori excludes the presence of singularities; these will be dealt with in later Sections. Let us remind that the post-Newtonian approximation in "standard" form (e.g. [31,13,14]) is plagued with some apparently inherent difficulties, which crop up at some high post-Newtonian order like 3PN. Up to the 2.5PN order the approximation can be worked out without problems, and at the 3PN order the problems can generally be solved for each case at hand; but the problems worsen at higher orders. Historically these difficulties, even appearing at higher approximations, have cast a doubt on the actual soundness, from a theoretical point of view, of the post-Newtonian expansion. Practically speaking, they posed the question of the reliability of the approximation, when comparing the theory's predictions with very precise experimental results. In this Section and the next one we assess the nature of these difficulties -are they purely technical or linked with some fundamental drawback of the approximation scheme? -and eventually resolve them.
We first distinguish the problem of divergences in the post-Newtonian expansion: in higher approximations some divergent Poisson-type integrals appear. Recall that the post-Newtonian expansion replaces the resolution of an hyperboliclike d'Alembertian equation by a perturbatively equivalent hierarchy of elliptic-like Poisson equations. Rapidly it is found during the post-Newtonian iteration that the right-hand-sides of the Poisson equations acquire a non-compact support (it is distributed all over space), and that the standard Poisson integral diverges because of the bound of the integral at spatial infinity, i.e. r ≡ |x| → +∞, with t = const.
The divergencies are linked to the fact that the post-Newtonian expansion is actually a singular perturbation, in the sense that the coefficients of the successive powers of 1/c are not uniformly valid in space, typically blowing up at spatial infinity like some positive powers of r. We know for instance that the post-Newtonian expansion cannot be "asymptotically flat" starting at the 2PN or 3PN level, depending on the adopted coordinate system [32]. The result is that the standard Poisson integrals are in general badly-behaving at infinity. Trying to solve the post-Newtonian equations by means of the Poisson integral does not a priori make sense. This does not mean that there are no solutions to the problem, but simply that the Poisson integral does not constitute the good solution of the Poisson equation in the context of post-Newtonian expansions. So the difficulty is purely of a technical nature, and will be solved once we succeed in finding the appropriate solution to the Poisson equation.
We shall now prove (following [33]) that the post-Newtonian expansion can be indefinitely iterated without divergences. 4 Let us denote by means of an overline the formal (infinite) post-Newtonian expansion of the field inside the source's nearzone, which is of the form The n-th post-Newtonian coefficient is naturally the factor of the n-th power of 1/c; however, we know [36] that the post-Newtonian expansion also involves some logarithms of c, included for convenience here into the definition of the coefficients h n . For the stress-energy pseudo-tensor (10) we have the same type of expansion, The expansion starts with a term of order c 2 corresponding to the rest mass-energy (τ αβ has the dimension of an energy density). Here we shall always understand the infinite sums such as ( The second term comes from the split of the d'Alembertian into a Laplacian and a second time derivative: = ∆ − 1 c 2 ∂ 2 t . This term is zero when n = 2 and 3. We proceed by induction, i.e. fix some post-Newtonian order n, assume that we succeeded in constructing the sequence of previous coefficients p h for p ≤ n − 1, and from this show how to infer the next-order coefficient n h.
To cure the problem of divergencies we introduce a generalized solution of the Poisson equation with non-compact support source, in the form of an appropriate finite part of the usual Poisson integral obtained by regularization of the bound at infinity by means of a specific process of analytic continuation. For any source term like n τ, we multiply it by the "regularization" factor where B ∈ C is a complex number and r 0 denotes an arbitrary length scale. Only then do we apply the Poisson integral, which therefore defines a certain function of B.
The well-definedness of that integral heavily relies on the behavior of the integrand at the bound at infinity. There is no problem with the vicinity of the origin inside the source because of the smoothness of the pseudo-tensor. Then one can prove [33] that the latter function of B generates a (unique) analytic continuation down to a neighborhood of the value of interest B = 0, except at B = 0 itself, at which value it admits a Laurent expansion with multiple poles up to some finite order. Then, we consider the Laurent expansion of that function when B → 0 and pick up the finite part, or coefficient of the zero-th power of B, of that expansion. This defines our generalized Poisson integral: The integral extends over all three-dimensional space but with the latter finite-part regularization at infinity denoted FP B=0 . The main properties of our generalized Poisson operator is that it does solve the Poisson equation, namely and that the so defined solution ∆ −1 n τ owns the same properties as its source n τ, i.e. the smoothness and the same type of behavior at infinity.
The most general solution of the Poisson equation (18) will be obtained by application of the previous generalized Poisson operator to the right-hand-side of (18), and augmented by the most general homogeneous solution of the Poisson equation. Thus, we can write The last term represents the general solution of the Laplace equation which is regular at the origin r ≡ |x| = 0. It can be written, using the symmetric-trace-free (STF) language, as a multipolar series of terms of the typex L , 5 and multiplied by some 5 Here L = i 1 · · ·i ℓ denotes a multi-index composed of ℓ multipolar spatial indices i 1 , · · · , i ℓ (ranging from 1 to 3); x L ≡ x i 1 · · ·x i ℓ is the product of ℓ spatial vectors x i ; ∂ L = ∂ i 1 · · · ∂ i ℓ is the product of ℓ partial derivatives ∂ i = ∂ /∂ x i ; in the case of summed-up (dummy) multi-indices L, we do not write the ℓ summations from 1 to 3 over their indices; the STF projection is indicated with a hat, i.e.x L ≡ STF[x L ] and similarly∂ L ≡ STF[∂ L ], or sometimes using brackets surrounding the indices, STF-tensorial functions of time n B L (t). These functions will be associated with the radiation reaction of the field onto the source; they will depend on which boundary conditions are to be imposed on the gravitational field at infinity from the source. It is now trivial to iterate the process. We substitute for n−2 h in the right-handside of (22) the same expression but with n replaced by n − 2, and similarly come down until we stop at either one of the coefficients 0 h = 0 or 1 h = 0. At this point n h is expressed in terms of the "previous" p τ's and p B L 's with p ≤ n − 2. To finalize the process we introduce what we call the operator of the "instantaneous" potentials −1 I . Our notation is chosen to contrast with the standard operators of the retarded and advanced potentials −1 R and −1 A , see (15). However, beware of the fact that unlike −1 R,A the operator −1 I will be defined only when acting on a post-Newtonian series such as τ. Indeed, we pose where ∆ −k−1 is the k-th iteration of the operator (20). It is readily checked that in this way we have a solution of the source-free d'Alembertian equation, On the other hand, the homogeneous solution in (22) will yield by iteration an homogeneous solution of the d'Alembertian equation which is necessarily regular at the origin. Hence it should be of the anti-symmetric type, i.e. be made of the difference between a retarded multipolar wave and the corresponding advanced wave. We shall therefore introduce a new definition for some STF-tensorial functions A L (t) parametrizing those advanced-minus-retarded free waves. It will not be difficult to relate the post-Newtonian expansion of A L (t) to the functions n B L (t) which were introduced in (22). Finally the most general post-Newtonian solution, iterated ad infinitum and without any divergences, is obtained into the form . (25) We shall refer to the A L (t)'s as the radiation-reaction functions. If we stay at the level of the post-Newtonian iteration which is confined into the near zone we cannot do more than (25); there is no means to compute the radiation-reaction functions A L (t). We are here touching the second problem faced by the standard post-Newtonian approximation.

Post-Newtonian expansion calculated by matching
As we now understand this problem is that of the limitation to the near zone. Indeed the post-Newtonian expansion assumes that all retardations r/c are small, so it can be viewed as a formal near-zone expansion when r → 0, valid only in the region surrounding the source that is of small extent with respect to the wavelength of the emitted radiation: r ≪ λ . As we have seen, a consequence is that the post-Newtonian coefficients blow up at infinity, when r → +∞. It is thus not possible, a priori, to implement within the post-Newtonian scheme the physical information that the matter system is isolated from the rest of the Universe. The no-incoming radiation condition imposed at past null infinity I − η cannot be taken into account, a priori, within the scheme.
The near-zone limitation can be circumvented to the lowest post-Newtonian orders by considering retarded integrals that are formally expanded when c → +∞ as series of "instantaneous" Poisson-like integrals [31]. This procedure works well up to the 2.5PN level and has been shown to correctly fix the dominant radiation reaction term at the 2.5PN order [13,14]. Unfortunately such a procedure assumes fundamentally that the gravitational field, after expansion of all retardations r/c → 0, depends on the state of the source at a single time t, in keeping with the instantaneous character of the Newtonian interaction. However, we know that the post-Newtonian field (as well as the source's dynamics) will cease at some stage to be given by a functional of the source parameters at a single time, because of the imprint of gravitational-wave tails in the near zone field, in the form of some modification of the radiation reaction force at the 1.5PN relative order [37,38]. Since the reaction force is itself of order 2.5PN this means that the formal post-Newtonian expansion of retarded Green functions is no longer valid starting at the 4PN order.
The solution of the problem resides in the matching of the near-zone field to the exterior field, a solution of the vacuum equations outside the source which has been developed in previous works using some post-Minkowskian and multipolar expansions [36,39]. In the case of post-Newtonian sources, the near zone, i.e. r ≪ λ , covers entirely the source, because the source's radius itself is such that a ≪ λ . Thus the near zone overlaps with the exterior zone where the multipole expansion is valid. Matching together the post-Newtonian and multipolar-post-Minkowskian solutions in this overlapping region is an application of the method of matched asymptotic expansions, which has frequently been applied in the present context, both for radiation-reaction [40,41,37,38] and wave-generation [42,43,44] problems.
In the previous Section we obtained the most general solution (25) for the post-Newtonian expansion, as parametrized by the set of unknown radiation-reaction functions A L (t). We shall now impose the matching condition telling that the multipole decomposition of the post-Newtonian expansion h of the inner field, agrees with the near-zone expansion of the multipole expansion M (h) of the external field. Here the calligraphic letter M stands for the multipole de-composition or far-zone expansion, while the overbar denotes the post-Newtonian or near-zone expansion. The matching equation results from the numerical equality h = M (h), clearly verified in the exterior part of the near-zone, namely our overlapping region a < r ≪ λ . The left-hand-side is expanded when r → +∞ yielding M (h) while the right-hand-side is expanded when r → 0 leading to M (h). The matching equation is thus physically justified only for post-Newtonian sources, for which the exterior near-zone exists. It is actually a functional identity; it identifies, term-by-term, two asymptotic expansions, each of them being formally taken outside its own domain of validity. In the present context, the matching equation insists that the infinite far-zone expansion (r → ∞) of the inner post-Newtonian field is identical to the infinite near-zone expansion (r → 0) of the exterior multipolar field.
Let us now state that (26), plus the condition of no-incoming radiation, permits determining all the unknowns of the problem: i.e., at once, the external multipolar decomposition M (h) and the radiation-reaction functions A L and hence the inner post-Newtonian expansion h. When applied to a multipole expansion such as that of the pseudo-tensor, i.e. M (τ αβ ), we have to define a special type of generalized inverse d'Alembertian operator, built on the standard retarded integral (15) Like (15) this integral extends over the whole three-dimensional space, but a regularization factor | x ′ | B given by (19) has been "artificially" introduced for application of the finite part operation FP B=0 . The reason for introducing such regularization is to cure the divergencies of the integral when |x ′ | → 0; these are coming from the fact that the multipolar expansion is singular at the origin. We notice that this regularization factor is the same as the one entering the generalized Poisson integral (20), however its role is different, as it takes care of the bound at |x ′ | = 0 rather than at infinity. We easily find that this new object is a particular retarded solution of the wave equation, Therefore M (h αβ ) should be given by that solution plus a retarded homogeneous solution of the d'Alembertian equation (imposing the no-incoming radiation condition), i.e. be of the type Now the matching equation (26) will determine both the A L 's in (25) and the F L 's in (29). We summarize the results which have been obtained in [33]. The functions F L (t) will play an important role in the following, because they appear as the multipole moments of a general post-Newtonian source as seen from its exterior near zone. Their closed-form expression obtained by matching reads Again the integral extends over all space but the bound at infinity (where the post-Newtonian expansion becomes singular) is regularized by means of the same finite part. The z-integration involves a weighting function δ ℓ (z) defined by The integral of that function is normalized to one: +1 −1 dz δ ℓ (z) = 1. Furthermore it approaches the Dirac function in the limit of large multipoles: lim ℓ→∞ δ ℓ (z) = δ (z). The multipole moments (30) are physically valid only for post-Newtonian sources. As such, they must be considered only in a perturbative post-Newtonian sense. With the result (30) the multipole expansion (29) is fully determined and will be exploited in the next Section.
Concerning the near-zone field (25) we find that the radiation-reaction functions A L are composed of the multipole moments F L which will also characterize "linearorder" radiation reaction effects starting at 2.5PN order, and of an extra contribution R L which will be due to non-linear effets in the radiation reaction and turn out to arise at 4PN order. Thus, where F L is given by (30) and where R L is defined from the multipole expansion of the pseudo-tensor as Here the regularization deals with the bound of the integral at |x| = 0. Since the variable z extends up to infinity these functions truly depend on the whole pasthistory of the source. The weighting function therein is simply given by γ ℓ (z) ≡ −2δ ℓ (z), this definition being motivated by the fact that the integral of that function is normalized to one: +∞ 1 dz γ ℓ (z) = 1. 6 The specific contributions due to R L in the post-Newtonian metric (25) are associated with tails of waves [37,38]. The fact that the external multipolar expansion M (τ) is the source term for the function R L , and therefore will enter the expression of the near-zone metric (25), is a result of the matching condition (26) and reflects of course the no-incoming radiation condition imposed at I − η . The post-Newtonian metric (25) is now fully determined. However let us now derive an interesting alternative formulation of it [45]. To this end we introduce still another object which will be made of the expansion of the standard retarded integral (15) when c → ∞, but acting on a post-Newtonian source term τ, (34) Each of the terms is regularized by means of the finite part to deal with the bound at infinity where the post-Newtonian expansion is singular. This regularization is crucial and the object should carefully be distinguished from the "global" solution −1 (14) and in which the pseudo-tensor is not expanded in post-Newtonian fashion. We emphasize that (34) constitutes merely the definition of a (formal) post-Newtonian expansion, each term of which being built from the post-Newtonian expansion of the pseudo-tensor. Such a definition is of interest because it corresponds to what one would intuitively think as the "natural" way of performing the post-Newtonian iteration, i.e. by Taylor expanding the retardations as in [31]. Moreover, each of the terms of the series (34) is mathematically well-defined thanks to the finite part, and can therefore be implemented in practical computations. The point is that (34) solves, in a post-Newtonian sense, the wave equation, so constitutes a good prescription for a particular solution of the wave equation -as legitimate a prescription as (23). Therefore (23) and (34) should differ by an homogeneous solution of the wave equation which is necessarily of the antisymmetric type. Detailed investigations yield in which the homogeneous solution is parametrized precisely by the multipolemoment functions F L (t). This formula is the basis of our writing of the new form of the post-Newtonian expansion. Indeed, by combining (25) and (36), we nicely get which is our final expression for the general solution of the post-Newtonian field in the near-zone of any isolated matter distribution. This expression is probably the most convenient and fruitful when doing practical applications. We recognize in the first term of (37) (notwithstanding the finite part therein) the old way of performing the post-Newtonian expansion as it was advocated by Anderson & DeCanio [31]. For computations limited to the 3.5PN order, i.e. up to the level of the 1PN correction to the radiation reaction force, such first term is sufficient. However, at the 4PN order there is a fundamental breakdown of this scheme and it becomes necessary to take into account the second term in (37) which corresponds to non-linear radiation reaction effects associated with tails.
Note that the post-Newtonian solution, in either form (25) or (37), has been obtained without imposing the condition of harmonic coordinates in an explicit way, see (7). We have simply matched together the post-Newtonian and multipolar expansions, satisfying the "relaxed" Einstein field equations (8) in their respective domains, and found that the matching determines uniquely the solution. An important check (carried out in [33,45]) is therefore to verify that the harmonic coordinate condition (7) is indeed satisfied as a consequence of the conservation of the pseudotensor (9), so that we really grasp a solution of the full Einstein field equations.

Multipole moments of a post-Newtonian source
The multipole expansion of the field outside a general post-Newtonian source has been obtained in the previous Section as 7 where the multipole moments are explicitly given by (30), and the second piece reflects the non-linearities of the Einstein field equations and reads To write the latter expression we have used the fact that since the matter tensor T αβ has a spatially compact support we have M (T αβ ) = 0. Thus u αβ is indeed generated by the non-linear gravitational source term (11). We notice that the divergence of this piece, say w α ≡ ∂ µ u α µ , is a retarded homogeneous solution of the wave equation, i.e. of the same type as the first term in (38). Now from w α we can construct a secondary object v αβ which is also a retarded homogeneous solution of the wave equation, and furthermore whose divergence satisfies ∂ µ v α µ = −w α , so that it will cancel the divergence of u αβ (see [44] for details). With the above construction of v αβ we are able to define the following combination, which will constitute the linearized approximation to the multipolar expansion M (h αβ ) outside the source. Then we have Having singled out such linearized part, it is clear that the sum of the second and third terms should represent the non-linearities in the external field. If we index those non-linearities by Newton's constant G, then we can prove indeed that u αβ + v αβ = O(G 2 ). More precisely we can decompose u αβ + v αβ as a complete nonlinearity or "post-Minkowskian" expansion of the type One can effectively define a post-Minkowskian "algorithm" [36,44] able to construct the non-linear series up to any post-Minkowskian order m. The post-Minkowskian expansion represents the most general solution of the Einstein field equations in harmonic coordinates valid in the vacuum region outside an isolated source. The above linearized approximation h (1) solves the linearized vacuum Einstein field equations in harmonic coordinates and it can be decomposed into multipole moments in a standard way [49]. Modulo an infinitesimal gauge transformation preserving the harmonic gauge, namely where the infinitesimal gauge vector ϕ α (1) satisfies ϕ α (1) = 0, we can decompose 8 This decomposition defines two types of multipole moments, both assumed to be STF: the mass-type I L (u) and the current-type J L (u). These moments can be arbitrary functions of the retarded time u ≡ t − r/c, except that the monopole and dipole moments (having ℓ ≤ 1) satisfy standard conservation laws, namely The gauge transformation vector admits a decomposion in similar fashion, The six sets of STF multipole moments I L , J L , W L , X L , Y L and Z L will collectively be called the multipole moments of the source. They contain the full physical information about any isolated source as seen from its exterior near zone. Actually it should be clear that the main moments are I L and J L because the other moments W L , · · · , Z L parametrize a linear gauge transformation and thus have no physical implications at the linearized order. However because the theory is covariant with respect to non-linear diffeomorphisms and not merely with respect to linear gauge transformations, the moments W L , · · · , Z L do play a physical role starting at the non-linear level. We shall occasionaly refer to the moments W L , X L , Y L and Z L as the gauge moments.
To express in the best way the source multipole moments, we introduce the following notation for combinations of components of the pseudo-tensor τ αβ , where τ ii ≡ δ i j τ i j . Here the overbar reminds us that we are exclusively dealing with post-Newtonian-expanded expressions, i.e. formal series of the type (17). Then the general expressions of the "main" source multipole moments I L and J L in the case of the time-varying moments for which ℓ ≥ 2, are The integrands are computed at the spatial point x and at time u + z|x|/c, where u = t −r/c is the retarded time at which are evaluated the moments. We recall that z is the argument of the function δ ℓ defined in (31). Similarly we can write the expressions of the gauge-type moments W L , · · · , Z L . Notice that the source multipole moments (48) have no invariant meaning; they are defined for the harmonic coordinate system we have chosen.
Of what use are these results for the multipole moments I L and J L ? From (44) these moments parametrize the linearized metric h (1) which is the "seed" of an infinite post-Minkowskian algorithm symbolized by (42). For a specific application, i.e. for a specific choice of matter tensor like the one we shall describe in Section 8, the expressions (48) have to be worked out up to a given post-Newtonian order. The moments should then be inserted into the post-Minkowskian series (42) for the computation of the non-linearities. The result will be in the form of a non-linear multipole decomposition depending on the source moments I L , J L , · · · , Z L , say In the next Section we shall expand this metric at (retarded) infinity from the source in order to obtain the observables of the gravitational radiation field.

Radiation field and polarisation waveforms
The asymptotic waveform at future null infinity from an isolated source is the transverse-traceless (TT) projection of the metric deviation at the leading order 1/R in the distance R = |X| to the source, in a radiative coordinate system X µ = (c T, X). 9 The waveform can be uniquely decomposed [49] into radiative multipole components parametrized by mass-type moments U L and current-type ones V L . We shall define the radiative moments in such a way thay they agree with the ℓ-th time derivatives of the source moments I L and J L at the linear level, i.e.
At the non-linear level the radiative moments will crucially differ from the source moments; the relations between these two types of moments will be discussed in the next Section. The radiative moments U L (U) and V L (U) are functions of the retarded time U ≡ T − R/c in radiative coordinates.
The asymptotic waveform at distance R and retarded time U is then given by 9 Radiative coordinates T and X, also called Bondi-type coordinates [50], are such that the metric coefficients admit an expansion when R → +∞ with U ≡ T − R/c being constant, in simple powers of 1/R, without the logarithms of R plaguing the harmonic coordinate system. Here U is a null or asymptotically null characteristic. It is known that the "far-zone" logarithms in harmonic coordinates can be removed order-by-order by going to radiative coordinates [51].
We denote by N = X/R = (N i ) the unit vector pointing from the source to the faraway detector. The TT projection operator reads P i jkl = P ik P jl − 1 2 P i j P kl where P i j = δ i j − N i N j is the projector orthogonal to the unit direction N. We project out the asymptotic waveform (51) on polarization directions in a standard way. We denote the two unit polarisation vectors by P and Q, which are orthogonal and transverse to the direction of propagation N (hence P i j = P i P j + Q i Q j ). Our conventions and choice for P and Q will be specified in Section 12. Then the two "plus" and "cross" polarisation states of the waveform are Although the multipole decomposition (51) entirely describes the waveform, it is also important, especially having in mind the comparison between the post-Newtonian results and numerical relativity [52], to consider separately the modes (ℓ, m) of the waveform as defined with respect to a basis of spin-weighted spherical harmonics. To this end we decompose h + and h × as (see e.g. [52,53]) where the spin-weighted spherical harmonics of weight −2 is a function of the spherical angles (Θ , Φ) defining the direction of propagation N and reads Here k 1 = max(0, m − 2) and k 2 = min(ℓ + m, ℓ − 2). Using the orthonormality properties of these harmonics we obtain the separate modes h ℓm from the surface integral (with the overline denoting the complex conjugate) On the other hand, we can also write h ℓm directly in terms of the radiative multipole moments U L and V L , with result where U ℓm and V ℓm are the radiative mass and current moments in non-STF guise. These are given in terms of the STF moments by Here α ℓm L denotes the STF tensor connecting together the usual basis of spherical harmonics Y ℓm to the set of STF tensorsN L ≡ STF(N L ), recalling that Y ℓm andN L represent two basis of an irreducible representation of weight ℓ of the rotation group. They are related byN with the STF tensorial coefficient being The decomposition in spherical harmonic modes is especially useful if some of the radiative moments are known to higher post-Newtonian order than others. In this case the comparison with the numerical calculation [52,53] can be made for these individual modes with higher post-Newtonian accuracy.

Radiative moments versus source moments
The basis of our computation is the general solution of the Einstein field equations outside an isolated matter system computed iteratively in the form of a post-Minkowskian or non-linearity expansion (49) (see details in [36,39]). Here we give some results concerning the relation between the set of radiative moments {U L , V L } and the sets of source moments {I L , J L } and gauge moments {W L , · · · , Z L }. Complete results up to 3PN order are available and have recently been used to control the 3PN waveform of compact binaries [54]. Armed with definitions for all those moments, we proceed in a modular way. We express the radiative moments {U L , V L } in terms of some convenient intermediate constructs {M L , S L } called the canonical moments. Essentially these canonical moments take into account the effect of the gauge transformation present in (43). Therefore they differ from the source moments {I L , J L } only at non-linear order. We shall see that in terms of a post-Newtonian expansion the canonical and source moments agree with each other up to 2PN order. The canonical moments are then connected to the actual source multipole moments {I L , J L } and {W L , · · · , Z L }. The point of the above strategy is that the source moments (including gauge moments) admit closed-form expressions as integrals over the stress-energy distribution of matter and gravitational fields in the source, as shown in (48).
The mass quadrupole moment U i j (having ℓ = 2) is known up to the 3PN order [55]. At that order it is made of quadratic and cubic non-linearities, and we have j a (u) Notice the quadratic tail integral at 1.5PN order, the cubic tail-of-tail integral at 3PN order, and the non-linear memory integral at 2.5PN order [56,57,58,39]. The tail is composed of the coupling between the mass quadrupole moment M i j and the mass monopole moment or total mass M; the tail-of-tail is a coupling between M i j and two monopoles M × M; the non-linear memory is a coupling M i j × M kl . All these "hereditary" integrals imply a dependence of the waveform on the complete history of the source, from infinite past up to the current retarded time U ≡ T − R/c. The constant u 0 in the tail integrals is defined by u 0 ≡ r 0 /c, where r 0 is the arbitrary length scale introduced in (19). Note that the dominant hereditary integral is the tail arising at 1.5PN order in all radiative moments. For general ℓ we have at that order where the constants κ ℓ and π ℓ are given by Now it can be proved that the retarded time U in radiative coordinates reads where (t, r) are the harmonic coordinates. Inserting U into (61) we obtain the radiative moments expressed in terms of local source-rooted coordinates (t, r), e.g.
This no longer depends on the constant u 0 -i.e. the u 0 gets replaced by the retardation time r/c. More generally it can be checked that u 0 always disappears from physical results at the end. On the other hand we can be convinced that the constant κ ℓ (and π ℓ as well) depends on the choice of source-rooted coordinates (t, r). For instance if we change the harmonic coordinate system (t, r) to some "Schwarzschildlike" coordinates (t ′ , r ′ ) such that t ′ = t and r ′ = r + GM/c 2 , we get a new constant κ ′ ℓ = κ ℓ + 1/2. Thus we have κ 2 = 11/12 in harmonic coordinates [as shown in (60)], but κ ′ 2 = 17/12 in Schwarzschild coordinates. We still have to relate the canonical moments {M L , S L } to the source multipole moments. As we said the difference between these two types of moments comes from the gauge transformation (46) and arises only at the small 2.5PN order. The consequence is that we have to worry about this difference only for high post-Newtonian waveforms. For the mass quadrupole moment M i j , the requisite correction is given by where I i j denotes the source mass quadrupole, and where W is the monopole corresponding to the gauge moments W L (i.e. the moment having ℓ = 0). Up to 3PN order, W is only needed at Newtonian order. The expression (65) is valid in a mass-centred frame defined by the vanishing of the conserved mass dipole moment: I i = 0. Note that closed-form formulas generalizing (65) and similar expressions to all post-Newtonian orders (and all multipole interactions), if they exist are not known; they need to be investigated anew for specific cases. Thus it is convenient in the present approach to systematically keep the source moments {I L , J L , W L , X L , Y L , Z L } as the fundamental variables describing the source.

B Inspiralling compact binaries 8 Stress-energy tensor of spinning particles
So far the post-Newtonian formalism has been developed for arbitrary matter distributions. We want now to apply it to material systems made of compact objects (neutron stars or black holes) which can be described with great precision by point masses. We thus discuss the modelling of point-particles possibly carrying some intrinsic rotation or spin. This means finding the appropriate stress-energy tensor which will have to be inserted into the general post-Newtonian formulas such as the expressions of the source moments (48).
In the general case the stress-energy tensor will be the sum of a "monopolar" piece, which is a linear combination of monopole sources, i.e. made of Dirac delta-functions, plus the "dipolar" or spin piece, made of gradients of Dirac deltafunctions. Hence we write The monopole part takes the form of the stress-energy tensor for N particles (labelled by A = 1, · · · , N) without spin, reading in a four-dimensional picture Here δ (4) is the four-dimensional Dirac function. The world-line of particle A, denoted y α A , is parametrized by the particle's proper time τ A . The four-velocity is given by u α denotes the metric at the particle's location. The four-vector p α A is the particle's linear momentum. For particles without spin we shall simply have p α A = m A u α A . However for spinning particles p α A will differ from that and include some contributions from the spins as given by (75) below.
The dipolar part of the stress-energy tensor depends specifically on the spins, and reads in the classic formalism of spinning particles (due to Tulczyjew [59,60], Trautman [61], Dixon [62], Bailey & Israel [63]), where ∇ µ is the covariant derivative, and the anti-symmetric tensor S αβ A represents the spin angular momentum of particle A. In this formalism the momentum-like quantity p α A [entering (67)] is a time-like solution of the equation where D/dτ A denotes the covariant proper time derivative. The equation of translational motion of the spinning particle, equivalent to the covariant conservation ∇ µ T α µ = 0 of the total stress-energy tensor (66), is the Papapetrou equation [64] involving a coupling to curvature, The Riemann tensor is evaluated at the particle's position A. This equation can also be derived directly from an action principle [63]. It is well-known that a choice must be made for a supplementary spin condition in order to fix unphysical degrees of freedom associated with an arbitrariness in the definition of the spin tensor S αβ . This arbitrariness can be interpreted, in the case of extended bodies, as a freedom in the choice for the location of the center-of-mass worldline of the body, with respect to which the angular momentum is defined (see e.g. [65]). An elegant spin condition is the covariant one, which allows a natural definition of a spin four-(co)vector S A α such that For the spin vector S A α itself, we can choose a four-vector which is purely spatial in the particle's instantaneous rest frame, where u α A = (1, 0). Therefore, we deduce that in any frame As a consequence of the covariant spin condition (71), we easily verify that the spin scalar is conserved along the trajectories, i.e.
Furthermore, we can check, using (71) and also the law of motion (70), that the mass defined by m 2 A c 2 = −p µ A p A µ is also constant along the trajectories: m A = const. Finally, the relation linking the four-momentum p α A and the four-velocity u α A is readily deduced from the contraction of (69) with the four-momentum, which results in where (pu) A ≡ p A µ u µ A . Contracting further this relation with the four-velocity one deduces the expression of (pu) A and inserting this back into (75) yields the desired relation between p α A and u α A . Focusing our attention on spin-orbit interactions, which are linear in the spins, we can neglect quadratic and higher spin corrections denoted O(S 2 ); drastic simplifications of the formalism occur in this case. Since the right-hand-side of (75) is quadratic in the spins, we find that the four-momentum is linked to the four-velocity by the simple proportionality relation Hence, the spin condition (71) becomes Also, the equation of evolution for the spin, sometimes called the precessional equation, follows immediately from the relationship (69) together with the law (76) as Hence the spin vector S A α satisfies the equation of parallel transport, which means that it remains constant in a freely falling frame, as could have been expected beforehand. Of course the norm of the spin vector is preserved along trajectories, In Section 13 we shall apply this formalism to the study of spin-orbit effects in the equations of motion and energy flux of compact binaries.

Hadamard regularization
The stress-energy tensor of point masses has been defined in the previous Section by means of Dirac functions, and involves metric coefficients evaluated at the locations of the point particles, namely (g αβ ) A . However, it is clear that the metric g αβ becomes singular at the particles. Indeed this is already true at Newtonian order where the potential generated by N particles reads where r B ≡ |x − y B | is the distance between the field point and the particle B. Thus the values of U and hence of the metric coefficients [recall that g 00 = −1 + 2U/c 2 + O(c −4 )], are ill-defined at the locations of the particles. What we need is a self-field regularization, i.e. a prescription for removing the infinite self-field of the point masses. Arguably the choice of a particular regularization constitutes a fully qualified element of our physical modelling of compact objects. At Newtonian order the regularization of the potential (80) should give the well-known result where r AB ≡ |y A − y B | and the infinite self interaction term has simply been discarded from the summation. At high post-Newtonian orders the problem is not trivial and the self-field regularization must be properly defined. The post-Newtonian formalism reviewed in Sections 2-7 assumed from the start a continuous (smooth) matter distribution. Actually this formalism will be applica-ble to singular point-mass sources, described by the stress-energy tensor of Section 8, provided that we suplement the scheme by the self-field regularization. Note that this regularization has nothing to do with the finite-part process FP B=0 extensively used in the case of extended matter sources. The latter finite part was an ingredient of the rigorous derivation of the general post-Newtonian solution [see (20)], while the self-field regularization is an assumption regarding a particular type of singular source.
Our aim is to compute up to 3PN order the metric coefficients at the location of one of the particles: (g αβ ) A . At this stage different self-field regularizations are possible. We first review Hadamard's regularization [66,67], which has proved to be very efficient for doing practical computations, but suffers from the important drawback of yielding some "ambiguity parameters", which cannot be determined within the regularization, starting at the 3PN order.
Iterating the Einstein field equations with point-like matter sources (delta functions with spatial supports localized on y A ) yields a generic form of functions representing the metric coefficients in successive post-Newtonian approximations. The generic functions, say F(x), are smooth except at the points y A , around which they admit singular Laurent expansions in powers and inverse powers of r A ≡ |x − y A |. When r A → 0 we have (say, for any P ∈ N) The coefficients A f p of the various powers of r A depend on the unit direction n A ≡ (x − y A )/r A of approach to the singular point A. The powers p are relative integers, and are bounded from below by p 0 ∈ Z. The Landau o-symbol for remainders takes its standard meaning. The A f p 's depend also on the (coordinate) time t, through their dependence on velocities v B (t) and relative positions y BC (t) ≡ y B (t) − y C (t); however the time t is purely "spectator" in the regularization process, and thus will not be indicated. The coefficients A f p for which p < 0 are referred to as the singular coefficients of F around A. The function F being given that way, we define the Hadamard partie finie as the following value of F at the location of the particle A, 10 where dΩ A denotes the solid angle element centred on y A and sustained by n A . The brackets mean the angular average. The second notion of Hadamard partie finie concerns the integral d 3 x F, which is generically divergent at the points y A . Its partie finie (in short Pf) is defined by The first term integrates over the domain R 3 \ N A=1 B A (s) defined as R 3 deprived from the N spherical balls B A (s) ≡ {x; r A ≤ s} of radius s and centred on the points y A . The second term is the opposite of the sum of divergent parts associated with the first term around each of the particles in the limit where s → 0. We have Since the divergent parts are cancelled (by definition) the Hadamard partie finie is obtained in the limit s → 0. Notice that as indicated in (84) the Hadamard partie-finie integral is not fully specified: it depends on N strictly positive and a priori arbitrary constants s 1 , · · · , s N parametrizing the logarithms in (85).
We have seen that the post-Newtonian scheme consists of breaking the hyperbolic d'Alembertian operator into the elliptic Laplacian ∆ and the retardation term c −2 ∂ 2 t considered to be small, and put in the right-hand-side of the equation where it can be iterated; see (18). We thus have to deal with the regularization of Poisson integrals, or iterated Poisson integrals, of some generic function F. The Poisson integral will be divergent 11 and we apply the prescription (84). Thus, This definition is valid for each field point x ′ different from the y A 's, and we want to investigate the singular limit when x ′ tends to one of the source points y A , so as to define the object (P) A . The definition (83) is not directly applicable because the expansion of the Poisson integral P(x ′ ) when x ′ → y A will involve besides the normal powers of r ′ A ≡ |x ′ − y A | some logarithms of r ′ A . The proper way to define the Hadamard partie finie in this case is to include the ln r ′ A into the definition (83) as if it were a mere constant parameter. With this definition we arrive at [68] The first term is given by a partie-finie integral following the definition (84); the second involves the logarithm of r ′ A . The constants s 1 , · · · , s N come from (86). Since r ′ A is actually tending to zero, ln r ′ A represents a formally infinite "constant", which will ultimately parametrize the final Hadamard regularized 3PN equations of motion. In the two-body case we shall find that the constants r ′ A are unphysical in the sense that they can be removed by a coordinate transformation [69]. Note that the apparent dependence of (87) on the constant s A is illusory. Indeed the dependence on s A cancels out between the first and the second terms in the right-hand-side of (87), so the result depends only on r ′ A and the s B 's for B = A. We thus have a simpler rewriting of (87) as Unfortunately, the constants s B for B = A remaining in the result (88) will be the source of a genuine ambiguity. This ambiguity can in fact be traced back to the socalled "non-distributivity" of the Hadamard partie finie, a consequence of the presence of the angular integration in (83), and implying that (FG) A = (F) A (G) A in general. The non-distributivity arises precisely at the 3PN order both in the equations of motion and radiation field of point-mass binaries. At that order we are loosing with Hadamard's regularization an elementary rule of ordinary calculus. Consequently we expect that some basic symmetries of general relativity such as diffeomorphism invariance will be lost. However Hadamard's regularization can still be efficiently used to compute most of the terms in the equations of motion and radiation field at 3PN order; only a few ambiguous terms will show up which have then to be determined by another method.

Dimensional regularization
Dimensional regularization is an extremely powerful regularization which is free of ambiguities (at least up to the 3PN order). The main reason is that it is able to preserve the symmetries of classical general relativity; in fact dimensional regularization was invented [70,71] as a means to preserve the gauge symmetry of perturbative quantum field theories. In the present context we shall show that dimensional regularization permits to resolve the problem of ambiguities arising at the 3PN order in Hadamard's regularization. We shall employ dimensional regularization not merely as a trick to compute some particular integrals which would otherwise be divergent or ambiguous, but as a fundamental tool for solving in a consistent way the Einstein field equations with singular sources. We therefore assume that the correct theory is general relativity in D = d + 1 space-time dimensions. As usual, any intermediate formulas will be interpreted by analytic continuation for a general complex spatial dimension d ∈ C. In particular we shall analytically continue d down to the value of interest 3 and pose d = 3 + ε .
The Einstein field equations in d spatial dimensions take the same form as presented in Section 2, with the exception that the explicit expression of the gravitational source term Λ αβ now depends on d. We find that only the last term in (11) acquires a dependence on d; namely the factor 1 2 in g ρσ g επ − 1 2 g σ ε g ρπ should now read 1 d−1 . In addition the d-dimensional gravitational constant is related to the usual three-dimensional Newton constant G by where ℓ 0 is a characteristic length associated with dimensional regularization.
In the post-Newtonian iteration performed in d dimensions we shall meet the analogue of the function F, which we denote by F (d) (x) where x ∈ R d . It turns out that in the vicinity of the singular points y A , the function F (d) admits an expansion richer than in (82), and of the type The coefficients A f (ε) p,q (n 1 ) depend on the dimension through ε ≡ d − 3 and also on the scale ℓ 0 . The powers of r A are now of the type p + qε where the two relative integers p, q ∈ Z have values limited as indicated. Because F (d) reduces to F when ε = 0 we necessarily have the constraints To proceed with the iteration we need the Green function of the Laplace operator in d dimensions. Its explicit form is It satisfies ∆ u where Γ is the usual Eulerian function. It reduces to one when d → 3. Note that the volume Ω d−1 of the sphere with d − 1 dimensions is related to K by With these results the Poisson integral of F (d) , constituting the d-dimensional analogue of (86), reads In dimensional regularization the singular behavior of this integral is automatically taken care of by analytic continuation in d. Next we evaluate the integral at the singular point x ′ = y A . In contrast with Hadamard's regularization where the result was given by (87), in dimensional regularization this is quite easy, as we are allowed to simply replace x ′ by y A into the explicit integral form (96). So we simply have The main step of our strategy [72,73] will now consist of computing the difference between the d-dimensional Poisson potential (97) and its 3-dimensional counterpart which is defined from Hadamard's regularization as (87). We shall then add this difference (in the limit ε = d − 3 → 0) to the result obtained by Hadamard regularization in order to get the corresponding dimensional regularization result. This strategy is motivated by the fact that as already mentionned most of the terms do not present any problems and have already been correctly computed using Hadamard's regularization. Denoting the difference between the two regularizations by means of the script letter D, we write We shall only compute the first two terms of the Laurent expansion of D(P) A when ε → 0, which will be of the form D(P) A = a −1 ε −1 + a 0 + O(ε). This is the information needed to determine the value of the ambiguity parameters. Notice that the difference D(P) A comes exclusively from the contribution of terms developing some poles ∝ 1/ε in the d-dimensional calculation. The ambiguity in Hadamard's regularization at 3PN order is reflected by the appearance of poles in d dimensions.
The point is that in order to obtain the difference D(P) A we do not need the expression of F (d) for an arbitrary source point x but only in the vicinity of the singular points y A . Thus this difference depends only on the singular coefficients of the local expansions of F (d) near the singularities. We find [73] We still use the bracket notation to denote the angular average but this time per- The above differences for all the Poisson and interated Poisson integrals composing the equations of motion (i.e. the accelerations of the point masses) are added to the corresponding results of the Hadamard regularization in the variant of it called the "pure Hadamard-Schwartz" regularization (see [73] for more details). In this way we find that the equations of motion in dimensional regularization are composed of a pole part ∝ 1/ε which is purely 3PN, followed by a finite part when ε → 0, plus the neglected terms O(ε). It has been shown (in the two-body case N = 2) that: 1. The pole part ∝ 1/ε of the accelerations can be renormalized into some shifts of the "bare" world-lines by y A → y A + ξ ξ ξ A , with ξ ξ ξ A containing the poles, so that the result expressed in terms of the "dressed" world-lines is finite when ε → 0; 2. The renormalized acceleration is physically equivalent to the result of the Hadamard regularization, in the sense that it differs from it only by the effects of shifts ξ ξ ξ A , if and only if the ambiguities in the Hadamard regularization are fully and uniquely determined (i.e. take specific values).
These results [73] provide an unambiguous determination of the equations of motion of compact binaries up to the 3PN order. A related strategy with similar complete results has been applied to the problem of multipole moments and radiation field of point-mass binaries [74]. This finally completed the derivation of the general relativistic prediction for compact binary inspiral up to 3PN order (and even to 3.5PN order). In later Sections we shall review some features of the 3.5PN gravitationalwave templates of inspiralling compact binaries. Why should the final results of the employed regularization scheme be unique, in agreement with our expectation that the problem is well-posed and should possess a unique physical answer? The results can be justified by invoking the "effacing principle" of general relativity [20] -namely that the internal structure of the compact bodies does not influence the equations of motion and emitted radiation until a very high post-Newtonian order. Only the masses m A of the bodies should drive the motion and radiation, and not for instance their "compactness" Gm A /(c 2 a A ). A model of point masses should therefore give the correct physical answer, which we expect to be also valid for black holes, provided that the regularization scheme is mathematically consistent.

Energy and flux of compact binaries
The equations of motion of compact binary sources, up to the highest known post-Newtonian order which is 3.5PN, will serve in the definition of the gravitationalwave templates for two purposes: 1. To compute the center-of-mass energy E appearing in the left-hand-side of the energy balance equation to be used for deducing the orbital phase, 2. To order-reduce the accelerations coming from the time derivatives of the source multipole moments required to compute the gravitational-wave energy flux F in the right-hand-side of the balance equation.
We consider two compact objects moving under purely gravitational mutual interaction. In a first stage we assume that the bodies are non-spinning so the motion takes place in a fixed plane, say the x-y plane. The relative position x = y 1 − y 2 , velocity v = dx/dt, and acceleration a = dv/dt are given by The orbital frequency ω is related in the usual way to the orbital phase φ by ω =φ (time derivatives are denoted with a dot). Here the vector λ λ λ =ẑ × n is perpendicular to the unit vectorẑ along the z-direction orthogonal to the orbital plane, and to the binary's separation unit direction n ≡ x/r. Through 3PN order, it is possible to model the binary's orbit as a quasi-circular orbit decaying by the effect of radiation reaction at the 2.5PN order. The restriction to quasi-circular orbits is both to simplify the presentation, 12 and for physical reasons because the orbit of inspiralling compact binaries detectable by current detectors should be circular (see the discussion in Section 1). The radiation-reaction effect at 2.5PN order yields 13 Gm where γ is defined as the small [i.e. γ = O(c −2 )] post-Newtonian parameter Substituting these results into (102), we obtain the expressions for the velocity and acceleration during the inspiral, v = r ω λ λ λ − 64 5 Gm Notice that whileṙ = O(c −5 ), we haver = O(c −10 ) which is of the order of the square of radiation-reaction effects and is thus zero with the present approximation.
A central result of post-Newtonian calculations is the expression of the orbital frequency ω in terms of the binary's separation r up to 3PN order. This result has been obtained independently by three groups. Two are working in harmonic coordinates: Blanchet & Faye [69,75,73] use a direct post-Newtonian iteration of the equations of motion, while Itoh & Futamase [76,77,78] apply a variant of the surface-integral approach (à la Einstein-Infeld-Hoffmann [5]) valid for compact bodies without the need of a self-field regularization. The group of Jaranowski & Schäfer [79,80,72] employs Arnowitt-Deser-Misner coordinates within the Hamiltonian formalism of general relativity. 14 The 3PN orbital frequency in harmonic coordinates is Note the logarithm at 3PN order coming from a Hadamard self-field regularization scheme, and depending on a constant r ′ 0 defined by m ln r ′ 0 ≡ m 1 ln r ′ 1 + m 2 ln r ′ 2 , where r ′ A ≡ |x ′ − y A | are arbitrary "constants" discussed in Section 9. We shall see that r ′ 0 disappears from final results -it can be qualified as a gauge constant. To obtain the 3PN energy we need to go back to the equations of motion for general non-circular orbits, and deduce the energy as the integral of the motion associated with a Lagrangian formulation of (the conservative part of) these equations [75]. Once we have the energy for general orbits we can reduce it to quasi-circular orbits. We find A convenient post-Newtonian parameter x = O(c −2 ) is now used in place of γ; it is defined from the orbital frequency as The interest in this parameter stems from its invariant meaning in a large class of coordinate systems including the harmonic and ADM coordinates. By inverting (106) we find at 3PN order Section 11, with the orbital phase φ measuring the direction of the unit separation vector n = x/r, then n =x cos φ +ŷ sin φ , wherex andŷ are the unit directions along x and y. Following [88,89] we choose the polarization vector P to lie along the x-axis and the observer to be in the y-z plane in the direction N =ŷ sin i +ẑ cos i , where i is the orbit's inclination angle. With this choice P lies along the intersection of the orbital plane with the plane of the sky in the direction of the ascending node, i.e. that point at which the bodies cross the plane of the sky moving toward the observer. Hence the orbital phase φ is the angle between the ascending node and the direction of body 1. The rotating orthonormal triad (n, λ λ λ ,ẑ) describing the motion of the binary and used in (102) is related to the fixed polarization triad (N, P, Q) by The 3.5PN expression of the orbital phase φ as function of the orbital frequency or equivalently the x-parameter is obtained from the energy balance equation (101) in which the binary's conservative center-of-mass energy E and total gravitationalwave flux F have been obtained in (110)-(111). For circular orbits the orbital phase is computed from Various methods (numerical or analytical) are possible for solving (115) given the expressions (110) and (111). This yields different waveform families all valid at the same 3.5PN order, but which may differ when extrapolated beyond the normal domain of validity of the post-Newtonian expansion, i.e. in this case very near the coalescence. Such differences must be taken into account when comparing the post-Newtonian waveforms to numerical results [52]. It is convenient to perform a change of phase, from the actual orbital phase φ to the new phase variable where M is the binary's total mass monopole moment 15 and ω 0 = 1 4u 0 exp[ 11 12 − C] is related to the constant u 0 ≡ r 0 /c entering the tail integrals in (60). The logarithmic term in ψ corresponds physically to some spreading of the different frequency components of the wave along the line of sight from the source to the detector, and expresses the tail effect as a small delay in the arrival time of gravitational waves. This effect, although of formal 1.5PN order in (116), represents in fact a very small modulation of the orbital phase: compared to the dominant phase evolution whose order is that of the inverse of 2.5PN radiation reaction, this modulation is of order 4PN and can thus be neglected with the present accuracy.
The spherical harmonic modes of the polarization waveforms can now be obtained at 3PN order using the angular integration formula (55). We start from the expressions of the wave polarizations h + and h × as functions of the inclination angle i and of the phase ψ. We use the known dependence of the spherical harmonics (54) on the azimuthal angle. Denoting h ≡ h + − ih × = h(i, ψ) we find that the latter angular integration becomes exhibiting the azimuthal factor e −i mψ appropriate for each mode. Let us introduce a normalized mode coefficient H ℓm starting by definition with one at the Newtonian order for the dominant mode having (ℓ, m) = (2, 2). This means posing All the modes have been given in [54] up to the 3PN order. The dominant mode (2, 2), which is primarily needed for the comparison between post-Newtonian calculations and numerical simulations, reads at 3PN order

Spin-orbit contributions in the energy and flux
To successfully detect the gravitational waves emitted by spinning, precessing binaries and to estimate the binary parameters, spin effects should be included in the templates. For maximally spinning compact bodies the spin-orbit coupling (linear in the spins) appears dominantly at the 1.5PN order, while the spin-spin one (which is quadratic) appears at 2PN order. The spin effect on the free motion of a test particle was obtained by Papapetrou [64] in the form of a coupling to curvature. Seminal later works by Barker & O'Connell [90,91] obtained the leading order spin-orbit and spin-spin contributions in the post-Newtonian equations of motion. Based on these works, the spin-orbit and spin-spin terms were obtained in the radiation field [92,65,93,94], enabling the derivation of the orbital phase evolution (the crucial quantity that determines the templates). Finding the 1PN corrections to the leading spin-orbit coupling in both the (translational) equations of motion and radiation field was begun in [95,96] and completed in [97,98]. The result [97] for the equations of motion was confirmed by an alternative derivation based on the ADM-Hamiltonian formalism [99]. In Section 8 we discussed a covariant formalism for spinning particles [59,60,61,62,63]. We want now to find a convenient three-dimensional variable for the spin. Restricting ourselves to spin-orbit effects, i.e. neglecting O(S 2 ), we can write the components of the spin tensor S αβ A as We have used the momentum-velocity relation (76) and have taken into account the spin condition (73). A first possibility is to adopt as the basic spin variable the contravariant components of the spin covector S A i in (120), which are obtained by raising the index by means of the contravariant spatial metric, viz where γ i j is the inverse of the covariant spatial metric γ i j ≡ g i j i.e. satisfies γ ik γ k j = δ i j . The choice of spin variable (121) has been adopted in [95,96]. However to express final results (to be used in gravitational-wave templates) it is better to use a different set of spin variables characterized by having some conserved Euclidean lengths. Such spins with constant Euclidean magnitude will be denoted by S c A . They can be computed in a straightforward way at a given post-Newtonian order in terms of the previous variables (121). For instance we find up to 1PN order, The (regularized) gravitational potential (U) A is defined by (81). The constantmagnitude spin variable S c A obeys a spin precession equation which is necessarily of the form dS c Indeed this equation implies that |S c A | = const. The precession angular frequency vector Ω Ω Ω A for two-body systems has been computed for the leading spin-orbit and spin-spin contributions [65] and for the 1PN correction to the spin-orbit [97,98,99].
For two bodies we conveniently use the following combinations (introduced in [65]) of the two spins: We recall that ∆ ≡ (m 1 − m 2 )/m; see the footnote 13. This expression is valid for quasi-circular orbits, and we neglect the spin-spin terms. Similarly the gravitationalwave flux will be modified at the same 1.5PN and 2.5PN orders. Posing F = F mono + F spin where F mono is given by (111), we find Having in hand the spin contributions in E and F , we can deduce the evolution of the orbital phase from the energy balance equation (101). In absence of precession of the orbital plane, e.g. for spins aligned or anti-aligned with the orbital angular momentum, the gravitational-wave phase will reduce to the "carrier" phase φ GW ≡ 2φ (keeping only the dominant harmonics), where φ is the orbital phase which is obtained by integrating the orbital frequency. However, in the general case of non-aligned spins, we must take into account the effect of precession of the orbital plane induced by spin modulations. Then the gravitational-wave phase is given by Φ GW = φ GW + δ φ GW , where the precessional correction δ φ GW arises from the changing orientation of the orbital plane, and can be computed by standard methods using numerical integration [100]. Thus, the carrier phase φ GW constitutes the main theoretical output to be provided for the gravitational-wave templates, and can directly be obtained numerically from using the integration formula (115).