About gravitational-wave generation by a three-body system

While the two-body problem is completely solvable in Newtonian theory, no general exact solution to it is known in general relativity (GR). As a result, the dynamics of a binary system in GR can only be obtained by solving perturbatively the field equations, or through numerical techniques on a computer ('numerical relativity' [1–5]). Perturbative schemes, valid in different regimes, include the post-Newtonian (PN) approximation [6], which consists in expanding the dynamics in powers of $v/c \ll 1$ ⁵ ($v$ being the relative velocity of the binary, c the speed of light in vacuum) and the self-force formalism, which instead relies upon an expansion in the binary's mass ratio, assumed to be small [7]. The detection of gravitational waves (GWs), indirectly from binary pulsars [8] and directly from systems of two merging black holes (BHs) [9–11], provides an excellent benchmark to test the general-relativistic two-body dynamics. As a matter of fact, these observations are to date in perfect agreement with the GR predictions [10–15].

Unsurprisingly, the three-body problem, which already in Newtonian theory does not admit any tractable closed form solution and, moreover, gives rise to chaotic dynamics, becomes considerably harder in GR. However, its general-relativistic dynamics can be obtained within the PN approximation scheme, just like in the two-body case. Indeed, one can write a PN-expanded, time-dependent Hamiltonian [16–19] that describes the conservative dynamics of a system of three non-spinning bodies up to the 2PN order (i.e. through order $(v/c){\hspace{0pt}}^4$ beyond the leading order Newtonian dynamics), as well as its dissipative dynamics (i.e. the back-reaction due to GW emission) to leading order in $v/c$, which corresponds to a contribution of 2.5PN order, or ${{\mathcal O}}(v/c){\hspace{0pt}}^5$, in the equations of motion.

The PN dynamics of three body systems is not just an academic curiosity. Kozai–Lidov resonances [20–23], first discovered in Newtonian triplets, are believed to be a relevant astrophysical mechanism for the formation of binary systems of stellar-mass BHs (observable by ground-based GW detectors) in dense stellar environments [24, 25] or in isolation [26]. It may also play an important role in the formation and evolution of binaries of supermassive BHs [27–31], whose GW signal is targeted by existing pulsar-timing arrays [32–41] and by the future space-borne interferometer LISA [42]. It turns out that the PN corrections are crucial to assess the efficiency of the Kozai–Lidov mechanism, as they can destroy the resonance on which it relies. Indeed, the coherent piling up of the perturbation induced by the third body may be disrupted due to relativistic precession effects appearing at 1PN order and beyond [22, 27, 43]. Finally, the GW emission from systems of three BHs has been studied in detail by means of numerical techniques [19, 31, 44] in the event that they form and radiate in sufficiently large number to provide a sizable population for GW detectors. For particular configurations, they have also been investigated analytically [47, 48] so as to gain some insight on their dynamics.

The GW emission from binary systems with relative velocities $v\ll c$ can be modeled, at leading order, through the Einstein quadrupole formula [6, 49, 50]. Next-to-leading order corrections are given by the mass-octupole and current-quadrupole contributions [6, 51]. A key requirement implicit in the derivation of the corresponding formulae is that the binary must be contained in its 'near coordinate zone' (NCZ), i.e. a region (centered on the origin of the coordinates) of radius comparable to (but smaller than) the GW wavelength λ. This requirement comes about because the PN formalism for GW generation, which can only be legitimately applied as long as the source is much smaller than λ, is based on a systematic multipole expansion of the gravitational field outside the source. In order to ensure an overlap between the domain of validity of this expansion (say $\vert {\boldsymbol x}\vert \gtrsim r_{\rm min}$) and the near zone (where the dynamics of the source is computed neglecting retardation effects), one must clearly have $r_{\rm min}\sim$ (size of the near zone) $\sim\lambda$, so the coordinate origin and the source cannot be more than one wavelength apart⁶. Indeed, these formulae are usually applied in the reference frame of the binary's CoM. In that frame, in the PN regime, the existence of a NCZ containing the binary is guaranteed, since the size of the system—its separation a—is negligible relative to the wavelength $\lambda\sim a/(v/c)$.

For a triple system with relative velocities $v \ll c$, it would seem natural to apply the very same formulae in the reference frame of the CoM of the three-body system. However, by doing so, one obtains unphysical results such as those reported in figure 18 of [19], as we will now explain. Indeed, we have reproduced the same behavior by applying the quadrupole and 'quadrupole–octupole' formulae in the CoM reference frame of a series of triple systems with mass ratios $m_2/m_1=0.5$ and $m_3/(m_1+m_2)=0.05$, whose trajectories are computed with the code of [52] (which includes the 1PN and 2PN conservative triple dynamics, and the leading order dissipative dynamics). The 'inner binary' (comprised of m₁ and m₂) of these hierarchical triplets has zero initial eccentricity and an initial separation $a_{\rm in}=150Gm_t/c^{2}$, where G is Newton's constant and m_t the total mass of the triplet (throughout the paper we instead reserve the symbol m to indicate the total mass of binary systems, i.e. $m=m_1+m_2$). The 'outer binary' (comprised of m₃ and the CoM of the inner binary) has instead initial separation varying in the range $a_{\rm out} \in [625, 10\, 000]Gm_t/c^{2}$, and zero initial eccentricity. The results are displayed in figure 1, where one can observe, paradoxically, that the effect induced by the third body grows as it gets farther away from the inner binary. We will analyze this situation in detail in this paper, and show that the problem is connected to the fact that a NCZ region centered on the CoM of the triplet and having size comparable to the minimum gravitational wavelength excited by the system does not include the whole triplet, unlike what happens for a binary system.

Zoom In Zoom Out Reset image size

This, however, is just one example of the subtleties one should be aware of when computing GW emission from binary or triple systems in a too naive fashion. Another interesting apparent paradox arises, e.g. if one tries to compute the gravitational waveforms of a binary (or triple) system by directly integrating the equations for the linear perturbations $h_{\mu\nu}$ over a background Minkowski space-time (endowed with a flat metric $\newcommand{\e}{{\rm e}} \eta_{\mu\nu}={\rm diag} (-1, 1, 1, 1)$ and coordinates $\{x^\mu\}_{\mu=0, 1, 2, 3}$).

In the harmonic gauge, which is defined by the condition $\partial_\mu \bar{h}^{\mu\nu}=0$, where $\partial_\mu$ is the flat four-dimensional derivative and $\newcommand{\e}{{\rm e}} \bar{h}^{\mu\nu} = h^{\mu\nu} - 1/2 \ \eta^{\mu\nu} h^{\alpha}_{~\alpha}$ represents the trace-reversed metric perturbation⁷, the linearized Einstein equations read (see, e.g. [53] section 1.1, [54] section 35.1)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:first_wave_eq} \Box_{\rm flat} \bar{h}^{\mu\nu}= -\frac{16 \pi G}{c^4} T^{\mu\nu}, \nonumber \end{align} \tag{ 1 }$$

where the d'Alembert operator $\newcommand{\e}{{\rm e}} \Box_{\rm flat}= \eta^{\mu\nu}\partial_\mu\partial_\nu$ is computed with the background Minkowski metric and $T^{\mu\nu}$ is the source stress-energy tensor. These equations can be integrated exactly by using the (retarded) Green function of $\Box_{\rm flat}$. The resulting waveforms (obtained from the transverse trace-free part of the spatial components) may then be compared to those predicted by the quadrupole formula (and its higher-order corrections that we have mentioned above).

The comparison between the GW amplitudes obtained with the two procedures for various binaries is shown in figure 2. As can be seen, there appears to be a factor  ∼2 discrepancy (this factor becomes exactly 2 for binary circular orbits). Similar discrepancies arise when integrating equation (1) for triple systems. This puzzling difference will be discussed in more details. It is related to the fact, often mentioned but rarely illustrated in introductory GR textbooks (see however [53, 54]), that a naive derivation of the quadrupole formula based on equation (1) is wrong. It is because that equation (via the harmonic gauge condition) implies that $\partial_\mu T^{\mu\nu}=0$, which is clearly not verified for a binary system since it entails that bodies move along straight lines.

Zoom In Zoom Out Reset image size

The focus of this paper is thus pedagogical. We will discuss the two problems mentioned above as well as other subtleties that we have encountered when computing gravitational waveforms from numerically-integrated orbits of triple systems. More precisely, the organization is as follows: in section 2 we will illustrate, tackle and solve the problems that arise when applying the standard quadrupole formula (and its higher-order corrections) to a triple system. This will provide a solution to the discrepancy demonstrated in figure 1, which, as already mentioned, will turn out to be due to the source not being contained in its NCZ when the latter is centered on the triplet's center of mass. In section 3 we will further comment about the inconsistent derivation of equation (1), highlighting the need to use the non-linear Einstein equations to compute GW emission self-consistently. Finally, in section 4 we will draw our conclusions. Throughout the paper we use the $(-, +, +, +)$ signature convention.

The leading-order contribution to the GW signal observed at space position ${\boldsymbol x}$ and time t is given, in an appropriate 'radiative' gauge, by the quadrupole formula (see [54] section 36.10, [53] section 3.3, [6] section 2.5)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:first_quad_formula} h_{ij}^{\rm TT}(t,{\boldsymbol x}) = \frac{2 G}{R c^4} \Lambda_{ijkl}({\boldsymbol n}) \ \frac{\ud^2 M^{kl}(t_{\rm ret})}{\ud t^2} + \mathcal{O}\Big(\frac{1}{c^5}\Big)\,, \nonumber \end{align} \tag{ 2 }$$

where $\newcommand{\e}{{\rm e}} R=\vert {\boldsymbol x}\vert \equiv \sqrt{x^i x_i}$ is the distance of the observer (assumed to be very far from the source compared with the wavelength λ of the emitted GWs), $t_{\rm ret}=t-R/c$ is the retarded time of the background space-time⁸, $\Lambda_{ijkl}({\boldsymbol n})$ denotes the projector on the transverse-traceless (TT) gauge (see appendix A for the explicit definition), while

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{Q} M^{ij}(t) = \int {\ud^3{\boldsymbol x}' \ c^{-2}T^{00}(t,{\boldsymbol x}^{\prime}) \Big(x^{{\prime i}}{x^{\prime j}}-\frac{1}{3}\delta^{ij}x^{{\prime}k} {x^{\prime}}\!_k \Big)} \nonumber \end{align} \tag{ 3 }$$

represents the mass quadrupole moment of the source. As mentioned in the introduction, implicit in the derivation of the quadrupole formula (equation (2)) is the assumption that the source be contained in its NCZ (see, e.g. [53, 54, 57]), i.e. the reference frame in which the quadrupole moment (equation (3)) is evaluated must be such that the source be contained within a region of size  ∼λ centered on the origin of the coordinates. Finding a frame satisfying this property is always possible for slowly moving binary systems, since λ is related to the system's typical (relative) velocity $v$ and its typical separation a by $\lambda\sim a/(v/c)$.

The most natural reference frame to describe the dynamics of an N-body system is that of the CoM, where the equations of motion take their simplest form. (Note that the usual Newtonian expression of the CoM position in terms of the body locations ${\boldsymbol x}_1$, ${\boldsymbol x}_2$ for a binary, namely ${\boldsymbol r}_0=(m_1 {\boldsymbol x}_1 + m_2 {\boldsymbol x}_2)/m$, is modified beyond the leading order⁹.) In the case under consideration here (i.e. a three-body system), a particularly interesting configuration is that of the so-called hierarchical triplet. The latter is comprised of an inner close binary $(m_1, m_2)$, supplemented by a third body m₃ at larger distance. In practice, the system in this configuration can be regarded as formed by two separate binaries: an inner binary on the one hand, and an outer binary, formed by the third outer body and the CoM of the inner binary, on the other hand. From the point of view of the dynamics, the choice of the reference frame is of course irrelevant. Indeed, the Hamiltonian (both at the Newtonian order and when including the PN corrections) depends only upon the relative separations of the three bodies, hence the dynamics of the system is frame-independent (see e.g. [60]). However, caution must be exercised when applying the quadrupole formula (equation (2)) (and its higher-order generalizations including octupolar corrections, etc) to the orbits resulting from numerical integrations of the equations of motion (see figure 1 and related discussion).

To illustrate this point, let us consider a circular binary located 'far' from the origin of the coordinates. In triple systems, this happens for the inner binary when $m_3 a_{\rm out} \gg (m_1+m_2+m_3) a_{\rm in}$ since, in this case, the inner binary is located far away from the CoM of the triplet chosen as the origin of the coordinates. Setting the GW source in the xy plane and the observer along the z axis, equation (2) takes the simplest possible form, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:matrix_pol} h_{ij}^{\rm TT}(t,{\boldsymbol x}) = \frac{2 G}{R c^4} \left(\begin{array}{@{}ccc@{}} \frac{1}{2}[\ddot{M}_{11}(t_{\rm ret}) - \ddot{M}_{22}(t_{\rm ret})] & \ddot{M}_{12}(t_{\rm ret}) & 0 \nonumber \\ \ddot{M}_{12}(t_{\rm ret}) & -\frac{1}{2}[\ddot{M}_{11}(t_{\rm ret}) - \ddot{M}_{22}(t_{\rm ret})] & 0 \nonumber \\ 0 & 0 & 0 \end{array}\right)\,, \nonumber \end{align} \tag{ 4 }$$

where the dots placed over the quadrupole components $M_{ij}$ represent time derivatives. The two independent polarizations of a propagating GW, referred to as the 'plus' and 'cross' polarizations, are (in the situation considered here) simply the diagonal and off-diagonal part of equation (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (h_+)_{\rm quad} &= \frac{G}{R c^4}(\ddot{M}_{11} - \ddot{M}_{22})\nonumber \\ (h_{\times})_{\rm quad} &= \frac{2 G}{R c^4} \ddot{M}_{12}\,, \label{eq:matrix_pol2} \nonumber \end{align} \tag{ 5 }$$

where all $M_{ij}$ components are evaluated at retarded time. Explicitly, for a binary these expressions become (see, e.g. Problem 3.2 of [53])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_circ_bin} (h_+)_{\rm quad} &= \frac{4G^2 m_1 m_2}{a R c^4}\cos(2\omega t_{\rm ret})\nonumber \\ (h_{\times})_{\rm quad} &= \frac{4G^2 m_1 m_2}{a R c^4}\sin(2\omega t_{\rm ret})\,, \nonumber \end{align} \tag{ 6 }$$

a denoting the separation of the binary and ω its orbital frequency.

Let us now consider two different circular binaries, both with $m_1 = 0.9~m$ and $m_2 = 0.1~m$ but representative of two different regimes: a rather relativistic binary with separation $a = 20Gm/c^{2}$, which corresponds to a relative orbital velocity $v/c \simeq 0.2$, and a mildly-relativistic one, with $a = 2000Gm/c^{2}$ (corresponding to $v/c\simeq 0.02$). For these two systems, equations (6) gives $Rc^2/({Gm})(h_+)_{\rm quad} \simeq 1.8\times 10^{-2}$ and $Rc^2/({Gm})(h_+)_{\rm quad} \simeq 1.8\times 10^{-4}$, respectively. We then evolve them numerically in two different frames: (1) one with the origin coinciding with the CoM, and (2) one with the origin shifted by 10⁵ gravitational radii (i.e. $10^5\times Gm/c^{2}$) from the CoM. Next, we compute the waveforms directly via equations (5) from the numerical trajectories. Results are reported in figure 3 as dashed lines. The mildly-relativistic case is consistent with the analytic predictions of equations (6). Instead, for the relativistic binary, $(h_{+})_{\rm quad}$ given by equations (5) is more than one order of magnitude higher than the prediction from equations (6) when the origin of coordinates is far away from the CoM. Indeed, we have checked that applying equations (5) directly to numerically-integrated trajectories yields results that are coordinate-dependent. The discrepancy with equations (6) grows with the binary's relative velocity.

Zoom In Zoom Out Reset image size

Let us consider now the next-to-leading order contributions to the waveform, comprised of a mass octupole and a current quadrupole term. When these terms are taken into account, equation (2) becomes (see [53] section 3.4)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:h_oct} h_{ij}^{\rm TT}(t,{\boldsymbol x}) = \frac{2 G}{R c^4}\Lambda_{ijkl}({\boldsymbol n}) \biggl[ \frac{\ud^2 M^{kl}}{\ud t^2} + \frac{n_m}{3c}\Big(\frac{\ud^3 M^{klm}}{\ud t^3} + 2 \frac{\ud^2 S^{klm}}{\ud t^2}\Big) \biggr] + \mathcal{O}\Big(\frac{1}{c^6}\Big)\,, \nonumber \end{align} \tag{ 7 }$$

where ${\boldsymbol n} = {\boldsymbol x}/R$, and $M^{klm}, S^{klm}$ are respectively the Newtonian octupole and current quadrupole moments (evaluated at t_ret (see appendix A)). Again, for a binary in the xy plane and an observer in the same plane¹⁰, the above equation implies (see problem 3.3 of [53] and the corresponding erratum)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_oct_bin} (h_+)_{\rm oc+cq} &= \frac{1}{R} \frac{G^{5/2} m_1 m_2 \delta m}{4 (m a^3)^{1/2} c^5} \left[5\cos(\omega t_{\rm ret}) - 9\cos(3\omega t_{\rm ret}) \right] \nonumber \\ (h_{\times})_{\rm oc+cq} &= 0\,, \nonumber \end{align} \tag{ 8 }$$

with $\delta m = m_1-m_2$. For the two binaries considered above, we obtain $Rc^2/({Gm}) (h_+)_{\rm oc+cq}\simeq$$2.2\times 10^{-3}$ in the relativistic case, and $Rc^2/({Gm}) (h_+)_{\rm oc+cq} \simeq 2.2\times 10^{-6}$ in the mildly relativistic one. From figure 4 (dashed lines), we can see that if the origin of the coordinates coincides with the CoM, then $(h_+)_{\rm oc+cq}$ as given by equation (7) applied to the numerically-integrated orbits agrees well with the analytic result of equations (8). Conversely, if one shifts significantly the origin of the coordinates, $(h_+)_{\rm oc+cq}$ computed from the numerical trajectories no longer agrees with equations (8), even for the mildly-relativistic binary.

Zoom In Zoom Out Reset image size

In the next two sections, we will analyze the reasons behind these discrepancies and explain how they can be avoided, first for the Newtonian quadrupole formula (section 2.1), next for the 0.5PN quadrupole formula with octupolar corrections (section 2.2).

2.1. Quadrupole waveform

If we explicitly compute the second-order time derivatives in equation (2), for a binary system, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_der_exp} h_{ij}^{\rm TT}(t,{\boldsymbol x}) =&\, \frac{2 G}{R c^4} \Lambda_{ijkl}({\boldsymbol n})\biggl(2 m_1 \dot{x}_1^k \dot{x}_1^l + 2 m_2 \dot{x}_2^k \dot{x}_2^l \nonumber \\ &+ m_1 (\ddot{x}_1^k x_1^l + x_1^k \ddot{x}_1^l) + m_2 (\ddot{x}_2^k x_2^l + x_2^k \ddot{x}_2^l) \biggr) + \mathcal{O}\Big(\frac{1}{c^5}\Big)\, . \nonumber \end{align} \tag{ 9 }$$

Now, the position vectors of the two masses can be expressed in terms of the CoM position ${\boldsymbol r}_{0}$ and the relative separation vector ${\boldsymbol r}={\boldsymbol x}_1-{\boldsymbol x}_2$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:cm_def} {\boldsymbol x}_1 = {\boldsymbol r}_{0} + \frac{m_2}{m}{\boldsymbol r}\,, \qquad \qquad {\boldsymbol x}_2 = {\boldsymbol r}_{0} - \frac{m_1}{m}{\boldsymbol r}, \nonumber \end{align} \tag{ 10 }$$

so that equation (9) takes the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_cm_sub} h_{ij}^{\rm TT}(t,{\boldsymbol x}) =&\, \frac{2 G}{R c^4} \Lambda_{ijkl}({\boldsymbol n}) \bigg[ \frac{2}{m} (m_1 \dot{x}_1^k+m_2 \dot{x}_2^k)(m_1 \dot{x}_1^l+m_2 \dot{x}_2^l) \nonumber \\ & + \frac{2m_1 m_2}{m} (\dot{x}_1^k-\dot{x}_2^k)(\dot{x}_1^l-\dot{x}_2^l) + r_0^{k}(m_1 \ddot{x}_1^{l} + m_2 \ddot{x}_2^{l}) + r_0^{l}(m_1 \ddot{x}_1^{k} + m_2 \ddot{x}_2^{k})\nonumber \\ & + \frac{m_1 m_2}{m} \Big(r^k(\ddot{x}_1^l-\ddot{x}_2^l) + r^l(\ddot{x}_1^k-\ddot{x}_2^k)\Big) \bigg] + \mathcal{O}\Big(\frac{1}{c^5}\Big)\, . \nonumber \end{align} \tag{ 11 }$$

Since the CoM absolute coordinates $\boldsymbol r_0$ explicitly appears in this expression, it would seem that the GW amplitude should depend on the choice of the origin of the coordinate system. However, because equation (11) is only correct at leading order in PN theory, it is actually sufficient to compute the accelerations $\ddot{\boldsymbol{x}}_1$ and $\ddot{\boldsymbol{x}}_2$ at leading (i.e. Newtonian) order. If one does so, the identity $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2=0$ holds for an isolated system, hence the dependence on the position of the center of mass (and thus on the location of the origin) disappears from equation (11). Similarly, $m_1 \dot{x}_1^k+m_2 \dot{x}_2^k$ is constant and independent of the location of the origin.

One may want, however, to integrate the binary's equations of motion to higher PN order, either analytically or numerically. For instance, the code of [52] integrates the PN Hamiltonian for binary or triple systems through the 2.5PN order including the dissipative effects of radiation reaction. Now, when one includes these PN corrections, $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2\neq0$ already at 1PN order, so that the dependence on the location of the origin does not disappear. This is the reason of the unphysical behavior visible in figure 3 (and partly in figure 1, see also the next section).

A first solution can be therefore to avoid the use of the numerical trajectories to compute the accelerations in equation (11), but instead evaluate them directly from the positions of the two bodies by using the Newtonian dynamics (i.e. Newton's second law, which ensures $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2=0$). Alternatively, one can note that the combination $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2$ is simply (at Newtonian order) the time derivative of the total linear momentum ${\boldsymbol P}_{\rm N}=m_1 \dot{\boldsymbol{x}}_1 + m_2 \dot{\boldsymbol{x}}_2$. Thus, the identity $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2=0$ just reflects the conservation of ${\boldsymbol P}_{\rm N}$, which is an automatic consequence of the Newtonian dynamics. Beyond it, when PN corrections are included, $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2$ does not vanish, because the Newtonian linear momentum ${\boldsymbol P}_{\rm N}=m_1 \dot{\boldsymbol{x}}_1 + m_2 \dot{\boldsymbol{x}}_2$ is no longer a conserved quantity. This is what causes the dependence on the choice of ${\boldsymbol r}_0$ observed in figure 3. However, one can exploit the fact that there exists a conserved PN linear momentum ${\boldsymbol P}_{n {\rm PN}}$ generalizing ${\boldsymbol P}_{\rm N}$ at the nPN order. In practice, replacing ${\boldsymbol P}_{\rm N}$ with ${\boldsymbol P}_{n {\rm PN}}$ is equivalent to computing the accelerations appearing in equation (11) as $\ddot{\boldsymbol{x}}_i={\boldsymbol \pi}_i/m_i$, with $i=1, 2$ and ${\boldsymbol \pi}_i$ denoting the conjugate momentum of each body entering the Hamilton equations. Then, the combination $m_1 \ddot{\boldsymbol{x}}_1 + m_2 \ddot{\boldsymbol{x}}_2$ always vanishes, even if PN corrections are included in the Hamiltonian dynamics.

In conclusion, either of these two workarounds (which give rise to the 'amended' waveforms represented by blue solid lines in figure 3) is sufficient to eliminate the unphysical dependence on the origin of the coordinates.

2.2. Octupole and Current Quadrupole waveforms

Let us now examine what happens to the contribution of the mass octupole and current quadrupole moments to the waveform under a change of reference frame of the form of equations (10). After expanding the time derivatives appearing in the term $n_m/(3c)\, (\mathop{M}\limits^{\ldots}{\hspace{0pt}}^{klm} + 2 \ddot{S}^{klm})$ of equation (7) by means of the Leibniz rule, the dependence on the CoM location does not cancel out in the waveform, but gives instead a contribution¹¹

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:offending_terms} \delta_{r_0} h_{ij}^{\rm TT} = \Lambda_{ijkl}({\boldsymbol n})\frac{n_m r_0^m}{3c}\sum_A m_A \left(3 \mathop{x}\limits^{\ldots}{\hspace{0pt}}\!_A^k x_A^l + 3 x_A^k \mathop{x}\limits^{\ldots}{\hspace{0pt}}\!_A^l + 9 \ddot{x}_A^k \dot{x}_A^l + 9 \dot{x}_A^k \ddot{x}_A^l\right)\, . \nonumber \end{align} \tag{ 12 }$$

This may be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:d3Q} \delta_{r_0} h_{ij}^{\rm TT}=\Lambda_{ijkl}({\boldsymbol n})\frac{\ud^3 M^{kl}}{\ud t^3}\frac{\delta R}{c}\,, \nonumber \end{align} \tag{ 13 }$$

with $\delta R = {\boldsymbol n} \cdot {\boldsymbol r}_0 = n_m r_0^m$. Hence, the terms proportional to ${\boldsymbol r}_0$ produced by the mass octupole and current quadrupole moments can be reabsorbed in a shift $\delta R/c$ of the retarded time at which the (quadrupole) waveform is evaluated:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:resum} \Lambda_{ijkl}({\boldsymbol n}) \bigg[ \frac{\ud^2 M^{kl}(t-R/c)}{\ud t^2} &+ \frac{\ud^3 M^{kl}(t-R/c)}{\ud t^3} \frac{\delta R}{c}\bigg] \nonumber \\ &= \Lambda_{ijkl}({\boldsymbol n}) \bigg[\frac{\ud^2 M^{kl}(t-R/c+\delta R/c)}{\ud t^2} \bigg] + \mathcal{O}\bigg(\frac{[\delta R]^2}{c}\bigg) \, . \nonumber \end{align} \tag{ 14 }$$

This time shift simply enforces the invariance of the waveform under translations of the reference frame in which it is computed. Indeed, the retarded time is always, by definition, $t_{\rm ret} = t-\vert \mathbf{x}\vert /c$ in the generic frame $(t, {\boldsymbol x})$ (assuming radiative coordinates), where $R=\vert {\boldsymbol x}\vert$ is the distance between the observer at position ${\boldsymbol x}$ and the origin. This implies in particular that $t_{\rm ret}-t_{\rm ret}^{\rm CoM} = \vert {\boldsymbol x}^{\rm CoM}\vert /c-\vert {\boldsymbol x}\vert /c = -\delta R/c + \mathcal{O}([\delta R]^2/c)$, where the quantities labeled with the superscript CoM are referring to the CoM frame. On the other hand, the expressions of the multipole moments when ${\boldsymbol r}_0\neq {\boldsymbol 0}$ differ from their standard forms in CoM coordinates. Obviously, the modifications of t_ret and those of the multipole moments must (and do!) compensate each other so that the GW signal remains invariant, irrespective of the choice of the origin of the coordinates. This main conclusion remains true for linearly propagating waves when other multipole moments are taken into account¹²:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:hinvariance} h_{ij}^{\rm TT} [\{M_{ab},M_{abc},...,S_{abc},...\},t,{\boldsymbol x}] = h_{ij}^{\rm TT} [\{M_{ab}^{\rm CoM},M_{abc}^{\rm CoM},...,S_{abc}^{\rm CoM},...\},t,{\boldsymbol x}^{\rm CoM}]\, . \nonumber \end{align} \tag{ 15 }$$

However, the analytic 'resummations' needed for the above argument to work, such as the one in equation (14), are based on a Taylor expansion. Therefore, one has to implicitly assume that terms like those in equation (12) are 'small' or, more precisely, that the displacement $\vert {\boldsymbol r}_0\vert$ of the CoM from the origin of the coordinates is much smaller than the wavelength of the quadrupole waveform, $\lambda=\pi c/\omega$. This is in fact very natural because, as already mentioned, one of the assumptions implicit in the derivation of the quadrupole/octupole formulae is that the source be well contained in a NCZ of size  ∼λ centered on the origin of the coordinates, where retardation effects are negligible. If the generalized quadrupole formula is applied to systems for which the source is not well contained in its NCZ, one will not be able to resum the terms of equation (12) into a time shift. This is the origin of the discrepancy shown in figure 4 for binary systems.

This observation also highlights the reason of the unphysical behavior shown in figure 5 (left panels) for triple systems. In fact, for a weakly/mildly relativistic binary ($v\lesssim c$) with separation a, one has $\lambda \sim a /(v/c) \gtrsim a$, i.e. if one chooses the origin of the coordinates to coincide with the binary's CoM, the NCZ will always contain the binary. For a hierarchical triple system, instead, if one sets the origin at the location of the triplet's CoM, the inner binary will be outside its NCZ provided that the separation of the outer binary is sufficiently large. This may be understood by noting that there are actually two NCZs for a hierarchical triplet, i.e. an inner-binary NCZ with size $\lambda_{\rm in}\sim a_{\rm in}/(v_{\rm in}/c)$, and an outer-binary NCZ with size $\lambda_{\rm out}\sim a_{\rm out}/(v_{\rm out}/c)$, as illustrated in figure 6. Clearly, while the outer binary will always be contained in its NCZ if $v_{\rm out}/c\lesssim 1$, the inner binary will eventually be outside its NCZ if a_out is sufficiently large. Indeed, as a_out increases, the inner binary's CoM ends up leaving its own NCZ (which is centered on the origin of the coordinates, i.e. on the triplet's CoM).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A simple fix to this issue, as shown in figure 6, is thus to evaluate the multipole moments in an inertial reference frame with origin instantaneously coinciding with the CoM of the inner binary¹³, which allows both the inner and outer binaries to lie within their respective NCZs. In the right panels of figure 5, we show that this eliminates the unphysical behavior of figure 1 and left panels of figure 5. Therefore, the problem exhibited by those figures (and the corresponding results of figure 18 in [19]) was simply that the waveforms were evaluated in the reference frame of the CoM of the triple system. The very same solution applies to the simpler binary cases reported in figure 4. Here, since the CoM does not move, only a single transformation is needed. Once that transformation has been performed, the waveforms are given by the solid blue lines and reproduce the correct predicted result.

As pointed out in the introduction, another often overlooked problem arises when gravitational waveforms are computed by direct integration of equation (1) with the help of the retarded Green function. To understand it, let us go back to the textbook derivation of the quadrupole formula for GW generation. We start from the Einstein equations relaxed by the condition of harmonic coordinates (see, e.g. [54] section 20.3, section 36.9, or [6] section 6.3 and [61] for the complete derivation):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:relaxed_eq} \Box_{\rm flat} H^{\alpha\beta} = -\frac{16 \pi G}{c^4}\tau^{\alpha\beta}\,, \nonumber \end{align} \tag{ 16 }$$

where the pseudo-tensor $H^{\alpha\beta}$ is defined in terms of the Minkowski metric $\newcommand{\e}{{\rm e}} \eta^{\alpha\beta}$ and the space-time (inverse) metric $g^{\alpha\beta}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def_H} H^{\alpha\beta} \equiv \eta^{\alpha\beta} - (-g)^{1/2}g^{\alpha\beta}, \nonumber \end{align} \tag{ 17 }$$

and satisfies the harmonic gauge condition, i.e. $\partial_{\beta} H^{\alpha\beta} = 0$. Moreover, the 'effective' stress-energy pseudo-tensor $\tau^{\alpha\beta}$ is comprised of a contribution from the stress-energy tensor of matter, and a contribution $\Lambda^{\alpha\beta}$ from the non-linearities of the gravitational field, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def_tau} \tau^{\alpha\beta} = (-g)T^{\alpha\beta} + \frac{c^4}{16\pi G} \Lambda^{\alpha\beta}\,, \nonumber \end{align} \tag{ 18 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda^{\alpha \beta} = \frac{16\pi G}{c^4} (-g) t_{\rm LL}^{\alpha \beta } + (\partial_\nu H^{\alpha \mu} \partial_\mu H^{\beta \nu} -\partial_{\mu}\partial_{\nu} H^{\alpha \beta} H^{\mu \nu})\, \nonumber \end{align} \tag{ 19 }$$

with $t_{\rm LL}^{\alpha \beta }$ denoting here the Landau–Lifshitz pseudo-tensor [50]. As a consequence of the harmonic gauge condition, $\tau^{\alpha\beta}$ is also flat-space conserved, i.e. $\partial_{\beta} \tau^{\alpha\beta} = 0$.

By using the retarded Green function, we can now integrate equation (16) and obtain its formal solution:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:H_local_ret} H^{\alpha\beta}(t,{\boldsymbol x}) = \frac{4G}{c^4}\int \frac{\ud^3{\boldsymbol x}^{\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert }\tau^{\alpha\beta} \left(t-\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert /c,{\boldsymbol x}^{\prime}\right). \nonumber \end{align} \tag{ 20 }$$

If the field point ${\boldsymbol x}$ lies very far from the source, $\newcommand{\e}{{\rm e}} \vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert \approx \vert {\boldsymbol x}\vert \equiv R$, and one can neglect the differences in retarded time among the source components by considering a single global retardation $t_{\rm ret}=t-R/c$. This yields¹⁴

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:H_global_ret} H^{\alpha\beta}(t,{\boldsymbol x}) \approx \frac{4G}{R c^4}\int \ud^3{\boldsymbol x}^{\prime}\tau^{\alpha\beta}\left(t_{\rm ret},{\boldsymbol x}^{\prime}\right)\, . \nonumber \end{align} \tag{ 21 }$$

Let us now expand the metric up to $1/c^2$ corrections,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:phi} g_{00} &= -1 -2\frac{\phi}{c^2} + \mathcal{O}\Big(\frac{1}{c^4}\Big)\nonumber \\ g_{0i} &= \mathcal{O}\Big(\frac{1}{c^3}\Big)\nonumber \\ g_{ij} &= \Big(1-2\frac{\phi}{c^2}\Big)\delta_{ij} + \mathcal{O}\Big(\frac{1}{c^4}\Big)\,, \nonumber \end{align} \tag{ 22 }$$

so that the ij components are given at the 1PN order while the 00 and 0i components are Newtonian. The gravitational potential ϕ in equations (22) must satisfy the Poisson equation $\nabla^2 \phi = 4\pi G \rho$, with ρ being the mass density, in order for ϕ to be a solution of the relaxed Einstein equations. We can see that, at this accuracy level, the metric is linear in the source, which implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H^{\mu\nu} = \bar{h}^{\mu\nu} + \mathcal{O}\Big(\frac{1}{c^4},\frac{1}{c^3},\frac{1}{c^4}\Big)\,, \nonumber \end{align} \tag{ 23 }$$

where the three remainders in the arguments of the Landau symbol refer to the 00, 0i and ij components, respectively. From the flat-space conservation of $\tau^{\mu\nu}$, it then follows that (see [54] section 36.10)¹⁵

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:MTW_identity} \partial_0\partial_0(\tau^{00} x^{\,j} x^k) =\partial_l\partial_m(\tau^{lm} x^{\,j} x^k) - 2\partial_l(\tau^{\,jl} x^k + \tau^{kl} x^{\,j}) + 2 \tau^{\,jk}\,, \nonumber \end{align} \tag{ 24 }$$

which allows one to recast the spatial part of equation (21) as (see [54] section 36.10)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:hij} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{2G}{R c^4} \frac{\ud^2}{\ud t^2}\int \ud^3{\boldsymbol x}^{\prime} c^{-2} \tau^{00}\left(t_{\rm ret},{\boldsymbol x}^{\prime}\right)x^{{\prime}i}{x^{\prime j}} + \mathcal{O}\Big(\frac{1}{c^5}\Big). \nonumber \end{align} \tag{ 25 }$$

A priori, the PN expansion of the stress-energy pseudo-tensor $\tau^{\alpha\beta}$ defined by equation (18) still involves both the matter term $T^{\alpha\beta}$ and the purely gravitational non-linear source term $\Lambda^{\alpha\beta}$. Modeling the matter system by point particles, i.e. taking (see [68] section 2.8, [53] section 3.3.5)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:T_point} T^{\alpha\beta}(t,{\boldsymbol x}) = \sum_A \frac{m_A u_A^{\alpha} u_A^{\beta}}{(u^0_A/c) \ \sqrt{-g}}\delta^3({\boldsymbol x}-{\boldsymbol x}_A(t))\,, \nonumber \end{align} \tag{ 26 }$$

where m_A is the mass of particle A and $u_A^\alpha$ is its four-velocity, one finds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:tau_expansion_a} \tau^{00}(t,{\boldsymbol x}) = \sum_A m_A c^2 \delta^3\!\left({\boldsymbol x}-{\boldsymbol x}_A(t)\right) + \mathcal{O}\Big(\frac{1}{c^0}\Big) \,, \nonumber \end{align} \tag{ 27 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:tau_expansion_b} \tau^{0i}(t,{\boldsymbol x}) = \sum_A m_A c \, \dot{x}_A^i \delta^3\!\left({\boldsymbol x}-{\boldsymbol x}_A(t)\right) + \mathcal{O}\Big(\frac{1}{c}\Big) \,, \nonumber \end{align} \tag{ 28 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:tau_expansion_c} \tau^{ij}(t,{\boldsymbol x}) =&\, \sum_A m_A \dot{x}_A^i \dot{x}_A^{\,j} \delta^3\!\left({\boldsymbol x}-{\boldsymbol x}_A(t)\right) \nonumber \\ & + \frac{1}{4\pi G}\Big(\partial^i\phi\, \partial^{\,j}\phi- \frac{1}{2}\delta^{ij}\partial_k\phi \, \partial^k\phi\Big) + \mathcal{O}\Big(\frac{1}{c^2}\Big)\,, \nonumber \end{align} \tag{ 29 }$$

where $\delta^3$ is the three-dimensional Dirac delta function and $\dot{x}_A^i$ represents the components of the three-dimensional velocity. We observe that $\tau^{00}\approx T^{00}$ at leading PN order, hence equation (25) coincides with the 'usual' quadrupole formula,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:4polo} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{2 G}{R c^4} \frac{\ud^2 M^{ij}}{\ud t^2}(t_{\rm ret}) + \mathcal{O}\Big(\frac{1}{c^5}\Big)\, . \nonumber \end{align} \tag{ 30 }$$

By contrast, $\tau^{ij}$ contains a direct contribution from the gravitational field at Newtonian order, i.e. a term arising at the same order as $T^{ij}$. Therefore, the direct integration of equation (1) will give an incorrect result, because the right-hand side is wrong already at the leading PN order.

To be more explicit, let us evaluate equation (25) in the special case of a binary system ($A, B=1, 2$) by using equation (27) and the Newtonian equations of motion

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ddot{x}_A^i = -\sum_{B \neq A} \frac{G m_B }{r_{AB}^2}n^i_{AB}\,, \nonumber \end{align} \tag{ 31 }$$

where ${\boldsymbol r}_{AB} = {\boldsymbol x}_A - {\boldsymbol x}_B$, $r_{AB} = \vert {\boldsymbol x}_A - {\boldsymbol x}_B\vert$, and ${\boldsymbol n}_{AB} = {\boldsymbol r}_{AB}/r_{AB}$. Straightforward algebra yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_2b_quadrupole} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{4 G}{R c^4}\Big(m_1 \dot{x}_1^i \dot{x}_1^{\,j} + m_2 \dot{x}_2^i \dot{x}_2^{\,j} -\frac{G m_1 m_2 n_{12}^i n_{12}^{\,j}}{r_{12}}\Big) + \mathcal{O}\Big(\frac{1}{c^5}\Big)\,, \nonumber \end{align} \tag{ 32 }$$

while direct integration of equation (1) would lead to the different expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h_green} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{4 G}{R c^4}\Big(m_1 \dot{x}_1^i \dot{x}_1^{\,j} + m_2 \dot{x}_2^i \dot{x}_2^{\,j} \Big) + \mathcal{O}\Big(\frac{1}{c^5}\Big)\, . \nonumber \end{align} \tag{ 33 }$$

The extra term in equation (32) that is missing in equation (33) is related to the purely gravitational part of the right-hand side of equation (29) (i.e. the part involving the Newtonian potential ϕ), which is prematurely neglected in the source $T^{ij}$ of equation (1). This is the origin of the factor  ∼2 discrepancy shown in figure 2.

Nonetheless, one can still obtain the correct expression without resorting to the identity (24). In fact, by substituting equation (29) into the spatial components of equation (21), one gets

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:H_spatial} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{4G}{c^4} \int \frac{\ud^3{\boldsymbol x}^{\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert }\biggl[ T^{ij} + \frac{1}{4\pi G} \Big(\partial^i\phi\, \partial^{\,j}\phi- \frac{1}{2}\delta^{ij}\partial_k\phi\, \partial^k\phi\Big)\biggr] + \mathcal{O}\Big(\frac{1}{c^5}\Big), \nonumber \end{align} \tag{ 34 }$$

which may be rewritten as (see appendix B for details)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:h_cont} \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{4G}{R c^4}\int \ud^3{\boldsymbol x}^{\prime} \ T^{ij}- \frac{2 G^2}{R c^4}\int \ud^3{\boldsymbol y}^{\prime} \ud^3{\boldsymbol y}^{\prime\prime}\rho({\boldsymbol y}^{\prime})\rho({\boldsymbol y}^{\prime\prime}) \frac{\hat{n}^i \hat{n}^{\,j}}{\vert {\boldsymbol y}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert } + \mathcal{O}\Big(\frac{1}{c^5}\Big)\,, \nonumber \end{align} \tag{ 35 }$$

where $\hat{n}^{k} = (y^{{\prime}{k}}-y^{{\prime\prime}{k}})/\vert {\boldsymbol y}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert$. Focusing again on the case of two point particles, equation (35) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{h}^{ij}(t,{\boldsymbol x}) \approx \frac{4G}{R c^4} \Big(m_1 \dot{x}_1^i \dot{x}_1^{\,j} + m_2 \dot{x}_2^i \dot{x}_2^{\,j} - \frac{G m_1 m_2 n_{12}^i n_{12}^{\,j}}{r_{12}}\Big) + \mathcal{O}\Big(\frac{1}{c^5}\Big). \nonumber \end{align} \tag{ 36 }$$

This expression agrees with equation (32). In conclusion, the discrepancy shown in figure 2 was simply due to neglecting the purely gravitational part of the source in the relaxed Einstein equations, which would be equivalent to assuming motion along straight lines for the binary components. In other words, the last term in equation (32) accounts for the fact that motion does not take place on rectilinear trajectories, but instead along curved-spacetime geodesics.

In this paper, we have highlighted several subtleties that occur when applying the quadrupole and quadrupole–octupole formulae to numerically-integrated binary and triple systems. We have shown that, as expected, applying these formulae to a binary in a reference frame whose origin is displaced from the CoM of the system leads to unphysical spurious results, if the displacement exceeds the wavelength λ of the emitted GWs. This simply happens because implicit in the derivation of the (generalized) quadrupole formula is the assumption that the multipoles are defined in the NCZ [6], i.e. a region of size  ∼λ centered on the origin in which the binary is supposed to be contained. The same problem manifests itself in hierarchical triple systems, when the quadrupole and quadrupole–octupole formulae are applied in a reference frame centered on the triplet's CoM. The resulting unphysical behavior, which to the best of our knowledge has gone unrecognized in the literature [19], is in contrast with what happens for a binary system (where it is safe to define the multipoles in the CoM reference frame), but can be understood by bearing in mind that a hierarchical triplet may be decomposed in an inner binary and an outer one. Indeed, as the separation of the outer binary grows, the CoM of the inner binary will eventually move out of its NCZ, thus violating the assumptions on which the derivation of the quadrupole and quadrupole–octupole formulae relies.

We have described two remedies to this problem. When using the leading order quadrupole formula, it suffices to express the waveforms in terms of appropriate conserved quantities (namely the total linear momentum) to eliminate the observed spurious behavior. When using the quadrupole–octupole formula, the simplest approach is instead to compute the multipoles in an inertial frame whose origin instantaneously coincides with the CoM of the inner binary. Remarkably, neither of these two fixes seems to work for a system of four bodies¹⁶, for which a more sophisticated approach should be developed.

Finally, we have shown that, if one were to compute the GW emission from a binary or triple system by integrating directly the equations for the linear metric perturbations over flat space, one would obtain GW amplitudes that are wrong by a factor  ∼2. We have found that this is related to the fact that the derivation of the quadrupole formula is quite subtle and actually requires one to use the non-linear Einstein equations.

The work of MB was supported in part by the ERC Project No. 267117 (DARK) hosted by Université Pierre et Marie Curie (UPMC)—Paris 6, PI J Silk. MB acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. For the numerical simulations, we have made use of the Horizon Cluster, hosted by the Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us. AS is supported by the Royal Society.

In this appendix we summarize the explicit leading order expressions for the quadrupole and octupole mass radiation, as well as for the current quadrupole radiation.

The second and third mass moments are defined from the time-time component of the matter stress-energy tensor as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle M^{ij}(t) &= \int \ud^3{\boldsymbol x} \ c^{-2} T^{00}(t,{\boldsymbol x}) x^{\langle i} x^{\,j \rangle}\nonumber \\ M^{ijk}(t) &= \int \ud^3{\boldsymbol x} \ c^{-2} T^{00}(t,{\boldsymbol x}) x^{\langle i} x^{\,j} x^{k \rangle}\,, \nonumber \end{align} \tag{ A.1 }$$

where $\langle \rangle$ represents the symmetric trace-free (STF) operator, i.e. (see [53] section 3.5.1)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x^{\langle i} x^{\,j \rangle} = x^i x^{\,j} - \frac{1}{3} \delta^{ij}r^2\,, ~ x^{\langle i} x^{\,j} x^{k \rangle} = x^i x^{\,j} x^k - \frac{1}{5} \left(\delta^{ij}r^2 x^k + \delta^{ik}r^2 x^{\,j} + \delta^{\,jk}r^2 x^i\right)\, . \nonumber \end{align} \tag{ A.2 }$$

The current quadrupole moment is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:curr_quad_mom} S^{ijk} = \int \ud^3{\boldsymbol x} \ [ x^i j^{\,kj} + x^{\,j} j^{\,ki}]\,, \nonumber \end{align} \tag{ A.3 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle j^{\,ij} = \frac{1}{c}\left(x^i T^{0j} - x^{\,j} T^{0i}\right) \nonumber \end{align} \tag{ A.4 }$$

is the angular momentum density tensor, which is connected to the angular momentum density vector by $j^{\,ij} = \varepsilon^{ij}_{~k\,}j^{\,k}$. Alternatively, one may define the STF quadrupole tensor

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle J^{ij} = \mathop{\rm STF}_{ij}\varepsilon_{abi} \int \ud^3{\boldsymbol x} \ x^{\langle\,j} x^{a\rangle} T^{0b}\,, \nonumber \end{align} \tag{ A.5 }$$

and replace $S^{abk}$ by $-2 \varepsilon_{~l}^{k~a} J^{bl}$ in the waveform, using the fact that $S^{abk} n_k$ and $-2 \varepsilon_{~l}^{k~a} J^{bl} n_k$ have the same transverse trace-free part (with respect to ${\boldsymbol n}$). The advantage of working with $J^{ij}$ is that it belongs to an irreducible representation of SO(3).

The expression of the GW waveform up to the next-to-leading order in the TT gauge is finally given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle h_{ij}^{\rm TT}(t,{\boldsymbol x}) &\approx \frac{2 G}{R c^4} \Lambda_{ijkl}({\boldsymbol n}) \Biggl[ \frac{\ud^2 M^{kl}(t_{\rm ret})}{\ud t^2} \nonumber \\ & \quad+ \frac{n_m}{3c} \frac{\ud^3 M^{ klm}(t_{\rm ret})}{\ud t^3} + \frac{2 n_m}{3c} \frac{\ud^2 S^{klm}(t_{\rm ret})}{{\rm d}t^2} \Biggr]+ \mathcal{O}\Big(\frac{1}{c^6}\Big)\,, \nonumber \end{align} \tag{ A.6 }$$

where the projector tensor $\Lambda_{ijkl}({\boldsymbol n})$ is defined in terms of the GW propagation direction ${\boldsymbol n}$ (see [54] section 36.10 and Box 35.1):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda_{ijkl}({\boldsymbol n}) = \mathcal{P}_{ik}\mathcal{P}_{jl} - \frac{1}{2}\mathcal{P}_{ij}\mathcal{P}_{kl}\,, \quad \mathcal{P}_{ij} = \delta_{ij} - n_i n_j\, . \nonumber \end{align} \tag{ A.7 }$$

In this appendix, we perform the explicit calculation of

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:H_generic} \bar{h}^{ij} (t,{\boldsymbol x}) = \frac{4G}{c^4}\int \frac{\ud^3{\boldsymbol x}^{\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert }\biggl[ T^{ij} + \frac{1}{4\pi G} \left(\partial^i\phi \, \partial^{\,j}\phi- \frac{1}{2}\delta^{ij}\partial_k\phi \, \partial^k\phi\right)\biggr] + \mathcal{O}\Big(\frac{1}{c^5}\Big)\, . \nonumber \end{align} \tag{ B.1 }$$

Since the first term of equation (B.1) is trivial to evaluate because it has a compact support, we focus on the second one, where one is not a priori allowed to approximate $\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert$ by R under the integral, since the integration extends up to spatial infinity. We perform the calculation by considering only the term $\partial^i\phi \, \partial^{\,j}\phi$. Indeed, $\delta^{ij}\partial_k\phi \, \partial^k\phi$, which is simply the trace of $\partial^i\phi\, \partial^{\,j}\phi$ multiplied by a Kronecker delta, disappears when taking the TT projection.

Inserting the expression for the Newtonian gravitational potential, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \label{eq:off_diag} \frac{1}{4\pi G}\int \frac{\ud^3{\boldsymbol x}^{\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert } \partial^i\phi\, \partial^{\,j}\phi &=\frac{G}{4\pi}\int \frac{\ud^3{\boldsymbol x}^{\prime} \ud^3{\boldsymbol y}^{\prime} \ud^3{\boldsymbol y}^{\prime\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert } \rho({\boldsymbol y}^{\prime})\rho({\boldsymbol y}^{\prime\prime}) \nonumber \\ &\quad\times \frac{\partial}{\partial x^{{\prime}{i}}} \left(\frac{1}{\vert {\boldsymbol x}^{\prime}-{\boldsymbol y}^{\prime}\vert }\right) \frac{\partial}{\partial x^{{\prime}{j}}}\left(\frac{1}{\vert {\boldsymbol x}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert }\right)\, . \nonumber \end{align} \tag{ B.2 }$$

After transforming the derivative $\partial/\partial x^{{\prime}i}$ that acts on $\vert {\boldsymbol x}^{\prime}-{\boldsymbol y}^{\prime}\vert ^{-1}$ into $-\partial/\partial y^{{\prime}i}$ by virtue of the translation invariance of ${\boldsymbol x}^{\prime}-{\boldsymbol y}^{\prime}$, and similarly for $\partial/\partial x^{{\prime}j}$, we may change the order of integration, so that $\partial^2/(\partial y^{{\prime}i} \partial y^{{\prime\prime}j})$ can be put outside the integral with respect to ${\boldsymbol x}^{\prime}$. With this trick, equation (B.2) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \frac{1}{4\pi G}\int \frac{\ud^3{\boldsymbol x}^{\prime}}{\vert {\boldsymbol x}-{\boldsymbol x}^{\prime}\vert } \partial^i\phi\, \partial^{\,j}\phi = G \int \ud^3{\boldsymbol y}^{\prime} \ud^3{\boldsymbol y}^{\prime\prime} \rho({\boldsymbol y}^{\prime})\rho({\boldsymbol y}^{\prime\prime}) \frac{\partial^2}{\partial y^{{\prime}{i}}\partial y^{{\prime\prime}{j}}} g({\boldsymbol x}^{\prime},{\boldsymbol y}^{\prime},{\boldsymbol y}^{\prime\prime}), \nonumber \end{align} \tag{ B.3 }$$

where g satisfies the Poisson equation $\Delta g({\boldsymbol x}, {\boldsymbol y}^{\prime}, {\boldsymbol y}^{\prime\prime}) = \vert {\boldsymbol x}-{\boldsymbol y}^{\prime}\vert ^{-1} \vert {\boldsymbol x}-{\boldsymbol y}^{\prime\prime}\vert ^{-1}$ in the sense of distributions. It is straightforward to check that the relevant solution is $g=\ln (\vert {\boldsymbol x}-{\boldsymbol y}^{\prime}\vert +\vert {\boldsymbol x}-{\boldsymbol y}^{\prime\prime}\vert +\vert {\boldsymbol y}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert)+{\rm constant}$ (see, e.g. p 355 in [69]), from which one infers the asymptotic behavior

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^2}{\partial y^{^{\prime}i}\partial y^{^{\prime\prime}j}}K = -\frac{\hat{n}^i \hat{n}^{\,j}-\delta^{ij}}{2 R \vert {\boldsymbol y}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert } + \mathcal{O}\Big(\frac{1}{R^2}\Big)\, . \nonumber \end{align} \tag{ B.4 }$$

At large distance R from the origin, equation (B.1) then reduces to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ud}{{\rm d}} \displaystyle \bar{h}^{ij} = \frac{4G}{c^4 R}\int \ud^3{\boldsymbol x}^{\prime} T^{ij} - \frac{2 G^2}{R c^4}\int \ud^3{\boldsymbol y}^{\prime} \ud^3{\boldsymbol y}^{\prime\prime}\rho({\boldsymbol y}^{\prime})\rho({\boldsymbol y}^{\prime\prime}) \frac{\hat{n}^i \hat{n}^{\,j}}{\vert {\boldsymbol y}^{\prime}-{\boldsymbol y}^{\prime\prime}\vert } + \mathcal{O}\Big(\frac{1}{R^2 c^5}\Big) . \nonumber \end{align} \tag{ B.5 }$$