Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries

Abstract

The article reviews the current status of a theoretical approach to the problem of the emission of gravitational waves by isolated systems in the context of general relativity. Part A of the article deals with general post-Newtonian 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 we obtain the explicit expressions of the source multipole moments. The relationships between the radiative and source moments involve many nonlinear multipole interactions, among them those associated with the tails (and tails-of-tails of gravitational waves. Part B of the article is devoted to the application to compact binary systems. 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 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.

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 French XVIIIth century 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 (denoted □), invented by d’Alembert in his Traité de dynamique of 1743.

Nowadays, the importance of the field lies in the exciting possibility of comparing the theory with contemporary astrophysical observations, made by a new generation of detectors — large-scale optical interferometers LIGO, VIRGO, GEO and TAMA — that should routinely observe the gravitational waves produced by massive and rapidly evolving systems such as inspiralling compact binaries. 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 waves from the source to the distant detector, as well as their back-reaction onto the source.

Gravitational-wave generation formalisms

The basic problem we face is to relate the asymptotic gravitational-wave form h_ij generated by some isolated source, at the location of some detector in the wave zone of the source, to the stress-energy tensor T^αβ of the matter fields¹. 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, keeping in mind that, sadly, such methods are often not related in a very precise mathematical way to the first principles of the theory. Therefore, a general wave-generation formalism must somehow manage the non-linearity of the field equations by imposing some suitable approximation series in one or several small physical parameters. Of ourse the ultimate aim of approximation methods is to extract from the theory some firm predictions for the outcome of experiments such as VIRGO and LIGO. 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), and the far-zone expansion (1/R-expansion in the distance). In particular, the post-Newtonian expansion has provided us in the past with our best insights into the problems of motion and radiation in general relativity. The most successful wave-generation formalisms make a gourmet cocktail of all these approximation methods. For reviews on analytic approximations and applications to the motion and the gravitational wave-generation see Refs. [143, 53, 54, 144, 150, 8, 13].

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 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 first to perform 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

$$ \in = \max \left\{ {\left| {\frac{{{T^{0i}}}}{{{T^{00}}}}} \right|,{{\left| {\frac{{{T^{ij}}}}{{{T^{00}}}}} \right|}^{1/2}},{{\left| {\frac{U}{{{c^2}}}} \right|}^{1/2}}} \right\},$$

((1))

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. [40, 42, 41], we shall henceforth write ε ≡ 1/c, even though ε is dimensionless whereas c has the dimension of a velocity. The small post-Newtonian remainders will be denoted $\hat L$. Thus, 1/c ≪ 1 in the case of post-Newtonian sources. We have |U/c²|^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. $\hat J$. 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 1/c can be as large as 30% in the last rotations, and whose description necessitates the control of high post-Newtonian approximations.

The lowest-order wave generation formalism, in the Newtonian limit 1/c → 0, is the famous quadrupole formalism of Einstein [68] and Landau and Lifchitz [97]. This formalism can also be referred to as Newtonian because the evolution of the quadrupole moment of the source is computed using Newton’s laws of gravity. It expresses the gravitational field h ^TT _ij in a transverse and traceless (TT) coordinate system, covering the far zone of the source², as

$$h_{ij}^{TT} = \frac{{2G}}{{{c^4}R}}{\mathcal{P}_{ijab}}(N)\left\{ {\frac{{{d^2}{Q_{ab}}}}{{d{T^2}}}(T - R/c) + \mathcal{O}\left( {\frac{1}{c}} \right)} \right\} + O\left( {\frac{1}{{{R^2}}}} \right),$$

((2))

where R = |X| is the distance to the source, N = X/R is the unit direction from the source to the observer, and $\frac{1}{{\sqrt {2\pi } }}$ is the TT projection operator, with $\hat {\mathcal{E}},\;\hat {\mathcal{C}}, \text{and} \hat {\mathcal{l}}_{z}$ being the projector onto the plane orthogonal to N. The source’s quadrupole moment takes the familiar Newtonian form

$${Q_{ij}}(t) = \int_{source} {d{}^3x\rho (x,t)\left( {{x_i}{x_j} - \frac{1}{3}{\delta _{ij}}{x^2}} \right)} ,$$

((3))

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

$$\mathcal{L} = \frac{G}{{5{c^5}}}\left\{ {\frac{{{d^3}{Q_{ab}}}}{{d{T^3}}}\frac{{{d^3}{Q_{ab}}}}{{d{T^3}}} = \mathcal{O}\left( {\frac{1}{{{c^2}}}} \right)} \right\}.$$

((4))

Our notation ${\mathcal{E}}, \text{and}\; \hat {\mathcal{l}}_{z}$ stands for the total gravitational “luminosity” of the source. The cardinal virtues of the Einstein-Landau-Lifchitz 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 Hulse-Taylor binary pulsar PSR 1913+16 [140, 141, 139]. Indeed the prediction of the quadrupole formalism for the waves emitted by the binary pulsar system comes from applying Eq. (4) to a system of two point masses moving on an eccentric orbit (the classic reference is Peters and Mathews [117]; see also Refs. [71, 148]). Then, relying on the energy equation

$$\frac{{dE}}{{dt}} = - \mathcal{L},$$

((5))

where E is the Newtonian binary’s center-of-mass energy, we deduce from Kepler’s third law the expression of the “observable”, that is, the change in the orbital period P of the pulsar, or P, as a function of P itself. From the binary pulsar test, we can say that the post-Newtonian corrections to the quadrupole formalism, which we shall compute in this article, have already received, in the case of compact binaries, strong observational support (in addition to having, as we shall demonstrate, a sound theoretical basis).

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, 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 formula (4) to include the current-quadrupole and mass-octupole moments [111, 110], and obtained the corresponding formulas for linear momentum [111, 110, 1, 124] and angular momentum [116, 46]. The general structure of the infinite multipole series in the linearized theory was investigated by several works [126, 127, 119, 142], 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 [142]. The explicit expressions of the multipole moments (for instance in STF guise) as integrals over the source, valid in the linearized theory but irrespective of a slow motion hypothesis, are completely known [101, 39, 38, 57].

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. [32] and Sachs [128]. The precise way under which such radiative space-times fall off asymptotically has been formulated geometrically by Penrose [114, 115] in the concept of an asymptotically simple space-time (see also Ref. [76]). 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 behaviour 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 and Damour [14, 3], following pioneering work by Bonnor and collaborators [33, 34, 35, 81] and Thorne [142]. 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 (whereas, the direct attack of the post-Minkowskian expansion, valid at once inside and outside the source, faces some calculational difficulties [147, 47]). In this “multipolar-post-Minkowskian” 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 [6, 11]. 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 [37, 36].

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 [14]. 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 farzone expansions of that general solution will be stated. Notably, we find [3] that the asymptotic behaviour of the solution at future null infinity is in agreement with the findings of the Bondi-Sachs-Penrose [32, 128, 114, 115, 76] 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 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? Only in the case of post-Newtonian sources has it been possible to answer this question. The general expression of the moments was obtained at the level of the second post-Newtonian (2PN) order in Ref. [6], and was subsequently proved to be in fact valid up to any post-Newtonian order in Ref. [11]. 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 ($\left( { - \frac{{2399}}{{56}}q\; - \;\frac{{773}}{3}\pi } \right)\;{x^5}$ 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. [16] for the mass-type moments (strangely enough, the mass moments admit a compact-support expression at 1PN order), and in Ref. [58] for the current-type ones.

The wave-generation formalism resulting from matching the exterior multipolar and post-Minkowskian field [14, 3] to the post-Newtonian source [6, 11] 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 3PN order in Refs. [17, 12, 10]. 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 [15, 5, 9]. 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, which means 4PN “absolute” order [15].

A different wave-generation formalism has been devised by Will and Wiseman [152] (see also Refs. [151, 112]), after earlier attempts by Epstein and Wagoner [70] and Thorne [142]. This formalism has exactly the same scope as ours, 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 [152]. 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 $\left( { - \frac{{1041349}}{{18144}} + \frac{{171}}{{16}}{q^2} - \frac{{243}}{8}q - \frac{{785}}{6}\pi } \right){x^5}$ enclosing the source, e.g., the radius of the near zone). By contrast, in our case [6, 11], the moments are given by some integrals covering the whole space and regularized by means of the finite part $y = \frac{C}{{{Q^2}}},\;{Q^2} = l_z^2 + {a^2}(1 - {E^2})$. We shall prove the complete equivalence, at the most general level, between the two formalisms. What is interesting about both formalisms is that the source multipole moments, which involve a whole series of relativistic corrections, are coupled together, in the true non-linear solution, in a very complicated way. These multipole couplings give rise to the many tail and related non-linear effects, which form an integral part of the radiative moments at infinity and thereby of the observed signal.

Part A of this article is devoted to a presentation of the post-Newtonian wave generation formalism. We try to state the main results in a form that is simple enough to be understood without the full details, but at the same time we 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.

Problem posed by compact binary systems

Inspiralling compact binaries, containing neutron stars and/or black holes, are promising sources of gravitational waves detectable by the detectors LIGO, VIRGO, GEO and TAMA. 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. For instance, the eccentricity of the orbit of the Hulse-Taylor binary pulsar is presently e₀ = 0.617. At the time when the gravitational waves emitted by the binary system will become visible by the detectors, i.e. when the signal frequency reaches about 10 Hz (in a few hundred million years from now), the eccentricity will be e = 5.3 × 10^-6 — a value calculated from the Peters [116] law, which is itself based on the quadrupole formula (2).

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 m₁ and m₂ (and possibly two spins), is sufficient. 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, and so on, are dominated by gravitational effects. However, the real 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 Q₁ and Q₂ induced by tidal interaction between two neutron stars. Let a₁ and a₂ be the radius of the stars, and L the distance between the two centers of mass. We have, for tidal moments,

$${Q_1} = {k_1}{m_2}\frac{{a_1^5}}{{{L^3}}},\;\;\;\;{Q_2} = {k_2}{m_1}\frac{{a_2^5}}{{{L^3}}},$$

((6))

where k₁ and k₂ are the star’s dimensionless (second) Love numbers [103], 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”, defined by the dimensionless ratios

$${K_1} = \frac{{G{m_1}}}{{{a_1}{c^2}}},\;\;\;\;{K_2} = \frac{{G{m_2}}}{{{a_2}{c^2}}},$$

((7))

which equal ∼ 0.2 for neutron stars (depending on their equation of state). The quadrupoles Q₁ and Q₂ will affect both sides of Eq. (5), i.e. the Newtonian binding energy E of the two bodies, and the emitted total gravitational flux $\mathcal{O}(\epsilon{^2})$ 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 [92] 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 we find that within such a model (in the case of two identical neutron stars) the orbital phase, deduced from Eq. (5), reads

$${\phi ^{finite\;size}} - {\phi _0} = - \frac{1}{{8{x^{5/2}}}}\left\{ {1 + const.\;k{{\left( {\frac{x}{K}} \right)}^5}} \right\},$$

((8))

where x — (Gmω/c³)^2/3 is a standard dimensionless post-Newtonian parameter ∼ 1/c² (ω is the orbital frequency), and where k is the Love number and K is the compactness of the neutron star. The first term in the right-hand side of (8) corresponds to the gravitational-wave damping of two point masses; the second term is the finite-size effect, which appears as a relative correction, proportional to (x/K)⁵, to the latter radiation damping effect. Because the finite-size effect is purely Newtonian, its relative correction ∼ (x/K)⁵ should not depend on c; and indeed the factors 1/c² cancel out in the ratio x/K. However, the compactness K of compact objects is by Eq. (7) of the order unity (or, say, one half), therefore the 1/c² it contains should not be taken into account numerically in this case, and so the real order of magnitude of the relative contribution of the finite-size effect in Eq. (8) is given by x⁵ alone. This means that for compact objects the finite-size effect should be comparable, numerically, to a post-Newtonian correction of order 5PN or 1/c¹⁰ (see Ref. [52] for the proof in the context of relativistic equations of motion). 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′ ≡ const. k ∼ 1 and K ∼ 0.2 for neutron stars (and the bandwidth of a VIRGO detector 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 point masses. The conclusion is that the finite-size effect can in general be neglected in comparison with purely gravitational-wave damping effects. 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 inspiralling compact binaries are ideally suited for application of a high-order post-Newtonian wave generation formalism. The main reason is that these systems are very relativistic, with orbital velocities as high as 0.3c in the last rotations (as compared to ∼ 10^-3c for the binary pulsar), and it is not surprising that the quadrupole-moment formalism (2, 3, 4, 5) constitutes a poor description of the emitted gravitational waves, since many post-Newtonian corrections play a substantial role. This expectation has been confirmed in recent years by several measurement-analyses [48, 49, 72, 50, 135, 121, 122, 96, 59], which have demonstrated that the post-Newtonian precision needed to implement successively the optimal filtering technique in the LIGO/VIRGO detectors corresponds grossly, in the case of neutron-star binaries, to the 3PN approximation, or 1/c⁶ 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 (∼ 16, 000 in the case of neutron stars), 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 the energy equation (5), on the center-of-mass binding energy E and the total gravitational-wave energy flux $\mathcal{O}\left( {{{(v/c)}^n}} \right)$.

In summary, the theoretical problem posed by inspiralling compact binaries is two-fold: On the one hand E, and on the other hand ${\psi _4} = - {C_{\alpha \beta \gamma \delta }}{n^\alpha }{\bar m^\beta }{n^\gamma }{\bar m^\delta }$, are to be deduced from general relativity with the 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 $_s\mathcal{O}$ it necessitates the application of a 3PN wave generation formalism (actually, things are more complicated because the equations of motion are also needed during the computation of the flux). It is quite interesting that such a high order approximation as the 3PN one should be needed in preparation for LIGO and VIRGO data analysis. As we shall see, the signal from compact binaries contains at the 3PN order the signature of several non-linear effects which are specific to general relativity. Therefore, we have here the possibility of probing, experimentally, some aspects of the non-linear structure of Einstein’s theory [28, 29].

Post-Newtonian equations of motion and radiation

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 term, at the 1PN order, was derived by Lorentz and Droste [98]. Subsequently, Einstein, Infeld and Hoffmann [69] obtained the 1PN corrections by means of their famous “surface-integral” method, in which the equations of motion are deduced from the vacuum field equations, and which 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 extended bodies, by Petrova [118] and Fock [73] (see also Ref. [109]).

The 2PN approximation was tackled by Otha et al. [105, 107, 106], who considered the post-Newtonian iteration of the Hamiltonian of N point-particles. We refer here to the Hamiltonian as the Fokker-type Hamiltonian, which is obtained from the matter-plus-field Arnowitt-Deser-Misner (ADM) Hamiltonian by eliminating the field degrees of freedom. The result for the 2PN and even 2.5PN equations of binary motion in harmonic coordinates was obtained by Damour and Deruelle [56, 55, 67, 51, 52], building on a non-linear iteration of the metric of two particles initiated in Ref. [2]. The corresponding result for the ADM-Hamiltonian of two particles at the 2PN order was given in Ref. [63] (see also Refs. [130, 131]). Kopeikin [93] derived the 2.5PN equations of motion for two extended compact objects. The 2.5PN-accurate harmonic-coordinate equations as well as the complete gravitational field (namely the metric g_αβ) generated by two point masses were computed by Blanchet, Faye and Ponsot [25], following a method based on previous work on wave generation [6].

Up to the 2PN level the equations of motion are conservative. Only at the 2.5PN order appears the first non-conservative effect, associated with the gravitational radiation reaction. The (harmonic-coordinate) equations of motion up to that level, as derived by Damour and Deruelle [56, 55, 67, 51, 52], have been used for the study of the radiation damping of the binary pulsar — its orbital P [52]. It is important to realize that the 2.5PN equations of motion have been proved to hold in the case of binary systems of strongly self-gravitating bodies [52]. This is via an “effacing” principle (in the terminology of Damour [52]) for the internal structure of the bodies. As a result, the equations depend only on the “Schwarzschild” masses, m₁ and m₂, of the compact objects. Notably their compactness parameters K₁ and K₂, defined by Eq. (7), do not enter the equations of motion, as has been explicitly verified up to the 2.5PN order by Kopeikin [93] and Grishchuk and Kopeikin [79], who made a “physical” computation, à la Fock, taking into account the internal structure of two self-gravitating extended bodies. The 2.5PN equations of motion have also been established by Itoh, Futamase and Asada [83, 84], who use a variant of the surface-integral approach of Einstein, Infeld and Hoffmann [69], that is valid for compact bodies, independently of the strength of the internal gravity.

The present state of the art is the 3PN approximation³. To this order the equations have been worked out independently by two groups, by means of different methods, and with equivalent results. On the one hand, Jaranowski and Schäfer [87, 88, 89], and Damour, Jaranowski and Schäfer [60, 62, 61], following the line of research of Refs. [105, 107, 106, 63], employ the ADM-Hamiltonian formalism of general relativity; on the other hand, Blanchet and Faye [21, 22, 20, 23], and de Andrade, Blanchet and Faye [66], founding their approach on the post-Newtonian iteration initiated in Ref. [25], compute directly the equations of motion (instead of a Hamiltonian) in harmonic coordinates. The end results have been shown [62, 66] to be physically equivalent in the sense that there exists a unique “contact” transformation of the dynamical variables that changes the harmonic-coordinates Lagrangian obtained in Ref. [66] into a new Lagrangian, whose Legendre transform coincides exactly with the Hamiltonian given in Ref. [60]. The 3PN equations of motion, however, depend on one unspecified numerical coefficient, ω_static, in the ADM-Hamiltonian formalism and Λ in the harmonic-coordinates approach, which is due to some incompleteness of the Hadamard self-field regularization method. This coefficient has been fixed by means of a dimensional regularization in Ref. [61].

So far the status of the post-Newtonian equations of motion is quite satisfying. There is mutual agreement between all the results obtained by means of different approaches and techniques, whenever it is possible to compare them: point particles described by Dirac delta-functions, extended post-Newtonian fluids, surface-integrals methods, mixed post-Minkowskian and post-Newtonian expansions, direct post-Newtonian iteration and matching, harmonic coordinates versus ADM-type coordinates, and different processes or variants of the regularization of the self field of point particles. In Part B of this article, we shall present the most complete results for the 3PN equations of motion, and for the associated Lagrangian and Hamiltonian formulations (from which we deduce the center-of-mass energy E).

The second sub-problem, that of the computation of the energy flux $\hat T$, has been carried out by application of the wave-generation formalism described previously. Following earliest computations at the 1PN level [149, 30], at a time when the post-Newtonian corrections in ${T_{nn}},{T_{\bar mn}},\;\text{and}\;{T_{\overline {mm} }}$ had a purely academic interest, the energy flux of inspiralling compact binaries was completed to the 2PN order by Blanchet, Damour and Iyer [18, 77], and, independently, by Will and Wiseman [152], using their own formalism (see Refs. [19, 27] for joint reports of these calculations). The preceding approximation, 1.5PN, which represents in fact the dominant contribution of tails in the wave zone, had been obtained in Refs. [153, 31] by application of the formula for tail integrals given in Ref. [17]. 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 by Blanchet [7, 10]. However, unlike the 1.5PN, 2.5PN and 3.5PN orders that are entirely composed of tail terms, the 3PN approximation also involves, besides the tails of tails, many nontail contributions coming from the relativistic corrections in the (source) multipole moments of the binary. These have been “almost” completed by Blanchet, Iyer and Joguet [26, 24], in the sense that the result still involves one unknown numerical coefficient, due to the use of the Hadamard regularization, which is a combination of the parameter λ in the equations of motion, and a new parameter θ coming from the computation of the 3PN quadrupole moment. In PartB of this article, we shall present the most up-to-date results for the 3.5PN energy flux and orbital phase, deduced from the energy equation (5), supposed to be valid at this order.

The post-Newtonian flux ${r_ \pm } = M \pm \sqrt {{M^2} - {a^2}}$, which comes from a “standard” post-Newtonian calculation, is in complete agreement (up to the 3.5PN order) with the result given by the very different technique of linear black-hole perturbations, valid in the “test-mass” limit where the mass of one of the bodies tends to zero (limit v → 0, where v = μ/m). Linear black-hole perturbations, triggered by the geodesic motion of a small mass around the black hole, have been applied to this problem by Poisson [120] at the 1.5PN order (following the pioneering work of Galt’sov et al. [75]), and by Tagoshi and Nakamura [135], using a numerical code, up to the 4PN order. This technique has culminated with the beautiful analytical methods of Sasaki, Tagoshi and Tanaka [129, 137, 138] (see also Ref. [102]), who solved the problem up to the extremely high 5.5PN order.

Part A: Post-Newtonian Sources

Einstein’s Field Equations

The field equations of general relativity form a system of ten second-order partial differential equations obeyed by the space-time metric g_αβ,

$${G^{\alpha \beta }}[g,\partial g,{\partial ^2}g] = \frac{{8\pi G}}{{{c^4}}}{T^{\alpha \beta }}[g],$$

((9))

where the Einstein curvature tensor $\bar R_{\ell m\omega }^{\text{in},\;\text{up}}$ is generated, through the gravitational coupling κ = 8πG/c⁴, by the matter stress-energy tensor T^αβ. Among these ten equations, four govern, via the contracted Bianchi identity, the evolution of the matter system,

$${\nabla _\mu }{G^{\alpha \mu }} \equiv 0\;\; \Rightarrow \;{\nabla _\mu }{T^{\alpha \mu }} = 0.$$

((10))

The space-time geometry is constrained by the six remaining equations, which place six independent constraints on the ten components of the metric g_αβ, leaving four of them to be fixed by a choice of a coordinate system.

In most of this paper we adopt the conditions of harmonic, or de Donder, coordinates. We define, as a basic variable, the gravitational-field amplitude

$${h^{\alpha \beta }} = \sqrt { - g} {g^{\alpha \beta }} - {\eta ^{\alpha \beta }},$$

((11))

where g^αβ denotes the contravariant metric (satisfying g^αμ g_μβ = δ ^α _β ), where g is the determinant of the covariant metric, g = det(g_αβ), and where η^αβ represents an auxiliary Minkowskian metric. The harmonic-coordinate condition, which accounts exactly for the four equations (10) corresponding to the conservation of the matter tensor, reads

$${\partial _\mu }{h^{\alpha \mu }} = 0.$$

((12))

The equations (11, 12) introduce into the definition of our coordinate system a preferred Minkowskian structure, with Minkowski metric η_αβ. Of course, this is not contrary to the spirit of general relativity, where there is only one physical metric g_αβ without any flat prior geometry, because the coordinates are not governed by geometry (so to speak), but rather are chosen by researchers when studying physical phenomena and doing experiments. Actually, the coordinate condition (12) is especially useful when we view the gravitational waves as perturbations of space-time propagating on the fixed Minkowskian manifold with the background metric η_αβ. This view is perfectly legitimate and represents a fructuous and rigorous way to think of the problem when using approximation methods. Indeed, the metric η_αβ, originally introduced in the coordinate condition (12), does exist at any finite order of approximation (neglecting higher-order terms), and plays in a sense the role of some “prior” flat geometry.

The Einstein field equations in harmonic coordinates can be written in the form of inhomogeneous flat d’Alembertian equations,

$$\square {h^{\alpha \beta }} = \frac{{16\pi G}}{{{c^4}}}{\tau ^{\alpha \beta }},$$

((13))

where □ ≡ □_η = η^μν∂_μ∂_ν. The source term, τ^αβ, can rightly be interpreted as the stress-energy pseudo-tensor (actually, τ^αβ is a Lorentz tensor) of the matter fields, described by T^αβ, and the gravitational field, given by the gravitational source term Λ^αβ, i.e.

$${\tau ^{\alpha \beta }} = |g|{T^{\alpha \beta }} + \frac{{{c^4}}}{{16\pi G}}{\Lambda ^{\alpha \beta }}.$$

((14))

The exact expression of Λ^αβ, including all non-linearities, reads

$$\begin{gathered} {\Lambda ^{\alpha \beta }} = - {h^{\mu v}}\partial _{\mu v}^2{h^{\alpha \beta }} + {\partial _\mu }{h^{\alpha v}}{\partial _v}{h^{\beta \mu }} + \frac{1}{2}{g^{\alpha \beta }}{g_{\mu v}}{\partial _\lambda }{h^{\mu \tau }}{\partial _\tau }{h^{v\lambda }} \hfill \\ \;\;\;\;\;\;\;\;\; - {g^{\alpha \mu }}{g_{v\tau }}{\partial _\lambda }{h^{\beta \tau }}{\partial _\mu }h{}^{v\lambda } - {g^{\beta \mu }}{g_{v\tau }}{\partial _\lambda }{h^{\alpha \tau }}{\partial _\mu }{h^{v\lambda }} + {g_{\mu v}}{g^{\lambda \tau }}{\partial _\lambda }{h^{\alpha \mu }}{\partial _\tau }{h^{\beta v}} \hfill \\ \;\;\;\;\;\;\;\;\; + \frac{1}{8}(2{g^{\alpha \mu }}{g^{\beta v}} - {g^{\alpha \beta }}{g^{\mu v}})(2{g_{\lambda \tau }}{g_{ \in \pi }} - {g_{\tau \in }}{g_{\lambda \pi }}){\partial _\mu }h{}^{\lambda \pi }{\partial _v}{h^{\tau \in }}. \hfill \\ \end{gathered} $$

((15))

As is clear from this expression, Λ^αβ is made of terms at least quadratic in the gravitational-field strength h and its first and second space-time derivatives. In the following, for the highest post-Newtonian order that we consider (3PN), we need the quadratic, cubic and quartic pieces of Λ^αβ. With obvious notation, we can write them as

$${\Lambda ^{\alpha \beta }} = {N^{\alpha \beta }}[h,h] + {M^{\alpha \beta }}[h,h,h] + {L^{\alpha \beta }}[h,h,h,h] + \mathcal{O}({h^5}).$$

((16))

These various terms can be straightforwardly computed from Eq. (15); see Eqs. (3.8) in Ref. [22] for explicit expressions.

As said above, the condition (12) is equivalent to the matter equations of motion, in the sense of the conservation of the total pseudo-tensor τ^αβ,

$${\partial _\mu }{\tau ^{\alpha \mu }} = 0\;\;\; \Leftrightarrow \;\;{\nabla _\mu }{T^{\alpha \mu }} = 0.$$

((17))

In this article, we look for the solutions of the field equations (13, 14, 15, 17) under the following four hypotheses:

• 1. The matter stress-energy tensor T^αβ is of spatially compact support, i.e. can be enclosed into some time-like world tube, say r ≤ a, where r = |x| is the harmonic-coordinate radial distance. Outside the domain of the source, when r > a, the gravitational source term, according to Eq. (17), is divergence-free,

  $${\partial _\mu }{\Lambda ^{\alpha \mu }} = 0\;\;\;(when\;r > a).$$

  ((18))

• 2. The matter distribution inside the source is smooth⁴: $\mathcal{E} = \mu \hat {\mathcal{E}},{l_z} = \mu {\hat l_z}$. We have in mind a smooth hydrodynamical “fluid” system, without any singularities nor shocks (a priori), that is described by some Eulerian equations including high relativistic corrections. In particular, we exclude from the start any black holes (however we shall return to this question when we find a model for describing compact objects).

• 3. The source is post-Newtonian in the sense of the existence of the small parameter defined by Eq. (1). For such a source we assume the legitimacy of the method of matched asymptotic expansions for identifying the inner post-Newtonian field and the outer multipolar decomposition in the source’s exterior near zone.

• 4. The gravitational field has been independent of time (stationary) in some remote past, i.e. before some finite instant $C = {\mu ^2}\hat C$ in the past, in the sense that

  $$\frac{\partial }{{\partial t}}[{h^{\alpha \beta }}(x,t)] = 0\;\;\;when\;t \leqslant - \mathcal{T}.$$

  ((19))

The latter condition is a means to impose, by brute force, the famous no-incoming radiation condition, ensuring that the matter source is isolated from the rest of the Universe and does not receive any radiation from infinity. Ideally, the no-incoming radiation condition should be imposed at past null infinity. We shall later argue (see Section 6) that our condition of stationarity in the past (19), although much weaker than the real no-incoming radiation condition, does not entail any physical restriction on the general validity of the formulas we derive.

Subject to the condition (19), the Einstein differential field equations (13) can be written equivalently into the form of the integro-differential equations

$${h^{\alpha \beta }} = \frac{{16\pi G}}{{{c^4}}}\square _{ret}^{ - 1}{\tau ^{\alpha \beta }},$$

((20))

containing the usual retarded inverse d’Alembertian operator, given by

$$(\square _{ret}^{ - 1}f)(x,t) \equiv - \frac{1}{{4\pi }}\smallint \smallint \smallint \frac{{{d^3}x'}}{{|x - x'|}}f(x',t - |x - x'|/c),$$

((21))

extending over the whole three-dimensional space $\dot t = dt/d\tau$.

Linearized Vacuum Equations

In what follows we solve the field equations (12, 13), in the vacuum region outside the compact-support source, in the form of a formal non-linearity or post-Minkowskian expansion, considering the field variable h^αβ as a non-linear metric perturbation of Minkowski space-time. At the linearized level (or first-post-Minkowskian approximation), we write:

$$h_{ext}^{\alpha \beta } = Gh_1^{\alpha \beta } + \mathcal{O}({G^2}),$$

((22))

where the subscript “ext” reminds us that the solution is valid only in the exterior of the source, and where we have introduced Newton’s constant G as a book-keeping parameter, enabling one to label very conveniently the successive post-Minkowskian approximations. Since h^αβ is a dimensionless variable, with our convention the linear coefficient h ^αβ₁ in Eq. (22) has the dimension of the inverse of G — a mass squared in a system of units where ħ = c = 1. In vacuum, the harmonic-coordinate metric coefficient h₁^αβ satisfies

$$\square h_1^{\alpha \beta } = 0,$$

((23))

$${\partial _\mu }h_1^{\alpha \mu } = 0.$$

((24))

We want to solve those equations by means of an infinite multipolar series valid outside a time-like world tube containing the source. Indeed the multipole expansion is the correct method for describing the physics of the source as seen from its exterior (r > a). On the other hand, the post-Minkowskian series is physically valid in the weak-field region, which surely includes the exterior of any source, starting at a sufficiently large distance. For post-Newtonian sources the exterior weak-field region, where both multipole and post-Minkowskian expansions are valid, simply coincides with the exterior r > a. It is therefore quite natural, and even, one would say inescapable when considering general sources, to combine the post-Minkowskian approximation with the multipole decomposition. This is the original idea of the “double-expansion” series of Bonnor [33], which combines the G-expansion (or m-expansion in his notation) with the a-expansion (equivalent to the multipole expansion, since the lth order multipole moment scales like a^l with the source radius).

The multipolar-post-Minkowskian method will be implemented systematically, using STF-harmonics to describe the multipole expansion [142], and looking for a definite algorithm for the approximation scheme [14]. The solution of the system of equations (23, 24) takes the form of a series of retarded multipolar waves⁵

$$h_1^{\alpha \beta } = \sum\limits_{l = 0}^{ + \infty } {{\partial _L}\left( {\frac{{K_L^{\alpha \beta }(t - r/c)}}{r}} \right),} $$

((25))

where r = |x|, and where the functions ${c_0} = (\ell - 1)\ell (\ell + 1)(\ell + 2) - 12iM\omega$ are smooth functions of the retarded time $R_{\ell m\omega }^{\text{in}}$, which become constant in the past, when $X_{\ell m\omega }^{\text{in}}$. It is evident, since a monopolar wave satisfies □(K_L(u)/r) = 0 and the d’Alembertian commutes with the multi-derivative ∂_L, that Eq. (25) represents the most general solution of the wave equation (23) (see Section 2 in Ref. [14] for a proof based on the Euler-Poisson-Darboux equation). The gauge condition (24), however, is not fulfilled in general, and to satisfy it we must algebraically decompose the set of functions K ⁰⁰ _L , K ^{0 i} _L , K ^ij _L into ten tensors which are STF with respect to all their indices, including the spatial indices i, ij. Imposing the condition (24) reduces the number of independent tensors to six, and we find that the solution takes an especially simple “canonical” form, parametrized by only two moments, plus some arbitrary linearized gauge transformation [142, 14].

Theorem 1 The most general solution of the linearized field equations (23, 24), outside some time-like world tube enclosing the source (r > a), and stationary in the past (see Eq. (19)), reads

$$h_1^{\alpha \beta } = k_1^{\alpha \beta } + {\partial ^\alpha }\varphi _1^\beta + {\partial ^\beta }\varphi _1^\alpha - {\eta ^{\alpha \beta }}{\partial _\mu }\varphi _1^\mu .$$

((26))

The first term depends on two STF-tensorial multipole moments, I_L (u) and J_L (u), which are arbitrary functions of time except for the laws of conservation of the monopole: I = const., and dipoles: I_i = const., J_i = const.. It is given by

$$\begin{gathered} k_1^{00} = - \frac{4}{{{c^2}}}\sum\limits_{l \geqslant 0} {\frac{{( - )l}}{{l!}}{\partial _L}} \left( {\frac{1}{r}{I_L}(u)} \right), \hfill \\ k_1^{0i} = - \frac{4}{{{c^3}}}\sum\limits_{l \geqslant 1} {\frac{{( - )l}}{{l!}}} \left\{ {{\partial _L}_{ - 1}\left( {\frac{1}{r}I_{iL - 1}^{(1)}(u)} \right) + \frac{l}{{l + 1}}{\varepsilon _{iab}}{\partial _{aL - 1}}\left( {\frac{1}{r}J{}_{bL - 1}(u)} \right)} \right\}, \hfill \\ k_1^{ij} = - \frac{4}{{{c^4}}}\sum\limits_{l \geqslant 2} {\frac{{( - )l}}{{l!}}} \left\{ {{\partial _L}_{ - 2}\left( {\frac{1}{r}I_{ijL - 2}^{(2)}(u)} \right) + \frac{{2l}}{{l + 1}}{\partial _{aL - 2}}\left( {\frac{1}{r}{\varepsilon _{ab}}({}_i{J_j})_{bL - 2}^{(1)}(u)} \right)} \right\}. \hfill \\ \end{gathered} $$

((27))

The other terms represent a linearized gauge transformation, with gauge vector ϕ ^α₁ of the type (25), and parametrized for four other multipole moments, say W_L(u), X_L(u), Y_L(u) and Z_L(u).

The conservation of the lowest-order moments gives the constancy of the total mass of the source, M ≡ I = const., center-of-mass position⁶, X_i ≡ I_i/I = const., total linear momentum P_i ≡ I ⁽¹⁾ _i = 0, and total angular momentum, S_i ≡ J_i = const.. It is always possible to achieve X_i = 0 by translating the origin of our coordinates to the center of mass. The total mass M is the Arnowitt-Deser-Misner (ADM) mass of the Hamiltonian formulation of general relativity. Note that the quantities M, X_i, P_i and S_i include the contributions due to the waves emitted by the source. They describe the “initial” state of the source, before the emission of gravitational radiation.

The multipole functions I_L(u) and J_L(u), which thoroughly encode the physical properties of the source at the linearized level (because the other moments W_L, ..., Z_L parametrize a gauge transformation), will be referred to as the mass-type and current-type source multipole moments. Beware, however, that at this stage the moments are not specified in terms of the stress-energy tensor T^αβ of the source: the above theorem follows merely from the algebraic and differential properties of the vacuum equations outside the source.

For completeness, let us give the components of the gauge-vector ϕ ^α₁ entering Eq. (26):

$$\begin{gathered} \varphi _1^0 = \frac{4}{{{c^3}}}\sum\limits_{l \geqslant 0} {\frac{{( - )l}}{{l!}}{\partial _L}} \left( {\frac{1}{r}{W_L}(u)} \right), \hfill \\ \varphi _1^i = - \frac{4}{{{c^4}}}\sum\limits_{l \geqslant 0} {\frac{{( - )l}}{{l!}}} {\partial _{iL}}\left( {\frac{1}{r}{X_L}(u)} \right) \hfill \\ \;\;\;\;\;\;\; - \frac{4}{{{c^4}}}\sum\limits_{l \geqslant 1} {\frac{{( - )l}}{{l!}}} \left\{ {{\partial _L}_{ - 1}\left( {\frac{1}{r}{Y_{iL - 1}}(u)} \right) + \frac{l}{{l + 1}}{\varepsilon _{iab}}{\partial _{aL - 1}}\left( {\frac{1}{r}{Z_{bL - 1}}(u)} \right)} \right\}. \hfill \\ \end{gathered} $$

((28))

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, in the following sense. If one takes these moments equal to zero and continues the calculations one ends up with a metric depending on I_L and J_L only, but that metric will not describe the same physical source as the one constructed from the six moments I_L, ..., Z_L. In other words, the two non-linear metrics associated with the sets of multipole moments {I_L, J_L, 0, ..., 0} and {I_L, J_L, W_L, ..., Z_L} are not isometric. We point out in Section 4.2 below that the full set of moments {I_L, J_L, W_L, ..., Z_L} is in fact physically equivalent to some reduced set {M_L, S_L, 0, ..., 0}, but with some moments M_L, S_L that differ from I_L, J_L by non-linear corrections (see Eq. (90)). All the multipole moments I_L, J_L, W_L, X_L, Y_L, Z_L will be computed in Section 5.

Non-linear Iteration of the Field Equations

By Theorem 1 we know the most general solution of the linearized equations in the exterior of the source. We then tackle the problem of the post-Minkowskian iteration of that solution. We consider the full post-Minkowskian series

$$h_{ext}^{\alpha \beta } = \sum\limits_{n = 1}^{ + \infty } {{G^n}h_n^{\alpha \beta },} $$

((29))

where the first term is composed of the result given by Eqs. (26, 27, 28). In this article, we shall always understand the infinite sums such as the one in Eq. (29) in the sense of formal power series, i.e. as an ordered collection of coefficients, e.g., $R_{\ell m\omega }^{\text{in}}\;\text{at}\;r \sim {r_0}$. We do not attempt to control the mathematical nature of the series and refer to the mathematical-physics literature for discussion (in the present context, see notably Refs. [45, 64]).

The post-Minkowskian solution

We insert the ansatz (29) into the vacuum Einstein field equations (12, 13), i.e. with τ^αβ = c⁴/(16πG)Λ^αβ, and we equate term by term the factors of the successive powers of our book-keeping parameter G. We get an infinite set of equations for each of the $R_{\ell m\omega }^{\text{in}}$,

$$\square h_n^{\alpha \beta } = \Lambda _n^{\alpha \beta }[{h_1},{h_2},...,{h_{n - 1}}],$$

((30))

$${\partial _\mu }h_n^{\alpha \mu } = 0.$$

((31))

The right-hand side of the wave equation (30) is obtained from inserting the previous iterations, up to the order n — 1, into the gravitational source term. In more details, the series of equations (30) reads

$$\square h_2^{\alpha \beta } = {N^{\alpha \beta }}[{h_1},{h_1}],$$

((32))

$$\square h_3^{\alpha \beta } = {M^{\alpha \beta }}[{h_1},{h_1},{h_1}] + {N^{\alpha \beta }}[{h_1},{h_2}] + {N^{\alpha \beta }}[{h_2},{h_1}],$$

((33))

$$\begin{gathered} \square h_4^{\alpha \beta } = {L^{\alpha \beta }}[{h_1},{h_1},{h_1},{h_1}] \hfill \\ \;\;\;\;\;\;\;\;\;\; + {M^{\alpha \beta }}[{h_1},{h_1},{h_2}] + {M^{\alpha \beta }}[{h_1},{h_2},{h_1}] + {M^{\alpha \beta }}[{h_2},{h_1},{h_1}] \hfill \\ \;\;\;\;\;\;\;\;\;\; + {N^{\alpha \beta }}[{h_2},{h_2}] + {N^{\alpha \beta }}[{h_1},{h_3}] + {N^{\alpha \beta }}[{h_3},{h_1}] \hfill \\ \;\;\;\;\;\;\;\;\;\; \vdots \hfill \\ \end{gathered} $$

((34))

The quadratic, cubic and quartic pieces of Λ^αβ are defined by Eq. (16).

Let us now proceed by induction. Some n being given, we assume that we succeeded in constructing, from the linearized coefficient h₁, the sequence of post-Minkowskian coefficients h₂, h₃, ..., h_n-1, and from this we want to infer the next coefficient h_n. The right-hand side of Eq. (30), Λ ^αβ _n , is known by induction hypothesis. Thus the problem is that of solving a wave equation whose source is given. The point is that this wave equation, instead of being valid everywhere in $X_{\ell m\omega }^{\text{in}} (z)$, is correct only outside the matter (r > a), and it makes no sense to solve it by means of the usual retarded integral. Technically speaking, the right-hand side of Eq. (30) is composed of the product of many multipole expansions, which are singular at the origin of the spatial coordinates r = 0, and which make the retarded integral divergent at that point. This does not mean that there are no solutions to the wave equation, but simply that the retarded integral does not constitute the appropriate solution in that context.

What we need is a solution which takes the same structure as the source term Λ_n^αβ, i.e. is expanded into multipole contributions, with a singularity at r = 0, and satisfies the d’Alembertian equation as soon as r > 0. Such a particular solution can be obtained, following the suggestion in Ref. [14], by means of a mathematical trick in which one first “regularizes” the source term And by multiplying it by the factor r^B, where $X_{\ell m\omega }^{\text{in}}$. Let us assume, for definiteness, that Λ ^αβ _n is composed of multipolar pieces with maximal multipolarity l_max. This means that we start the iteration from the linearized metric (26, 27, 28) in which the multipolar sums are actually finite⁷. The divergences when r → 0 of the source term are typically power-like, say 1/r_k (there are also powers of the logarithm of r), and with the previous assumption there will exist a maximal order of divergency, say k_max. Thus, when the real part of B is large enough, i.e. $X_{\ell}$, the “regularized” source term r^BΛ ^αβ _n is regular enough when r → 0 so that one can perfectly apply the retarded integral operator. This defines the B-dependent retarded integral

$${I^{\alpha \beta }}(B) \equiv \square _{ret}^{ - 1}[{\tilde r^B}\Lambda _n^{\alpha \beta }],$$

((35))

where the symbol □ ^-1_ret stands for the retarded integral (21). It is convenient to introduce inside the regularizing factor some arbitrary constant length scale r₀ in order to make it dimensionless. Everywhere in this article we pose

$$\tilde r \equiv \frac{r}{{{r_0}}}.$$

((36))

The fate of the constant r₀ in a detailed calculation will be interesting to follow, as we shall see, because it provides some check that the calculation is going well. Now the point for our purpose is that the function I^αβ(B) on the complex plane, which was originally defined only when ${\xi _\ell }$, admits a unique analytic continuation to all values of ${X _\ell }$ except at some integer values. Furthermore, the analytic continuation of I^αβ (B) can be expanded, when B → 0 (namely the limit of interest to us) into a Laurent expansion involving in general some multiple poles. The key idea, as we shall prove, is that the finite part, or the coefficient of the zeroth power of B in that expansion, represents the particular solution we are looking for. We write the Laurent expansion of I^αβ(B), when B → 0, in the form

$${I^{\alpha \beta }}(B) = \sum\limits_{p = {p_0}}^{ + \infty } {\iota _p^{\alpha \beta }{B^p}} ,$$

((37))

where ${\xi _\ell }$, and the various coefficients ι ^αβ _p are functions of the field point (x, t). When p₀ ≤ 1 there are poles; -p₀, which depends on n, refers to the maximal order of the poles. By applying the box operator onto both sides of (37), and equating the different powers of B, we arrive at

$$\begin{gathered} {p_0} \leqslant p \leqslant - 1\;\;\; \Rightarrow \;\;\;\square \iota _p^{\alpha \beta } = 0, \hfill \\ \;\;\;\;\;\;\;\;\;\;p \geqslant 0\;\; \Rightarrow \;\;\;\square \iota _p^{\alpha \beta } = \frac{{{{(\ln r)}^p}}}{{p!}}\Lambda _n^{\alpha \beta }. \hfill \\ \end{gathered} $$

((38))

As we see, the case p = 0 shows that the finite-part coefficient in Eq. (37), namely ι ^αβ₀ , is a particular solution of the requested equation: □ι ^αβ₀ = Ant. Furthermore, we can prove that this term, by its very construction, owns the same structure made of a multipolar expansion singular at r=0.

Let us forget about the intermediate name ι ^αβ₀ , and denote, from now on, the latter solution by u ^αβ _n ≡ ι ^αβ₀ , or, in more explicit term

$$u_n^{\alpha \beta } = \mathcal{F}{\mathcal{P}_{B = 0}}\square _{ret}^{ - 1}[{\tilde r^B}\Lambda _n^{\alpha \beta }],$$

((39))

where the finite-part symbol $\mathcal{O}(\epsilon)$ means the previously detailed operations of considering the analytic continuation, taking the Laurent expansion, and picking up the finite-part coefficient when B → 0. The story is not complete, however, because u ^αβ _n does not fulfill the constraint of harmonic coordinates (31); its divergence, say w ^α _n — ∂_μu ^αμ _n , is different from zero in general. From the fact that the source term is divergence-free in vacuum, ∂_μΛ ^αμ _n = 0 (see Eq. (18)), we find instead

$$\omega _n^{\alpha \beta } = \mathcal{F}{\mathcal{P}_{B = 0}}\square _{ret}^{ - 1}[B{\tilde r^B}\tfrac{{{n_i}}}{r}\Lambda _n^{\alpha i}].$$

((40))

The factor B comes from the differentiation of the regularization factor $\mathcal{O}(\epsilon^2)$. So, w ^α _n is zero only in the special case where the Laurent expansion of the retarded integral in Eq. (40) does not develop any simple pole when B → 0. Fortunately, when it does, the structure of the pole is quite easy to control. We find that it necessarily consists of a solution of the source-free d’Alembertian equation, and, what is more (from its stationarity in the past), the solution is a retarded one. Hence, taking into account the index structure of w ^α _n , there must exist four STF-tensorial functions of the retarded time u = t — r/c, say N_L(u), P_L (u), Q_L(u) and R_L(u), such that

$$\begin{gathered} \omega _n^0 = \sum\limits_{l = 0}^{ + \infty } {{\partial _L}} [{r^{ - 1}}{N_L}(u)], \hfill \\ \omega _n^i = \sum\limits_{l = 0}^{ + \infty } {{\partial _{iL}}} [{r^{ - 1}}{P_L}(u)] + \sum\limits_{l = 1}^{ + \infty } {\{ {\partial _{L - 1}}} [{r^{ - 1}}{Q_{iL - 1}}(u)] + {\varepsilon _{iab}}{\partial _{aL - 1}}[{r^{ - 1}}{R_{bL - 1}}(u)]\} . \hfill \\ \end{gathered} $$

((41))

From that expression we are able to find a new object, say v_n, which takes the same structure as w ^α _n (a retarded solution of the source-free wave equation) and, furthermore, whose divergence is exactly the opposite of the divergence of u ^αβ _n , i.e. ∂_μv ^αμ _n = -w ^α _n . Such a v ^αβ _n is not unique, but we shall see that it is simply necessary to make a choice for v ^αβ _n (the simplest one) in order to obtain the general solution. The formulas that we adopt are

$$\begin{gathered} \upsilon _n^{00} = {r^{ - 1}}{N^{( - 1)}} + {\partial _a}[{r^{ - 1}}( - N_a^{( - 1)} + C_a^{( - 2)} - 3{P_a})], \hfill \\ \upsilon _n^{0i} = {r^{ - 1}}( - {Q_i}^{( - 1)} + 3P_i^{(1)}) - {\varepsilon _{iab}}{\partial _a}[{r^{ - 1}}R_b^{( - 1)}] - \sum\limits_{l = 2}^{ + \infty } {{\partial _{L - 1}}} [{r^{ - 1}}{N_{iL - 1}}], \hfill \\ \upsilon _n^{ij} = {\delta _{ij}}{r^{ - 1}}P + \sum\limits_{l = 2}^{ + \infty } \{ 2{\delta _{ij}}{\partial _{L - 1}}[{r^{ - 1}}{P_{L - 1}}] - 6{\partial _{L - 2(i}}[{r^{ - 1}}{P_{j)L - 2}}] \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {\partial _{L - 2}}[{r^{ - 1}}(N_{ijL - 2}^{(1)} + 3P_{ijL - 2}^{(2)} - {Q_{ijL - 2}})] - 2{\partial _{aL - 2}}[{r^{ - 1}}{\varepsilon _{ab(i}}{R_{j)bL - 2}}]\} . \hfill \\ \end{gathered} $$

((42))

Notice the presence of anti-derivatives, denoted, e.g., by $\mathcal{O}(\epsilon)$; there is no problem with the limit v → — ∞ since all the corresponding functions are zero when $\mathcal{O}(\epsilon^{2\ell+2})$. The choice made in Eqs. (42) is dictated by the fact that the 00 component involves only some monopolar and dipolar terms, and that the spatial trace ii is monopolar: v ⁱⁱ _n = -3r^-1P. Finally, if we pose

$$h_n^{\alpha \beta } = u_n^{\alpha \beta } + \upsilon _n^{\alpha \beta },$$

((43))

we see that we solve at once the d’Alembertian equation (30) and the coordinate condition (31). That is, we have succeeded in finding a solution of the field equations at the nth post-Minkowskian order. By induction the same method applies to any order n, and, therefore, we have constructed a complete post-Minkowskian series (29) based on the linearized approximation h ^αβ₁ given by (26, 27, 28). The previous procedure constitutes an algorithm, which could be implemented by an algebraic computer programme.