Analytic Approximations in GR and Gravitational Waves

TL;DR

This paper reviews analytic approximation methods in general relativity, focusing on post-Newtonian, gravitational self-force, and post-Minkowskian approaches, highlighting their roles in modeling gravitational waves from compact binaries.

Contribution

It provides an overview of recent developments and interfaces among key approximation methods used in gravitational wave analysis within general relativity.

Findings

Post-Newtonian approximation effectively models early inspiral phases.

Interfaces between PN, GSF, and PM methods enhance gravitational wave modeling.

Recent developments improve the accuracy of gravitational wave signal predictions.

Abstract

Analytic approximation methods in general relativity play a very important role when analyzing the gravitational wave signals recently discovered by the LIGO and Virgo detectors. In this contribution, we present the state-of-the-art and some recent developments in the famous post-Newtonian (PN) or slow-motion approximation, which has successfully computed the equations of motion and the early inspiral phase of compact binary systems. We discuss also some interesting interfaces between the PN and the gravitational self-force (GSF) approach based on black-hole perturbation theory, and between PN and the post-Minkowskian (PM) approximation, namely a non-linearity expansion valid for weak field and possibly fast-moving sources.

Figures

Equations

$$\frac{\mathrm{d}{\mathbf{v}}_a}{\mathrm{d}t}$$

$$+\mathcal{O}\left(\frac{1}{c^9}\right).$$

$${\mathbf{A}}_a^{\text{N}}=-\sum _{b\ne a}\frac{G{m}_b}{r_{ab}^2}{\mathbf{n}}_{ab}.$$

$${\mathbf{A}}_a^{\text{1PN}}=$$

$$\displaystyle\qquad\qquad-4\sum_{c\neq a}\frac{Gm_{c}}{r_{ac}}-\sum_{c\neq b}\frac{Gm_{c}}{r_{bc}}\left(1-\frac{r_{ab}}{2r_{bc}}(n_{ab}n_{bc})\right)\biggr{]}$$

$$\displaystyle+\sum_{b\neq a}\frac{Gm_{b}}{r_{ab}^{2}}\bm{v}_{ab}\bigl{[}4(n_{ab}v_{a})-3(n_{ab}v_{b})\bigr{]}-\frac{7}{2}\sum_{b\neq a}\sum_{c\neq b}\frac{G^{2}m_{b}m_{c}}{r_{ab}r_{bc}^{2}}\bm{n}_{bc}\,,$$

$$S=\frac{c^{3}}{16\pi G}\int\mathrm{d}^{4}x\,\sqrt{-g}\Bigl{[}R-\frac{1}{2}g_{\mu\nu}\Gamma^{\mu}\Gamma^{\nu}\Bigr{]}-\sum_{a=1}^{N}m_{a}c^{2}\int\mathrm{d}t\,\sqrt{-(g_{\mu\nu})_{a}\,v_{a}^{\mu}v_{a}^{\nu}/c^{2}}\,.$$

$${\bar{g}}_{\mu \nu }(\mathbf{x},t)=g_{\mu \nu }[\mathbf{x};{\mathbf{y}}_a(t),{\mathbf{v}}_a(t),\cdots ].$$

$$\frac{\delta S_{\text{F}}}{\delta {\mathbf{y}}_a}=0.$$

$$S_{\text{F}}^{\text{tail}}=\frac{G^{2}M}{5c^8}{\int }_{-\infty }^{+\infty }\mathrm{d}t{I}_{ij}^{(3)}(t){\int }_0^{+\infty }\mathrm{d}\tau \ln \left(\frac{\tau }{{\tau }_0}\right)\left[I_{ij}^{(4)}(t-\tau )-I_{ij}^{(4)}(t+\tau )\right].$$

$$S_{\text{F}}^{\text{tail}}=\frac{G^{2}M}{5c^8}{\mathrm{Pf}}_{{\tau }_0}\int \int \frac{\mathrm{d}t\mathrm{d}{t}^{\prime }}{|t-t^{\prime }|}{I}_{ij}^{(3)}(t)I_{ij}^{(3)}(t^{\prime }).$$

$$\delta L_{\text{F}}=\frac{G^{4}(m_{1}+m_{2})\,m_{1}^{2}m_{2}^{2}}{c^{8}r_{12}^{4}}\Bigl{(}\delta_{1}(n_{12}v_{12})^{2}+\delta_{2}v_{12}^{2}\Bigr{)}\,.$$

$$L_{\text{F}}^{\text{HR}}=\mathop{\text{{\rm FP}}}_{B=0}\int_{r>\mathcal{R}}\mathrm{d}^{3}\mathbf{x}\,\Bigl{(}\frac{r}{r_{0}}\Bigr{)}^{B}F(\mathbf{x})\,.$$

$$L_{\text{F}}^{\text{DR}}={\int }_{r>\mathcal{R}}\frac{{\mathrm{d}}^d\mathbf{x}}{{\ell }_0^{d-3}}{F}^{(d)}(\mathbf{x}),$$

$$L_{\text{F}}^{\text{DR}}-L_{\text{F}}^{\text{HR}}=\sum_{q}\biggl{[}\frac{1}{(q-1)\varepsilon}-\ln\left(\frac{r_{0}}{\ell_{0}}\right)\biggr{]}\int\mathrm{d}\Omega_{2+\varepsilon}\,\varphi^{(\varepsilon)}_{3,q}(\mathbf{n})+\mathcal{O}\left(\varepsilon\right)\,.$$

$$\tau_{0}^{\text{DR}}=\frac{2\ell_{0}}{c\sqrt{4\pi}}\,\text{exp}\Bigl{[}\frac{1}{2\varepsilon}-\frac{1}{2}\gamma_{\text{E}}-\frac{41}{60}\Bigr{]}\,,$$

$${\delta }_1=-\frac{2179}{315},{\delta }_2=\frac{192}{35},$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=-\mathcal{F}.$$

$$\varphi =\int \omega \mathrm{d}t=-\int \frac{\omega }{\mathcal{F}}\frac{\mathrm{d}E}{\mathrm{d}\omega }\mathrm{d}\omega .$$

$$E$$

$$\displaystyle\qquad\quad+\left(-\frac{675}{64}+\biggl{[}\frac{34445}{576}-\frac{205}{96}\pi^{2}\biggr{]}\nu-\frac{155}{96}\nu^{2}-\frac{35}{5184}\nu^{3}\right)x^{3}$$

$$+(-\frac{3969}{128}+\left[-\frac{123671}{5760}+\frac{9037}{1536}{\pi }^2+\frac{896}{15}{\gamma }_{\text{E}}+\frac{448}{15}\ln (16x)\right]\nu$$

$$\displaystyle\qquad\qquad\quad\left.+\left[-\frac{498449}{3456}+\frac{3157}{576}\pi^{2}\right]\nu^{2}+\frac{301}{1728}\nu^{3}+\frac{77}{31104}\nu^{4}\right)x^{4}\biggr{\}}\,,$$

$$I_{i_1{i}_2}=m\nu \left(A{x}_{⟨i_1{i}_2⟩}+B\frac{r^2}{c^2}{v}_{⟨i_1{i}_2⟩}+\frac{G^2{m}^2\nu }{c^{5}r}C{x}_{⟨i_1}{v}_{i_2⟩}\right),$$

$$A$$

$$+{\gamma }^3(\frac{395899}{13200}-\frac{428}{105}\ln \left(\frac{r}{r_0}\right)+\left[\frac{3304319}{166320}-\frac{44}{3}\ln \left(\frac{r}{{r^{\prime }}_0}\right)\right]\nu$$

$$+\frac{162539}{16632}{\nu }^2+\frac{2351}{33264}{\nu }^3),$$

$$B$$

$$+{\gamma }^2\left(-\frac{357761}{19800}+\frac{428}{105}\ln \left(\frac{r}{r_0}\right)-\frac{92339}{5544}\nu +\frac{35759}{924}{\nu }^2+\frac{457}{5544}{\nu }^3\right),$$

$$C$$

$$U_{i_1{i}_1}(t)$$

$$+\frac{G^2{M}^2}{c^6}{\int }_0^{+\infty }\mathrm{d}\tau I_{i_1{i}_2}^{(5)}(t-\tau )[2{\ln }^2\left(\frac{c\tau }{2b_0}\right)+\frac{11}{3}\ln \left(\frac{c\tau }{2b_0}\right)$$

$$-\frac{214}{105}\ln \left(\frac{c\tau }{2r_0}\right)+\frac{124627}{22050}]$$

$$+\frac{G^3{M}^3}{c^9}{\int }_0^{+\infty }\mathrm{d}\tau I_{i_1{i}_2}^{(6)}(t-\tau )[\frac{4}{3}{\ln }^3\left(\frac{c\tau }{2b_0}\right)+\frac{11}{3}{\ln }^2\left(\frac{c\tau }{2b_0}\right)$$

$$+\frac{124627}{11025}\ln \left(\frac{c\tau }{2b_0}\right)-\frac{428}{105}\ln \left(\frac{c\tau }{2b_0}\right)\ln \left(\frac{c\tau }{2r_0}\right)$$

$$-\frac{1177}{315}\ln \left(\frac{c\tau }{2r_0}\right)+\frac{129268}{33075}+\frac{428}{315}{\pi }^2].$$

$$\mathcal{F}$$

$$+\left(-\frac{44711}{9072}+\frac{9271}{504}\nu +\frac{65}{18}{\nu }^2\right)x^2+\left(-\frac{8191}{672}-\frac{583}{24}\nu \right)\pi x^{5/2}$$

$$+(\frac{6643739519}{69854400}+\frac{16}{3}{\pi }^2-\frac{1712}{105}{\gamma }_{\text{E}}-\frac{856}{105}\ln (16x)$$

Full text

\bodymatter

1 Methods to generate gravitational wave templates

The LIGO & Virgo detectors have opened up a fantastic new avenue in Astronomy with the discovery of gravitational waves (GWs) generated by the orbital motion and merger of binary black hole and neutron star systems.[1, 2] This also highlights the crucial role played by analytic approximation methods in general relativity (GR), since they permit an accurate description of the two-body problem in GR, which is of direct use in the data analysis of the detectors.[44]

The most important method in this respect is the post-Newtonian (PN) approximation, which is an expansion when the slowness parameter $\epsilon_{\text{PN}}=v/c$ of the compact binary system tends to zero, where $v$ is the relative orbital velocity and $c$ the speed of light. For gravitationally bound systems such as compact binaries on quasi-circular orbits, the PN approximation comes along with the post-Minkowskian (PM) one, namely a non-linearity expansion around the Minkowski background, with small expansion parameter $\gamma_{\text{PM}}=Gm/(rc^{2})$, where $r$ is the size of the orbit and $m$ the total mass of the source. Indeed, in the bounded case we have $\gamma_{\text{PM}}\sim\epsilon_{\text{PN}}^{2}$. However, the most important physical application of the PM approximation is for unbound orbits, when $\gamma_{\text{PN}}$ and $\epsilon_{\text{PM}}$ are unrelated, i.e., the problem of scattering of ultra-relativistic particles ($\epsilon_{\text{PN}}\lesssim 1$) and small deviation angle. The PM approximation is sometimes called the weak-field fast-moving approximation.

Black hole perturbation theory constitutes another large body of analytic approximations in GR. In the context of compact binary systems, this approximation is important, first, for analyzing the post-merger waveform of two black holes (BHs) during the so-called ringdown phase, when the newly formed BH emits quasi-normal mode radiation, and, secondly, for describing the dynamics and GWs of asymmetric compact binaries, i.e., endowed with an extreme mass ratio, $\nu=m_{1}m_{2}/(m_{1}+m_{2})^{2}\ll 1$. In the latter case the perturbation method takes the more suggestive name of gravitational self-force (GSF), since it is concerned with the modifications of the background geometry of the larger BH and of the geodesic motion of the particle, due to the self field generated by the particle itself.

The domains of validity of these approximation methods, together with interesting mutual overlaps, are illustrated in Figs. 1 and 2. Also shown in Fig. 1 is some comparison with numerical relativity (NR), which succeeded at solving the long standing binary black hole challenge.[93, 45, 4] At first sight it could seem that NR would be able to tackle and solve the complete problem of the inspiral, merger and ringdown for two compact objects. However, in order to monitor the early inspiral of two neutron stars, thousands of orbital cycles have to be computed with high precision. Then the computing times of NR become prohibitively long, and the precision of the NR simulation will never be competitive with that provided by the PN approximation. A fact that is of uttermost importance for building GW templates (and was not a priori obvious several years ago[42]), is that the overlap between PN and NR exists and is quite significant. On the other hand, when the mass ratio between the two bodies is extreme, the full NR approach is unfeasible, due to the different length scales corresponding to the very different physical sizes of the compact bodies. We shall discuss in this article the important intersections between PN and GSF as shown in Fig. 1, and between PN and PM, see Fig. 2.

The GW templates are defined as the theoretical prediction from GR, and weighted in the Fourier domain by the detector’s spectral density of noise. The templates are cross-correlated with the detector’s output, and the correlation builds up when a good match occurs between a particular template and the real signal.[44] This technique is highly sensitive to the phase evolution of the signal, which, in PN templates of compact binary coalescence, is computed from the energy balance between the decay of the binary’s energy and (minus) the GW flux.

In principle, as there is a significant overlap between the PN and NR regimes, the templates are obtained by matching together the best PN waveform for the inspiral phase (currently known to order 3.5PN) to a highly accurate numerical waveform for the merger and ringdown phases.[67] For low mass compact binaries, such as double neutron star systems, the detectors are mostly sensitive to the inspiral phase prior to the final coalescence, and the currently known analytical PN templates are accurate enough for detection, at least when the compact bodies have moderate spins. Thus, the data analysis of neutron star binaries events such as GW170817, is essentially based on the 3.5PN templates.

For larger masses, like BH binary events such as GW150914, the merger occurs at lower frequency, right in the middle of the detector’s frequency band. Since only a few orbital cycles before the merger are seen, the match between NR and the PN is not very good and the GW templates are essentially based on the NR results. Nevertheless it is worth pointing out that the mere zeroth-order Newtonian waveform, i.e., based on the Einstein quadrupole formula, gives a reasonably satisfying physical interpretation of the signal even for GW150914!

In practice, for the data analysis of large mass BH binary events, it is important to dispose of analytic rather than numerical templates, since the cross correlations must be performed with many templates associated by many trial parameters (masses and spins). In this case the templates are based on effective analytical methods that interpolate somehow between the initial PN and final NR phases. Two such techniques play a key role in the data analysis of the BH events. One is called the hybrid inspiral-merger-ringdown (IMR or IMR-Phenom) waveform and is constructed by matching together the PN and NR waveforms in an overlapping time interval described phenomenologically.[3] The other technique consists of recasting the real two-body dynamics, as given by PN theory, into a simpler one-body dynamics described in a non perturbative analytic way. The so-called effective-one-body (EOB) waveform obtained in that way extends the domain of validity of the PN approximation (because it is non perturbative) and can therefore be compared and matched to the NR waveform.[43] The IMR and EOB waveforms are extensively used in the LIGO & Virgo data analysis of the recent binary BH events.

2 State-of-the-art on equations of motion

The equations of motion (EOM) of a self-gravitating $N$-body system are written in PN like form, up to 4PN order, using one’s favorite coordinate system in GR, as

[TABLE]

The first term is of course, the usual Newtonian acceleration of $N$ “planets”,

[TABLE]

The historical works in GR started in the early days of this theory, and solved the problem of the EOM at the 1PN level beyond the Newtonian term.[79, 58] This famous Lorentz-Droste-Einstein-Infeld-Hoffmann 1PN correction is fully given as

[TABLE]

where we denote $\bm{v}_{a}=\mathrm{d}\bm{y}_{a}/\mathrm{d}t$, $\bm{v}_{ab}=\bm{v}_{a}-\bm{v}_{b}$, $r_{ab}=|\bm{y}_{a}-\bm{y}_{b}|$, $\bm{n}_{ab}=(\bm{y}_{a}-\bm{y}_{b})/r_{ab}$ and the parenthesis indicate the usual Euclidean scalar product, e.g. $(n_{ab}v_{b})=\bm{n}_{ab}\cdot\bm{v}_{b}$.

Up to the 2PN level the system is conservative, i.e., admits the ten invariants associated with the symmetries of the Poincaré group. The first dissipative effect appears at the 2.5PN order and features the radiation reaction damping of the system by GW emission. The 2PN and 2.5PN approximations were motivated by the Hulse-Taylor binary pulsar and worked out at the time of its discovery.[87, 50, 47, 98, 33, 72] Later the motivation for the 3PN EOM came from the development of the LIGO & Virgo detectors and the need of accurate GW templates for inspiralling compact binaries. The 3PN EOM took some time to be fully understood and completed,[73, 30, 54, 51, 37, 23, 71, 69, 60] together with the relatively easier dissipative 3.5PN term.[88, 76, 86, 70]

Three techniques have been undertaken to obtain the 4PN EOM. One is based on the Arnowitt-Deser-Misner (ADM) Hamiltonian formalism of GR in ADM coordinates,[74, 75, 52, 53] and has led to complete results but for the appearance of one “ambiguity” parameter. The second technique is based on the Fokker action of GR in harmonic coordinates,[7, 9, 8, 83, 10] and has obtained complete results, i.e., free of any ambiguity parameter. The third one is the effective field theory (EFT),[61, 62, 65, 59] which yielded partial results up to now (the terms $\propto G^{4}$ still being in progress), but is expected to also be free of any ambiguity parameter.[90]

In this section we describe the approach based on the Fokker action in harmonic coordinates. We start with the gravitation-plus-matter action of GR, in which the gravitational piece includes the usual harmonic gauge-fixing term, and the matter term is that of $N$ particles without spins, and with negligible internal structure:

[TABLE]

Here $\Gamma^{\mu}=g^{\rho\sigma}\Gamma^{\mu}_{\rho\sigma}$ and we use, for practical calculations, the Landau-Lifshitz form of the action (i.e., modulo a total divergence).

The Fokker action is obtained when we insert back into (4) an explicit PN solution of the corresponding gauge-fixed Einstein field equations, and given by an explicit functional of the particle’s trajectories, i.e.,

[TABLE]

The ellipsis indicate that the metric also depends on accelerations $\bm{a}_{a}=\mathrm{d}\bm{v}_{a}/\mathrm{d}t$, derivatives of accelerations $\bm{b}_{a}=\mathrm{d}\bm{a}_{a}/\mathrm{d}t$, etc., since we do not perform any replacements of accelerations when iterating the Einstein field equations, the EOM being considered off-shell at this stage. Substituting $\bar{g}_{\mu\nu}$ into Eq. (4) defines the Fokker action $S_{\text{F}}[\bm{y}_{a},\bm{v}_{a},\cdots]$, and the EOM of the self-gravitating system of particles are obtained as the (generalized Euler-Lagrange) equations

[TABLE]

Once they have been constructed, the EOM can be order reduced by replacing all the higher-order accelerations by their expressions coming from the lowest-order PN equations. The Fokker action describes only the conservative dynamics, and the dissipative effects have to be added separately in the EOM. Note that the Fokker action is equivalent to the EFT action in the “tree-level” approximation, in which we neglect quantum loops.

An interesting feature of the local (near zone) 4PN dynamics is that there is an imprint of GW tails propagating at infinity. The tails are secondary non-linear waves caused by backscattering of linear waves onto the space-time curvature generated by the total mass $M$ of the source. Part of the effect can be seen as a tail-induced modification of the leading 2.5PN radiation reaction force at the relative 1.5PN order.[21, 16, 18] However, associated with this dissipative piece, there exists also a conservative effect which thus enters into the 4PN conservative dynamics. Its contribution to the Fokker action reads as[62, 52, 7, 65]

[TABLE]

Here $I_{ij}=\sum_{a}m_{a}y_{a}^{\langle i}y_{a}^{j\rangle}$ is the Newtonian quadrupole moment of the system (the angular brackets refer to the symmetric-trace-free projection), the superscript $(n)$ denotes multiple time derivatives, and $\tau_{0}$ is an arbitrary constant. To leading order the total mass $M$ reduces to $\sum_{a}m_{a}$, but at higher order it should involve the contribution of the gravitational binding energy of the particles. An elegant rewriting of Eq. (7) is with the Hadamard “Partie finie” (Pf) integral,

[TABLE]

Due to this conservative tail contribution, the 4PN dynamics is non-local in time, and this entails subtleties in the derivation of the invariants of motion, which have been recently fully elucidated.[52, 53, 9, 41]

The calculation crucially relies on the systematic use of dimensional regularization (DR), to cure both ultra-violet (UV) divergences due to the model of point particles adopted to describe the compact objects, and infra-red (IR) divergences that start appearing precisely at the 4PN order and are associated with GW tails. We are here borrowing DR from EFT and quantum field theory, and we use it in the classical $N$-body problem as a mean to preserve the diffeomorphism invariance of GR. For this reason we conjecture that DR is the only known regularization technique able to successfully solve the problem of EOM in high PN approximations.

In an initial calculation (valid for two particles, $N=2$), we used DR for UV divergences but a variant of the Hadamard regularization (HR) for IR divergences.[7, 9] Based on some trial calculations, using various types of regularizations, we conjectured that the results of different IR regularizations will physically differ by at most two parameters called ambiguities and denoted $\delta_{1}$ and $\delta_{2}$.[9] This finding was in agreement with an earlier suggestion.[53] Modulo unphysical shifts of the trajectories, the two offending ambiguity terms in the Fokker Lagrangian ($S_{\text{F}}=\int\mathrm{d}tL_{\text{F}}$) turn out to appear at the difficult $G^{4}$ level and be of the form

[TABLE]

To determine what the values of these ambiguities are we embarked on the DR treatment of the IR divergences. Consider a typical term in the Fokker Lagrangian with non-compact support and generic function $F$, which diverges at infinity. With HR such term is treated as

[TABLE]

Since we focus on IR divergences we consider only the far zone contribution $r>\mathcal{R}$, where $\mathcal{R}$ denotes an arbitrary large radius, typically the inner radius of the wave zone. The Finite Part (FP) operation is closely related to the Hadamard partie finie Pf, and consists of applying analytical continuation in $B\in\mathbb{C}$, expanding the integral when $B$ tends to zero, and keeping only the coefficient of the zero-th power of $B$ in that expansion (discarding any strictly positive or negative power of $B$). On the other hand, with DR the same term is treated as

[TABLE]

where $F^{(d)}$ is the $d$-dimensional analogue of the generic function $F$ in (10), and where $\ell_{0}$ is the characteristic length scale associated with DR. We find that the difference between the two prescriptions is given by[8]

[TABLE]

The functions $\varphi^{(\varepsilon)}_{p,q}$ represent the coefficients of $r^{-p-q\varepsilon}$ (with $p,q\in\mathbb{Z}$) in the expansion of $F^{(d)}$ when $r\to+\infty$ along the direction $\mathbf{n}$. We pose $\varepsilon=d-3$ and neglect the terms dying with $\varepsilon\to 0$. The angular integration in (12) is over the sphere in $d-1=2+\varepsilon$ dimensions. Notice that the result (12) depends only on the singular coefficients $\varphi^{(\varepsilon)}_{3,q}$ of the expansion of $F^{(d)}$ at infinity, and that the arbitrary scale $\mathcal{R}$ has disappeared from it.

The formula (12) contains an IR pole $\propto 1/\varepsilon$. In the language of the EFT, this pole comes from the “potential mode” contribution. However, we have also to take into account the 4PN tail effect given by Eq. (8) but in $d$ dimensions. In the EFT language this will correspond to the “radiation” contribution and should also contain a pole $\propto 1/\varepsilon$, but this time of UV type. Our explicit calculation has shown that the two IR and UV poles exactly cancel out (modulo unphysical shifts of the trajectories),[8, 83] in complete agreement with general arguments within the EFT.[90] The 4PN tail term in $d$ dimensions takes the same form as (8) but with the arbitrary scale $\tau_{0}$ determined by the DR scale $\ell_{0}$ as

[TABLE]

with $\gamma_{\text{E}}$ denoting the Euler constant, and the UV pole $\propto 1/\varepsilon$ cancelling the IR one in (12). The result (13) has been obtained thanks to a “matching” equation relating the near zone which is the domain of validity of the PN approximation, to the far zone where GW tails propagate. Finally we find that the rational fraction $-\frac{41}{60}$ in (13), is just the one necessary and sufficient to determine the values of the two ambiguity parameters in (9) as

[TABLE]

therefore resolving the problem of ambiguities. The values (14) are consistent with numerical and analytical GSF calculations of the energy and periastron advance for circular orbits in the small mass ratio limit.[52, 53, 8] Remarkably, the result (13) agrees with that of Galley et al.,[65] obtained by means of a diagrammatic evaluation of the tail term in $d$ dimensions with EFT methods. On the other hand, the lack of a consistent matching between the near and far zones in the ADM Hamiltonian formalism,[74, 75, 52, 53] and therefore a complete control of the tail term (8) including the final determination of Eq. (13), forces this formalism to be still plagued by one ambiguity parameter, denoted $C$ in[52].

3 State-of-the-art on GW generation

The two basic ingredients in the theoretical PN analysis correspond to the two sides of the energy balance equation obeyed by the binary’s orbital frequency and phase. Since the orbit will have circularized by radiation reaction at the time when the signal enters the detectors’ bandwidth there is no need to invoke the balance equation for the orbital angular momentum. Thus we just impose

[TABLE]

The energy $E$ is nothing but the Noetherian conserved energy $E$ associated with the Fokker Lagrangian computed in the previous section. On the other hand, the GW energy flux $\mathcal{F}$ on the right-hand side is obtained from a GW generation formalism. From Eq. (15) one deduces the time evolution of the binary’s orbital frequency $\omega$ and orbital phase $\phi$ by solving

[TABLE]

At the 4.5PN order, for circular orbits, the conserved energy function is given by

[TABLE]

where $m=m_{1}+m_{2}$ is the total mass, $\nu=m_{1}m_{2}/m^{2}$ is the symmetric mass ratio, $\gamma_{\text{E}}$ is Euler’s constant, and we employ for convenience the PN ordering parameter $x=(\frac{Gm\omega}{c^{3}})^{2/3}$ defined from the orbital frequency $\omega$ of the circular orbit, and which constitutes an invariant in a large class of coordinate systems. Notice that Eq. (3) is valid up to the 4.5PN order included, as there is no term at the 4.5PN order in the conserved energy for circular orbits.

The most complete formula for the GW flux is valid at the 3.5PN order beyond the Einstein quadrupole formula. However, the 4.5PN coefficient is also known,[84] while the 4PN coefficient is in progress. This formula has been obtained by application of a GW generation formalism based on a Multipolar-Post-Minkowskian (MPM) expansion for the external field of an isolated source,[20, 15, 21, 22] and followed by a matching to the inner (near zone) PN field of that source.[17, 19, 91, 32] The first important step in this computation is the obtention of the multipole moments of the source, $I_{i_{1}\cdots i_{\ell}}$ (mass type) and $J_{i_{1}\cdots i_{\ell}}$ (current type). The most difficult of these moments (because it necessitates the highest PN precision) is the mass quadrupole moment $I_{i_{1}i_{2}}$, given at 3.5PN order for quasi-circular orbits as

[TABLE]

where the terms are explicitly given by[39, 38]

[TABLE]

Here the PN ordering parameter is $\gamma=\frac{Gm}{rc^{2}}$ where $r$ is the separation distance in harmonic coordinates. Note the two constant scales entering the logarithmic terms at the 3PN order, one being the length scale $r_{0}$ coming from the MPM algorithm,[20] while the other one $r^{\prime}_{0}$ comes from the 3PN EOM in harmonic coordinates.[30]

The second step is the relationship between the multipole moments of the source, and the so-called “radiative” multipole moments parametrizing the observable GW at future null infinity. Such relationship involves in particular the well-known tail effects and their iterations. At the 4.5PN order the radiative mass quadrupole moment $U_{i_{1}i_{2}}$ is related to the mass quadrupole moment of the source $I_{i_{1}i_{2}}$ by

[TABLE]

For simplicity, we have not included here the non-linear memory effect which arises at 2.5PN order,[46, 106, 22, 102] as well as many instantaneous (non-tails) terms, that are relatively easy to compute. The terms at 1.5PN, 3PN and 4.5PN orders shown in (3) correspond to what can rightly be called the “tail”, the “tail-of-tail”, and the “tail-of-tail-of-tail”, respectively.[84] The expression (3) contains still another arbitrary scale $b_{0}$, parametrizing the coordinate transformation between harmonic coordinates and radiative coordinates. We find that the scale $b_{0}$ as well as the two previous scales $r_{0}$ and $r^{\prime}_{0}$ in Eq. (18) cleanly cancel out in the GW flux, expressed in terms of the invariant PN parameter $x$, which is finally given by[27, 26, 40, 39, 31, 38, 24, 25]

[TABLE]

This is valid up to 4.5PN order, with the notable exception that the 4PN coefficient, denoted $F_{\text{4PN}}$ in (3), is not yet known. However, from BH perturbation theory we know already the test mass limit of this coefficient, i.e., in the small mass ratio limit $\nu\to 0$:[96, 100, 101, 63, 64]

[TABLE]

Of course, this nice result from BH perturbation theory will have to be confirmed by PN theory, which will also be able to provide the mass ratio corrections $\mathcal{O}(\nu)$.

4 PN theory versus GSF theory

The conservative dynamics and GWs of compact binary systems in the extreme mass ratio limit, is the realm of the perturbative gravitational self force (GSF) theory.[55, 85, 94, 57, 66, 92] For the conservative dynamics, a comparison between GSF computations and traditional PN calculations was initiated some years ago,[56] applying to systems that are at once slowy moving and with extreme mass ratio, see the overlap region between PN and perturbation theory in Fig. 1. In recent years the possibility for this comparison has been dramatically extended. Such progress is due in large part due to high precision numerical and analytical computations from a self force perspective,[81, 80, 82, 95, 12, 14, 13, 99] and to extensive analytical computations within the PN approximation.[29, 28, 77, 34, 35]

For a particle moving on an exact circular orbit around a Schwarzschild BH (neglecting radiation reaction), one disposes of a very interesting quantity, which is the invariant associated with the helical Killing symmetry, appropriate for exact circular orbits.[56] This invariant, denoted $u_{1}^{t}$, is defined by

[TABLE]

where $u_{1}^{\mu}$ is the normalized four-velocity of the particle 1 (with mass $m_{1}\ll m_{2}$), $K^{\mu}$ is the helical Killing vector (HKV) and $K_{1}^{\mu}$ the HKV at the location of the particle. Adopting a coordinate system in which the HKV reads $K^{\mu}\partial_{\mu}=\partial_{t}+\omega\,\partial_{\varphi}$, where $\omega$ is the orbital frequency of the circular orbit, the invariant reduces to the time component of the four-velocity hence its name $u_{1}^{t}$, and we have

[TABLE]

where $(g_{\mu\nu})_{1}$ is the metric evaluated at the particle’s location, following a certain regularization (here $v_{1}^{\mu}=\mathrm{d}y_{1}^{\mu}/\mathrm{d}t$ denotes the coordinate velocity, i.e., $y_{1}^{0}=ct$ and $v_{1}^{0}=c$). The inverse of $u_{1}^{t}$ appears to be a redshift $z_{1}$, and sometimes $u_{1}^{t}$ itself is called the redshift. [For a generalization of the notion of redshift to eccentric orbits, see[5].] In the exact test mass limit $\nu=0$, the invariant reduces to the one appropriate to a Schwarzschild BH,

[TABLE]

where $y=(\frac{Gm_{2}\omega}{c^{3}})^{2/3}$ is the frequency-related PN parameter associated with the larger BH mass $m_{2}$. The GSF part is then defined as the coefficient of the mass ratio correction beyond Eq. (25),

[TABLE]

It is clear that with this approximation the symmetric mass ratio $\nu=\frac{m_{1}m_{2}}{m^{2}}$ can be replaced by the ordinary mass ratio $q=\frac{m_{1}}{m_{2}}$. In the PN approximation, the GSF part of the redshift factor appears to be an infinite PN series of the type

[TABLE]

We have included terms linear in the logarithm of $y$, but we neglect (just for this discussion) the higher powers of $\ln y$, that occur at large PN orders. Recall that the most general structure of the PN expansion involves any (integer) powers of the logarithm, $\sim(\ln y)^{k}y^{j+1}$.[20]

Here we report the results that have been obtained so far using the “traditional” PN approach. Recall that the PN method heavily relies on dimensional regularization (DR) to treat both UV and IR divergences, see Sec. 2. Another feature of the PN calculation is that it requires a machinery of tails and related non-linear effects, see Eq. 3. In the PN approach one computes $u_{1}^{t}$ as a redshift in harmonic coordinates using Eq. (24), and evaluates the metric at the particle’s location with DR. In that way the GSF redshift has been obtained up to 4PN order as[29, 28, 77, 12]

[TABLE]

In addition, PN theory has been able to fix the logarithmic term at the 5PN order, which is due to a subdominant tail effect, namely

[TABLE]

while the coefficient $\alpha_{5}$ is known only from GSF methods but has not yet been checked with PN theory. The results (4)–(29) are in full agreement with direct GSF computations. This constitutes a strong confirmation of the adequation of DR for traditional PN calculations, as well as of the procedure of subtraction of the singular field which is employed by GSF theory.

A feature of the PN expansion of the redshift factor at high orders is the appearance of half-integral PN approximations, say $\frac{n}{2}$PN. At first sight these terms sound surprising because the dynamics is purely conservative (exactly circular orbits with a HKV), and we are used to the fact that half-integral PN approximations like 2.5PN are associated with dissipative radiation reaction. The existence of such terms in the conservative redshift factor, starting at 5.5PN order, has been pointed out with numerical GSF methods,[99] and later it was proved that these terms originate from iterated non-linear tail effects, called “tails-of-tails”.[34, 35] The leading 5.5PN, next-to-leading 6.5PN and next-to-next-to-leading 7.5PN coefficients in this category have been found to be

[TABLE]

while the corresponding $\beta_{j}$’s are zero. Notice that 7.5PN is arguably the highest order ever reached by traditional PN methods. Again the PN results (30) are in full agreement with numerical and analytical results derived by GSF methods.

5 PN versus PM

The post-Minkowskian approximation has been developed in many pionneering works.[11, 104, 68, 89, 105, 6] Notably the gravitational scattering angle of two relativistic particles has been obtained up to 2PM order (quadratic in $G$).[103] Recently there has been a renewal of interest in the PM approximation. Ledvinka, Schäfer and Bičák[78] obtained a closed-form expression for the Hamiltonian of $N$ particles in the 1PM approximation, and new works appeared on the gravitational scattering angle and the link between the PM expansion and the EOB formalism.[48, 49] Here we outline our own contribution,[36] which concerns the comparison between 1PM and the recent 4PN calculation of the EOM of compact binaries. The Fig. 2 showed the domain of validity of the PM approximation versus that of the PN expansion.

At the 1PM approximation the field equations for a system of $N$ particles in harmonic coordinates, deduced from the gauge-fixed action (4), read

[TABLE]

where $\Box$ denotes the flat space-time d’Alembertian operator, $\delta^{(4)}$ is the four-dimensional Dirac function, $y_{a}^{\mu}$ are the particle’s worldlines and $u_{a}^{\mu}=\mathrm{d}y_{a}^{\mu}/(c\mathrm{d}\tau_{a})$ their normalized four-velocities, with the special-relativistic proper time $\mathrm{d}\tau_{a}=\sqrt{-\eta_{\mu\nu}\mathrm{d}y_{a}^{\mu}\mathrm{d}y_{a}^{\nu}/c^{2}}$. We solve Eqs. (31) by means of the standard Lienard-Wiechert procedure. Adopting a parametrization by the coordinate time $t$, i.e., such that $y_{a}^{\mu}=(ct,\bm{y}_{a})$, the retarded time $t_{a}^{\text{ret}}$ on the trajectory $a$ associated with the propagation from $a$ to the field point $x^{\mu}=(ct,\mathbf{x})$, is given by the implicit retardation equation $t_{a}^{\text{ret}}=t-r_{a}^{\text{ret}}/c$ with $r_{a}^{\text{ret}}=|\mathbf{x}-\bm{y}_{a}(t_{a}^{\text{ret}})|$. The solution of (31) is then

[TABLE]

where $k_{a}^{\mu}=[x^{\mu}-y_{a}^{\mu}(t_{a}^{\text{ret}})]/r_{a}^{\text{ret}}=(1,\bm{n}_{a}^{\text{ret}})$ is the Minkowski null vector between $a$ and the field point, $(ku)_{a}^{\text{ret}}=k^{a}_{\mu}u_{a}^{\mu}=\gamma_{a}(-1+\bm{n}_{a}^{\text{ret}}\cdot\bm{v}_{a}/c)$ is the usual redshift factor, with $\gamma_{a}=u_{a}^{0}$ the Lorentz factor and $v_{a}^{\mu}=cu_{a}^{\mu}/\gamma_{a}=(c,\bm{v}_{a})$, the velocities being computed at the retarded time $t_{a}^{\text{ret}}$.

However, we repeatedly use the fact that the accelerations are of order $G$ and therefore their contributions in (32) will be of order $G^{2}$, hence negligible with the 1PM approximation. Thus, we are allowed to assume that the four velocities $u_{a}^{\mu}$ and Lorentz factors $\gamma_{a}$ are constant. Furthermore, neglecting terms of order $G^{2}$, we can solve the retardation equation to get the retarded time $t_{a}^{\text{ret}}$, the distance $r_{a}^{\text{ret}}$, and the direction $\bm{n}_{a}^{\text{ret}}$, in terms of their current values at time $t$, i.e., the “instantaneous” distance $r_{a}=|\mathbf{x}-\bm{y}_{a}(t)|$ and direction $\bm{n}_{a}=[\mathbf{x}-\bm{y}_{a}(t)]/r_{a}$. In this way, Eq. (32) becomes equivalent to

[TABLE]

which is valid at any field point except at the singular locations of the particles. Nevertheless, we can easily extend its validity to the particles by using a self-field regularization. For this purpose, it is sufficient to discard the self-field contribution from the sum of particles. An explicit self-field regularization process yielding the same result was implemented in[6]. Therefore, at the location of particle $a$, we have

[TABLE]

where the sum runs over all particles except $a$, we pose $r_{ab}=|\bm{y}_{a}-\bm{y}_{b}|$, and denote $n_{ab}^{0}=0$ and $n_{ab}^{i}=[\bm{y}_{a}-\bm{y}_{b}]/r_{ab}$.

The EOM of the particles is just the geodesic equation, computed at the linearized order consistent with our approximation, and we obtain

[TABLE]

We use $\epsilon_{ab}=-(u_{a}u_{b})$ as a shorthand notation. Equivalently we have also the non-covariant form of the EOM (i.e., PN like form), in which we introduce the ordinary velocities and accelerations and the relevant Lorentz factors (with $\bm{v}_{ab}=\bm{v}_{a}-\bm{v}_{b}$),

[TABLE]

At the 1PM order the EOM are conservative, thus admit conserved integrals of energy, angular momentum and linear momentum. Indeed, the radiation reaction dissipative effects are at least 2PM, i.e., second order in $G$. The closed-form expression of the energy through 1PM reads $E=\sum_{a}m_{a}c^{2}\gamma_{a}+V$, with the first term being the usual special-relativistic expression and

[TABLE]

We have verified[36] that Eqs. (5) and (37) perfectly reproduce the PN results in harmonic coordinates, in the case of two particles ($N=2$) at the linear order in $G$ and up to the 4PN order.[8, 10]

Next we consider the Lagrangian associated with the EOM (35)–(5), for any $N$, in harmonic coordinates. The Lagrangian will be given by the special-relativistic expression plus terms of order $G$, and again, we neglect higher-order terms in $G$. In PN theory, it is known that the Lagrangian in harmonic coordinates is a generalized one, depending not only on positions and velocities $\bm{y}_{a},\bm{v}_{a}$ but also on accelerations $\bm{a}_{a}=\mathrm{d}\bm{v}_{a}/\mathrm{d}t$.[50] Such accelerations are contained in terms at least linear in $G$, so that, replacing the accelerations by the EOM would yield negligible terms of order $G^{2}$ at least. However, it is not allowed to replace accelerations in a Lagrangian while remaining in the same coordinate system. Such replacement is equivalent to a shift in the particles’ trajectories (or so-called “contact” transformation), i.e., the new Lagrangian is physically equivalent to the original one but written in a different coordinate system.[97] Furthermore, by employing the technique of double-zero (or multiple-zero) terms, it is sufficient to consider a Lagrangian that is linear in accelerations (since the procedure can work for any PN order, and is thus formally valid at the 1PM order). Therefore, we look for a Lagrangian of the form

[TABLE]

We symbolize the functional dependence by $L[y,v,a]\equiv L[\{\bm{y}_{a},\bm{v}_{a},\bm{a}_{a}\}]$. The terms $\lambda$ and $q_{a}^{i}$ are of order $G$ and depend only on positions and velocities, i.e., $\lambda[y,v]$ and $q_{a}^{i}[y,v]$. Denoting by $p_{a}^{i}$ and $q_{a}^{i}$ the conjugate momenta associated with the positions $y_{a}^{i}$ and velocities $v_{a}^{i}$, i.e.,

[TABLE]

the EOM take the ordinary Euler-Lagrange form

[TABLE]

while the conserved energy $E$ is given by the generalized Legendre transformation

[TABLE]

In both (40) and (41) we are allowed to replace the accelerations by the EOM. For instance, the term $q_{a}^{i}a_{a}^{i}$ in $E$ will be second-order in $G$ and can be neglected at 1PM order. With (38) we obtain the EOM

[TABLE]

where $f_{a}^{i}=m_{a}\frac{\mathrm{d}}{\mathrm{d}t}(\gamma_{a}v_{a}^{i})$ and $\frac{\delta\lambda}{\delta y_{a}^{i}}=\frac{\partial\lambda}{\partial y_{a}^{i}}-\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\partial\lambda}{\partial v_{a}^{i}})$, and the dots refer to time derivatives. The potential $V$ (such that $E=\sum_{a}m_{a}c^{2}\gamma_{a}+V$) reads

[TABLE]

The left-hand sides of (42) and (43) are known from Eqs. (35)–(5) and (37). However, the two equations are not independent, since $f_{a}^{i}$ and $V$ satisfy the constraint

[TABLE]

In order to find $L$, our strategy is to determine first a particular Lagrangian $\hat{L}$, characterized by $(\hat{\lambda},\hat{q}_{a}^{i})$, such that the conjugate momenta $\hat{q}_{a}^{i}$ obey

[TABLE]

To order $G$, the same equation is also satisfied by the time derivative $\dot{\hat{q}}_{a}^{i}$. Therefore, for the particular solution $(\hat{\lambda},\hat{q}_{a}^{i})$, the equation (43) reduces to an ordinary Legendre transformation,

[TABLE]

To determine $\hat{\lambda}$, we note that the potential $V$ given by (41) reduces in the limit $c\to+\infty$ to the Newtonian approximation, namely $V=U+\mathcal{O}(1/c^{2})$ where

[TABLE]

If we subtract its Newtonian limit $U$ to $V$, we get a quantity which tends to zero when $c\to+\infty$ like $\mathcal{O}(1/c^{2})$. Then, it is straightforward to show that a well-behaved solution of Eq. (46) is

[TABLE]

Namely, we have to insert into Eq. (37) all the relevant factors $c$ and make the replacement of $c$ by $s$, then integrate over the “speed of light” $s$ from the physical value $c$ up to infinity. The bound $s\to+\infty$ of the integral corresponds to the Newtonian limit and we see from the definition of the Newtonian potential (47) that the integral is convergent. The first term in Eq. (48) represents the Newtonian approximation with the correct minus sign for a Lagrangian, and the integral represents formally the complete series of PN corrections, but resummed in the PM approximation. The result (48) can be rewritten in a simpler way as the Hadamard “partie finie” (Pf) of the integral, in the same sense as was used in Eq. (8), for taking care of the divergence at infinity:

[TABLE]

For this very simple type of divergence $\sim s^{0}+\mathcal{O}(s^{-2})$ the Pf does not depend on an arbitrary constant, unlike in (8). The expressions (48)–(49) give a particular solution of the equation (43) but we still have to adjust $\hat{q}_{a}^{i}$ in order to satisfy the EOM, see (42). Thus, we look for $\hat{q}_{a}^{i}$ satisfying

[TABLE]

where the right-hand side is known. To order $G$ we have been able to integrate twice this relation to determine $\hat{q}_{a}^{i}$, and that solution automatically satisfies the constraint (45) by virtue of (44).

Finally we have found a particular Lagrangian $(\hat{\lambda},\hat{q}_{a}^{i})$. Now the general solution $(\lambda,q_{a}^{i})$ can be obtained by adding an arbitrary total time-derivative $\mathrm{d}F/\mathrm{d}t$, where $F$ is a function of the positions $y_{a}^{i}$ and velocities $v_{a}^{i}$. Hence the general solution (for the class of harmonic-coordinate Lagrangians that are linear in accelerations) reads

[TABLE]

At 1PM order the Lagrangian in harmonic coordinates irreducibly depends on accelerations, i.e., it is impossible to determine $F$ such that $q_{a}^{i}=0$. However, we know that the accelerations in a Lagrangian can be eliminated by appropriate shifts of the trajectories. In fact, it can be shown that the particular solution $\hat{\lambda}$ found in (48)–(49) represents an ordinary Lagrangian which is physically equivalent but expressed in some shifted (non harmonic) variables.[36]

Given the complicated structure of $V$ in Eq. (37), we could not find a closed form expression for the 1PM harmonic coordinate Lagrangian in the general case. However, we could easily work out the integral (48) in the PN approximation $c\to+\infty$ to any order. We start from the known 4PN expansion of the potential $V$ following from (37), and explicitly perform the integration (48) term by term, to obtain the corresponding 4PN expansion of $\hat{\lambda}$. Then, we derive the coefficient of accelerations $\hat{q}_{a}^{i}$ at 4PN order, see[36] for details. Finally, we find a unique total time-derivative, with some function $F_{\text{PN}}$ given in the form of a PN expansion, so that the Lagrangian satisfyingly agrees up to order $G$ with the published 4PN Lagrangian.[8, 10] Furthermore, we have pushed the analysis to the next order and obtained all the terms of order $G$ in the harmonic coordinates Lagrangian up to the 5PN order.[36]

In another application, we worked out the case of equal masses for which it is possible to find a closed form expression for the Lagrangian, and we have verified that the associated Hamiltonian differs from the one obtained by Ledvinka, Schäfer and Bičák [78] in the ADM Hamiltonian by a mere canonical transformation.