Radiated momentum in the Post-Minkowskian worldline approach via reverse unitarity

We compute the four-momentum radiated during the scattering of two spinless bodies, at leading order in the Newton's contant $G$ and at all orders in the velocities, using the Effective Field Theory worldline approach. Following \cite{Mougiakakos:2021ckm}, we derive the conserved stress-energy tensor linearly coupled to gravity generated by localized sources, at leading and next-to-leading order in $G$, and from that the classical probability amplitude of graviton emission. The total emitted momentum is obtained by phase-space integration of the graviton momentum weighted by the modulo squared of the radiation amplitude. We recast this as a two-loop integral that we solve using techniques borrowed from particle physics, such as reverse unitarity, reduction to master integrals by integration-by-parts identities and canonical differential equations. The emitted momentum agrees with recent results obtained by other methods. Our approach provides an alternative way of directly computing radiated observables in the post-Minkowskian expansion without going through the classical limit of scattering amplitudes.


Introduction
Gravitational waves from black hole and neutron star binaries will provide an unprecedented source of information about astrophysics, cosmology and fundamental physics. This will be possible thanks to an increase in sensitivity of future detectors, such as LISA [2], Einstein Telescope [3] and Cosmic Explorer [4], which will require improving the modeling of the relativistic two-body dynamics over the current state-of-the-art [5].
Currently, waveform templates are modeled using semi-analytical approaches such as the effective-one-body formalism [6]. A crucial ingredient of this approach is the energy map between the two-body and the effective one-body systems, originally established using the post-Newtonian (PN) approximation [7,8]. This consists in expanding for small gravitational potential Gm/r, with m the typical mass of the two objects and r their relative distance, and small relative velocity between the two bodies v. In [9,10], it has been suggested that this mapping can be improved by using the post-Minkowskian (PM) approximation method, i.e. by expanding in the gravitational constant G without assuming a small velocity (see [11][12][13][14][15][16][17] for early works on the PM approximation). This has sparked fervent activity in the application of PM methods to the two-body relativistic gravitational dynamics.
While the PN expansion is natural for gravitationally bound systems-after all, the gravitational potential and velocity are related by the virial theorem-the PM expansion applies naturally to unbound systems, such as the scattering of two massive black holes. Information about their dynamics can be extracted by modeling the scattering black holes as a quantum mechanical system of two scattering particles. This can be studied on the basis of the long-time well-established quantum field theory description of gravity [18][19][20][21][22][23][24][25], by taking its appropriate classical limit. For a bound binary system, the classical limit consists in the orbital angular momentum L = mvr much larger than or, in natural units, L ≫ 1 [26]. For a scattering process L = pb, where p is the asymptotic center-ofmass momentum of the particles and b the impact parameter, so that the classical limit is obtained for q ∼ 1/b ≪ p, where q is the momentum exchanged in the scattering process. Such a small q expansion is analogous to the so-called soft expansion familiar from the method of regions [27].
Scattering bodies are accompanied by emission of gravitational Bremsstrahlung radiation [17,[59][60][61][62][63][64], which is the unbound analog of the gravitational waves emitted by inspiral binaries and is suppressed by three powers of G. Although radiation reaction effects were thoroughly investigated in the past in the Regge limit (i.e., when the center-of-mass energy is much larger than the momentum transfer) [65][66][67][68] or in association to the loss of angular momentum in the collision [69][70][71][72][73][74], the full leading-order emitted momentum has been obtained only very recently in [75,76] via the formalism of [77], which derives classical observables from quantum scattering (see also [78,79] for extensions of the formalism of [77] to spin and classical observables in Yang-Mills theories), and in [80] using the eikonal approach to classical gravitational scattering. These calculations require evaluating the classical limit of relevant two-loop Feynman integrals, that can be solved by combining different techniques borrowed from particle physics, as shown in [81], namely reduction to master integrals by Integration-by-Parts (IBP) identities [82][83][84] and differential equations [85][86][87][88] to solve the latter, using the near-static regime as initial conditions.
In this paper we focus on the computation of the emitted momentum at leading PM order using an EFT worldline approach inspired by Non-Relativistic-General-Relativity (NRGR) [26]. In the traditional NRGR approach (see [89][90][91][92][93] for reviews) one builds an EFT of classically radiating gravitons by exploiting the separation of the three relevant scales in the system: the size of the inspiraling bodies, the orbital radius and the wavelength of the emitted gravitational radiation. The potential-graviton propagators are expanded in the small velocity limit in a PN expansion and then integrated out. The bodies are treated as static background sources for the graviton dynamics, i.e., their recoil due to graviton interaction is neglected because suppressed by q/p ∼ 1/L ≪ 1, where q is the momentum of the exchanged graviton and p the typical momentum of the bodies. This approximation amounts to impose that the system is classical from the onset, dispensing one from the (sometimes tedius) counting. Quantum corrections, described by graviton loops, are suppressed by powers of 1/L relative to tree diagrams and can be ignored.
The worldline approach has been extended to the PM approximation [94][95][96][97][98][99]. In this scheme, one expands the body trajectories around rectilinear motion, each order in the expansion carrying an additional power of G. Combined with the powerful methods of IBP reduction and differential equations to solve the integrals involved in the calculations, this approach has allowed to quickly reach most of the state-of-the-art results achieved by scattering-amplitude techniques [100][101][102][103][104].
In [1] (see also [105]) we applied this approach to compute the conserved stress-energy tensor linearly coupled to gravity generated by the two bodies at leading and next-toleading order in G and, from that, the classical probability amplitude of graviton emission. The radiated four-momentum is given by a phase-space integral of the graviton momentum, weighted by the modulo squared of the radiation amplitude. At leading order, the radiation amplitude is just a static piece that does not contribute to the emitted energy. The nextto-leading amplitude, instead, contributes at leading order. It is given by an integral in the graviton momentum exchanged between the two bodies, but we were unable to perform this integral and write it in terms of known functions. (We note, in passing, that the Fourier transform of the amplitude, which is simply the waveform in time domain, can instead be performed, leading to the expression originally computed in [62,63] (see also [105]).) We could express it in compact form as a one-dimensional integral over a Feynman parameter involving Bessel functions. Using this, we recovered the leading-order radiated angular momentum [70] and, upon expansion of the integrand in v, the total four-momentum radiated into gravitational waves up to order v 8 , finding agreement with [75].
We review these developments in Sec. 2, spelling out details missing in [1]. In particular, we lay out the formalism of the EFT worldline approach for the PM calculation and introduce the Feynman rules. We then compute the stress-energy tensor and the classical amplitude of graviton emission at leading-and next-to-leading order. Appendix A contains the expression of the conserved stress-energy tensor in arbitrary spacetime dimensions. Contrarily to the more compact expression used in the main text, this stress-energy tensor is conserved both on-shell and off-shell. Appendix B details the treatment of the Feynman integrals involved in the amplitude computation.
In Sec. 3 we bypass the problem of not having a solution for the amplitude by rewriting the phase-space integral of the four-momentum as a (cut) two-loop integral. Using reverse unitarity [106][107][108][109], we treat the phase-space delta function as a cut propagator. Following [81], in Sec. 4 we solve this integral using techniques developed in QCD and high-energy physics. In particular, we organize our calculations in terms of four topologies that come out naturally from our Feynman rules for the gravitons and we solve each topology, one by one. 1 This shows that these techniques can be applied to the wordline approach to derive classical observables, without going through the classical limit of scattering amplitudes. Appendix C contains a thorough discussion on the boundary conditions necessary to solve the ordinary differential equations involved by the master integrals and details how to compute their solutions. Finally, we conclude in Sec. 5.
We use the mostly minus convention for the signature of the metric, the following notation for the integration symbol, . Unless otherwise specified, we will use natural units with = c = 1 and define the Planck mass as 2 From the action to the emission amplitude Following [1], we outline here the general set up of the PM worldline approach and derive the stress-energy tensor and, finally, the amplitude of classical gravitational emission.

Post-Minkowskian effective field theory
We consider a system of two spinless massive bodies with masses m 1 and m 2 , interacting via gravity, described by the Einstein-Hilbert action. As discussed in the introduction, we rely on the separation between the relevant scales in the system. We will assume that the impact parameter b involved in the collision between the two bodies is much larger than their typical size. We thus treat the two bodies as point particles. Finite size effects can be incorporated systematically as higher-derivative operators along the worldline (see e.g. [97,101]). The objects are considered as external (i.e. non-propagating) sources of the gravitational field. The resulting action, describing the dynamics of the system, is then where, for each body a = 1, 2, τ a is the proper time, U µ a (τ a ) = dx µ a /dτ a is the four-velocity and m Pl is defined in (1.2). Note that for the above action we have used a Polyakovlike parametrization, which has the advantage of simplifying the coupling of matter with gravity [97,110,111].
We want to compute the classical pseudo-stress-energy tensor T µν (x), defined as the linear terms sourcing the gravitational field in the effective action [20,26,112], i.e., where h µν (x) = m Pl (g µν − η µν ). This includes all contribution coming from both the external sources, i.e. the point-particles, and the gravitational self-interaction. We can compute T µν via a matching procedure. In particular, we expand the action (2.1) for small h µν and use this to compute the one-point expectation value h µν , considering all Feynman diagrams that involve one external graviton. We do the same using an effective action composed by the quadratic action of h µν plus an interaction term (2.2). We can then find T µν by matching the two results. Pictorially, we can depict this procedure as follows where the left-hand side stands for all the possible Feynman diagrams with one external graviton. On the right-hand side we have denoted the Fourier transform by a tilde, This procedure can be done order by order in the perturbative expansion in G.
Once T µν is known, we can use it to compute the classical probability amplitude of emitting one graviton with helicity λ and momentum k µ , defined by where ǫ λ µν (k) is the helicity-2 polarization tensor, with normalization ǫ * λ µν (k)ǫ µν λ ′ (k) = δ λ λ ′ . This is directly related to the asymptotic waveform (see e.g. [113]) by where r is a distance much larger than the interaction region, and with n the unitary vector pointing along the direction of propagation of the emitted graviton. Equation (2.5) can also be used to compute radiated observables such as the radiated linear momentum, as we will see in Sec. 3. Let us introduce, then, the Feynman rules relevant for the computation of the leadingorder and next-to-leading-order stress-energy tensor. As usual, for the gravitational sector of the action (2.1) we need to introduce a gauge-fixing term in order to define a propagator for h µν (x). We choose to work in the so-called de Donder gauge, which consists in adding the following gauge-fixing action, where h ≡ h µν η µν . Since we want only classical contributions, in (2.3) we exclude diagrams involving closed graviton loops, as they are purely quantum [26,92]. For the same reason, ghost terms are unnecessary. From the expansion of the gravitational action at quadratic order, we then define the usual graviton free propagator in de Donder gauge, for h µν , i.e., µν ρσ As usual, one must specify the contour of integration in the complex plane k 0 , and this can be done by choosing the suitable i0 + prescription in the denominator. Since we want to take into account only outgoing graviton, one should impose retarded boundary conditions, i.e. (k 0 + i0 + ) 2 − |k| 2 −1 . However, this is not relevant for the current computation, because, as we shall see, all the integrated graviton momenta are off-shell, so that we do not need to specify the i0 + prescription for these propagators. This will no longer be true at higher orders, where hereditary effects contribute to the computation of T µν [114,115]. Expanding the action at higher orders gives the self-interaction vertices. For the current computation we will just need the cubic vertex where we have expanded the action using the Mathematica packages xTensor and xPert [116,117] to compute V α 1 β 1 α 2 β 2 α 3 β 3
Finally, we need to write down the Feynman rules coming from the interaction of gravity with the external sources. As one can see from eq. (2.1), in principle we have just one linear interaction vertex. However, in order to completely isolate the powers of G, we expand perturbatively the worldline around straight motion [97,100], i.e.
Here u a is the (constant) asymptotic incoming velocity and b a is the body displacement orthogonal to it, b a · u a = 0. With this expansion, at leading order we obtain the following Feynman rules, where a bullet stands for the point particle and the cross attached to the wiggly line is there to remind us that there is no propagator attached to the straight worldline. At first order in G we have (2.14) For the computation of this paper we can stop at this order. In order to compute the first-order deviations from straight trajectories, δ (1) x µ a and δ (1) u µ a , one has to compute the effective action by integrating out the graviton from eq. (2.1), i.e. Varying S eff [x a ], one can then derive and solve the equations of motion for the two point particles. This procedure can be done perturbatively in the Newton constant G as explained in [97]. In general, one must carefully include all the contributions from both the potential and the radiation modes of h µν . However, the leading-order corrections δ (1) x µ a and δ (1) u µ a are determined by only potential gravitons. In particular, for the first order corrections we just need the following action in the de Donder gauge [97] S (1) One can then vary this action to find the equations of motion for the two worldlines, and then expand as in eqs. (2.11) and (2.12) to solve them perturbatively in G. For particle 1 one eventually gets and b ≡ b 1 − b 2 . An analogous expression holds for particle 2. The +i0 + in the above equations ensures to recover straight motion in the asymptotic past, i.e. δ (1) u µ 1 (−∞) = 0 and δ (1) x µ 1 (−∞) = 0. Eqs. (2.17) and (2.18) implies that the integrated graviton momentum is orthogonal to the timelike vector u 2 . Therefore, q can never be on-shell and hit the pole q 2 = 0, allowing us to ignore the i0 + prescription in the above three equations at this order.

Stress-energy tensor
As we explained in the previous section, we can compute the stress-energy tensor via a matching procedure, as depicted in eq. (2.3). At leading order in G, particles move along straight trajectories, generating a static term. Using the Feynman rule written in eq. (2.13), for body 1 we have Therefore, adding the symmetric contribution, we immediately find that The non-radiating nature of this contribution is manifest by the presence of the delta functions δ − (k · u a ).
At the next order, the stress-energy tensorT µν NLO is given by the sum of three contributions. The first is obtained when the worldline of the first body is deflected by the second one. Using the rule (2.14), we obtain The second contribution toT µν NLO is analogous to the first one, with the roles of the two bodies exchanged. The last contribution involves the cubic gravitational vertex and comes from evaluating the following diagram, Note that this separation in three contributions is only for convenience and depends on the chosen gauge.
Summing everything together, we can thus write the next-to-leading-order stressenergy tensor in Fourier space as 24) where t µν , t µν and t µν ⊢ come respectively from eq. (2.22), its symmetric under 1 ↔ 2 exchange and (2.23). They are explicitly given by and The integral in eq. (2.24) is over the momentum of the graviton exchanged by the two bodies. The two delta functions arise from the fact that we are taking the two bodies as non-propagating external sources. Note that similar integrals and delta functions appear when taking the classical limit of quantum observables in the scattering process between two massive particles. In this case, the integration variable q is the difference between the momentum within the wavefunction and that in its conjugate (the so called momentum mismatch [77]) while the delta functions arise from the on-shell constraints on the momenta of the scattering particles. Also note that we have simplified the expression of the stressenergy tensor (particularly that in (2.27)) by using momentum conservation as well as on-shell and harmonic gauge conditions, i.e. k 2 = 0 and k µh µν (k) = (1/2)k ν η αβh αβ (k). This simplification implies that the total stress-energy tensor in eq. (2.24) is transverse only for on-shell momenta, i.e. k µT µν ∝ k 2 = 0 [1]. In App. A we report the complete expression ofT µν NLO that is transverse also for off-shell momenta. Either expressions can be used to compute the emitted four-momentum, obviously leading to the same result.
Finally, we stress that we have left implicit all the i0 + prescriptions in the denominators appearing either from the graviton propagator to specify the contour of integration in the complex k 0 plane or from the corrections to the straight motion of the two bodies, eqs. (2.17) and (2.18). Concerning the gravitons, in order to take into account only outgoing radiation one should impose retarded boundary conditions, e.g. (q 0 + i0 + ) 2 − |q| 2 −1 .( 2 ) In our case, however, these prescriptions are irrelevant because, as displayed by eq. (2.24), the two delta functions ensure that the momenta q and q − k are orthogonal to one of the two fourvelocities. Thus, these two momenta can hit the pole only in the trivial case q = q − k = 0. The same holds true for the matter poles [k · u a ± i0 + ] −1 : since k is on-shell k · u a can never vanish.

Amplitude and waveform
We can now compute the classical amplitude A λ perturbatively in G using eq. (2.5). The leading-order contribution is obtained from the static contribution of the stress-energy tensor, eq. (2.21), The next-to-leading-order amplitude is of order G 3/2 . Analogously to what we did for the stress-energy tensor in eq. (2.24), we can separate it in three contributions where the labels refer to the contribution with the same name given in eqs. (2.25), (2.26) and (2.27). Introducing the following set of integrals we have explicitly that We stress that in eqs (2.33) and (2.34) one needs to exchange the label 1 ↔ 2 also inside the definition of the integrals I µ 1 ...µn . At this point, we are left with solving the integrals in eqs. (2.30) and (2.31). As discussed in [1], for the set I µ 1 ...µn (n) it is possible to find an analytic solution in terms on known function while for the set J µ 1 ...µn (n) the best we can do is to write them as one dimensional integrals over a Feynman parameter. See App. B for details of the calculations.
The next-to-leading-order amplitude takes a rather compact form if we consider the polarization tensor in the transverse-traceless (TT) gauge, i.e., and we choose the reference frame in which one of the two bodies, say body 2, is at rest, i.e.
where e v and e b are mutually orthogonal unitary vectors, With this choice A λ (k) = 0 and all but one term in the symmetric contribution in eq. (2.34) drop. Finally, parametrizing the graviton four-momentum as k µ = ωn µ , with n µ given in eq. (2.7), and defining z ≡ γ|b|ω one can write the next-to-leading-order amplitude in a compact form as where I, J = v, b and the coefficients A IJ are explicitly given by 3 where the coefficients c and d are given by In Ref. [1] we have compared this amplitude with previous results [63], even if only in particular limits. But eq. (2.38) can be integrated in the frequency using (2.6) and this results in the next-to-leading-order waveform in direct space, which fully agrees with that derived long ago by Kovacs and Thorne [63] (see also [105]). The emission amplitude can be also used to compute radiated observables, such as the linear and angular momentum loss by the system. The emitted angular momentum starts at order O(G 2 ) and involves the leading-order amplitude (2.28) and only the soft limit of the next-to-leading order amplitude. In the soft limit, the integrals in eqs. (2.39)-(2.41) can be performed and the waveform can be given analytically and used to compute the angular momentum [1], which agrees with previous results [70]. However, the leading-order radiated four-momentum involves the full next-to-leading order amplitude. Due to the involved structure of the integrals over the Feynman parameter y, the waveform can only be analytically computed expanding (2.38) for small velocities. In [1] we have used this expansion to derive the emitted four-momentum up to O(v 8 ). In the next section, we will take another route which dispenses with the need for an analytical expression of A NLO λ and which leads directly to the full emitted momentum.

Radiated 4-momentum as 2-loop integral
In term of the classical amplitude of graviton emission A λ (k), the radiated total momentum P µ rad is given by [1,95] is the Lorentz-invariant graviton on-shell phase-space measure. The differential probability of emission of one graviton with polarization λ is This quantity is not well-defined classically: if we interpret k and k 0 in these expressions as classical wave-vector and frequency, respectively, and we restore = 1, we must add a factor of −1 in front of the right-hand side, 4 which shows that the number of emitted gravitons is divergent in the classical limit → 0. However, inserting the four-momentum of the graviton k µ gives a finite quantity in the classical limit and integrating over all gravitons we obtain the total classical emitted momentum above. Pictorially, eq. (3.1) can be represented by where the on-shell amplitude on the right-hand side is non perturbative in G. Here we focus on the leading-order emitted momentum and therefore we expand the amplitude in powers of G as The first two diagrams on the right-hand side, of order G 1/2 , are static (they are proportional to δ − (k · u a )) and when multiplied by k µ they do not contribute to the emitted power.
Therefore, the leading order contribution to the radiated power comes from squaring the last three diagrams, of order G 3/2 , As explained in the previous section, we were unable to solve the integral in eq. (2.24) in the momentum of the graviton exchanged between the particles, q, and express the full amplitude (and waveform) in terms of known functions. In the following we adopt a different strategy to compute the right-hand side of eq. (3.1). Another pictorial interpretation of this equation is where we interpreted δ − + (k 2 ) as a cut propagator so that the modulo squared of the amplitude has been replaced by a vacuum-to-vacuum diagram with a cut. We also depict explicitly the flow of the momentum k µ as dictated by the positive energy theta function in δ − + (k 2 ). From (3.5), at leading order we expect four different cut topologies on the right-hand side, coming from the different ways of combining the three contributions in the modulo squared of the amplitude at leading order, denoted here by M, N, IY and H type, i.e., (3.7) Then, following a technique employed in high-energy physics computations known under the name of reverse unitarity [106][107][108][109], we can differentiate the cut propagator of eq. (3.6) like normal virtual propagators and apply the standard procedure of IBP reduction [82][83][84] to the resulting integral.
In practice, we rewrite the modulo squared of the amplitude as where we have used eq. (2.5) and Expanding the right-hand side of eq. (3.8) using (2.24), we obtain where the numerator N can be organized in terms of the contributions from the four topologies above and is explicitly defined as N (q 1 , q 2 , k) ≡ t µν (q 1 , k)+t µν (q 1 , k)+t µν ⊢ (q 1 , k) P µνρσ t ρσ (q 2 , k)+t ρσ (q 2 , k)+t ρσ ⊢ (q 2 , k) * . (3.11) -13 -Finally, replacing the modulo squared of the amplitude in eq. (3.1) using (3.10) and renaming we obtain (3.13) At this stage, we have rewritten the total four-momentum emitted as a cut 2-loop integral, followed by a Fourier transform from q to b-space. The advantage of this procedure is that we can now solve the 2-loop integral all at once, making use of the powerful computational tools routinely employed in high-energy physics-IBP reduction into master integrals [82][83][84] and differential equation methods [85][86][87][88] to solve the latter-without the need of deriving the Fourier-space gravitational waveform. This is analogous to the calculations recently performed in [75,76,80]. However, here we do not have to consider any intermediate quantum or super-classical contributions in our integrals: the amplitude A λ is a classical observable from the start. Before computing the contribution from each of the topologies in eq. (3.7) we need to discuss the master integrals that we will need to solve the associated two-loop integrals. This is what we turn to now.

Master integrals
As explained in [1], since the modulo squared of the amplitude is symmetric under the exchanged k · b → −k · b, the four-momentum cannot depend on the spatial direction b µ . Moreover, the energy measured in the frame of one body is the same as the one measured in the frame of the other one, hence the final result must be proportional to u µ 1 + u µ 2 . We can write this result as [75] P µ rad = We focus on the computation of E(γ), which we can be extracted by multiplying both eqs. (3.13) and (4.1) by u µ 1 + u µ 2 and comparing their right-hand sides. Therefore, we get with Notice that both E and I are dimensionless and only dependent on γ = u 1 · u 2 . Indeed, the 2-loop integral on the right-hand side of eq. (4.3) has dimension one. It can only depend on q 2 and γ = u 1 · u 2 because q · u 1 = q · u 2 = 0 by the delta functions in eq. (4.2) and no poles are expected. Since only q is dimensionful, it must scale as −q 2 , which is compensated by the prefactor. The integral in eq. (4.2) has dimension three and since the only dimensionful parameter is b, it must scale like |b| −3 . This removes the |b|-dependence on the right-hand side making E dimensionless. We now discuss how to simplify and solve the 2-loop integral I. Use the notation [76,80,81] and and rewrite it as We use dimensional regularization and extend the 4-dimensional integration to d spacetime dimensions, i.e.
The set of propagators in eqs. (4.4) and (4.5) and the four master integrals above are enough to solve our four topologies in eq. (3.7). At this point, we can use the differential equation methods [85][86][87][88] to solve these integrals. It is convenient to replace the dependence on γ of the master integrals by that on the kinematic variable x, defined by x ≡ γ − γ 2 − 1 [81]. The following relations derived from this definition will be useful later, Differentiating with respect to x, one realizes that the above integrals satisfy a system of differential equations of the form where f ≡ {f 1 , f 2 , f 3 , f 4 } and F (x, ε) is a matrix of rational coefficients. The properties of Feynman integrals ensure that the above system has only regular singularities, i.e., it is a Fuchsian system of differential equations. To solve this equation, it is convenient to find a basis g = {g 1 , g 2 , g 3 , g 4 } such that the differential equation is in canonical form [88,120,121], i.e., ∂ x g(x, ε) = εA(x) g(x, ε) . (4.12) A system of this form can be trivially solved in terms of polylogarithms as a Laurent series in ε. The transformation between the basis f and g can be obtained with the help of the package Fuchsia [122,123], implementing the Lee algorithm [124]. The canonical basis of master integrals reads g 1 = −q 2 G 2,0,0,1,0,1,1,0,1 , (4.13) g 2 = −q 2 G 2,0,0,1,0,0,1,1,1 , (4.14) g 3 = ε −q 2 γ 2 − 1G 1,0,1,1,0,0,1,1,1 , (4.15) g 4 = −q 2 5 γ − 1 8 G 2,0,0,1,1,1,1,1,1 + −q 2 1 − 2ε(2 + 3γ) 12(1 + 2ε) G 2,0,0,1,0,0,1,1,1 + 2ε (1 + 2ε)(1 + γ) −q 2 G 2,0,0,1,0,1,1,0,1 , (4.16) which satisfies the following canonically normalized differential equation, This can be equivalently written as with This differential equation can be solved perturbatively in ε [80], i.e., for each j = 1, . . . , 4, (4.20) The last ingredient are the boundary conditions of the differential equation (4.18). These can be found by solving the master integrals in the near static limit, i.e. for γ → 1 (or x → 1). We give an explicit derivation of these boundary conditions in App. C and we report here the results: where For the following computation we just need the solutions up to order ε. These are explicitly given by , g where γ E is the Euler-Mascheroni constant. Up to a different normalization of the loop integrals, 5 these agree with [80]. To conclude this section, we stress that before solving the master integrals g in the static limit x → 1 one must recast the cut propagators as delta functions. In practice, the resulting integrals are hard to solve. It is more convenient to relate these cut integrals to their non-cut versions using the so-called Cutkosky's rules [125]. The latter can then be solve using standard loop integral techniques. We explain all this in details in App C.1.

Computing the four topologies
We have now all the ingredients to compute the leading-order radiated momentum P µ rad . As mentioned in the previous section, we will focus on computing E(γ) defined eq. (4.2), splitting the computation in four contributions coming from the four topologies in eq. (3.7), i.e., where and  The numerators for each topology, N I , are defined below. The details of the calculation can be found in the ancillary files accompaning the arXiv submission of this article. In particular, using xTensor [117] the Mathematica notebook Contractions.nb computes the integrand of eq. (4.30) using the stress-energy tensor and prints the results in four different text files. These files are then imported in IBP-Basis1.nb, which performs the needed IBP reductions using LiteRed [118,119] and computes E I (γ) for each topology.

M topology
We start from the M topology, i.e. we solve eq. (4.29) with N M = P µνρσ t µν t * ρσ . (4.31) Performing the contractions and IBP reduction with LiteRed [118,119], one can eventually write this contribution in terms of a single master integral, where C M is a (not very illuminating) function of ǫ and γ. Here we are interested in the limit ε → 0. Since C M starts at ε 0 , we just need g 1 at order ε 0 , see eq. (4.23). After performing the Fourier transform in q in eq. (4.29) using we obtain The final result is unaltered by the exchange 1 ↔ 2, therefore the symmetric contribution gives exactly the same result.
Note that the M topology does not contain contributions from the graviton cubic vertex and the involved Fourier-space waveform (the amplitude) can be computed exactly. In this case we can also compute the above contribution in a more "direct" way from eq. (3.1), by taking the relevant part of the amplitude from Sec. 2.3. Specifically, E M can be computed from eq. (3.1) as with A λ (k) defined in (2.32). Working in 4 dimensions and solving explicitly the integral in q of eq. (2.32) as we did in App. B, we find the following expression in terms of modified Bessel functions of the second kind K n , where we have defined, for a = 1, 2, Using again eq. (3.9), we can rewrite eq. (4.35) as follows where we used that Working again in the frame (2.36), one can first solve the integrals in the azimuthal angle φ and the frequency ω, then finally in the polar angle θ, eventually recovering eq. (4.34). We stress again that such a direct procedure is unavailable when the amplitude involves A ⊢ λ (k) (see eq. (2.34)).

N topology
For the N topology the numerator in eq. (4.29) is Performing again the IBP reduction procedure and using the symmetry u 1 ↔ u 2 , we find that I N can be rewritten in terms of two master integrals, g 2 and g 3 , 41) where C N,2 and C N,3 are functions starting at order ε 0 . The coefficient in front of g 3 diverges for ε → 0 but this is compensated by g 3 that starts at order ε, see eqs. (4.23) and (4.26). Inserting the leading order solutions for g 2 and g 3 , we eventually find where we have used eq. (4.10) to replace − log (x) = 2 arcsinh γ − 1 2 . (4.43)

IY topology
For the IY topology the numerator in eq. (4.29) is After the IBP reduction we find that I IY is given in terms of g 1 , g 2 and g 3 , i.e.
where the dependence on ε and γ of the functions above is understood. Both C IY,1 and C IY,2 start at order ε −1 , leading to a seemingly divergent term for ε → 0, However, this is finite because both g 1 − g 2 and g 3 start at order ε. Inserting the solutions for g 1 , g 2 and g 3 given in eqs. (4.23)-(4.26), we eventually obtain where we used again eq. (4.43) and log (x + 1) 2 4x = log γ + 1 2 . (4.48)

H topology
Finally, we need to compute the contribution of the H topology, for which IBP reducing one last time, we find Once again, the cancellation of divergencies for ε → 0 is non-trivial. Before expanding g 1 , g 2 and g 4 we obtain a seemingly divergent term, which however is finite once we use the solutions for g 1 , g 2 and g 4 , eqs. (4.23)-(4.25) and (4.27). Using these, we obtain (4.52)

Full result
Summing up all the above contributions, i.e. eqs. (4.34), (4.42),(4.47) and (4.52), and taking into account also the symmetric ones, we eventually obtain This result agrees with the one recently derived via other methods [75,76,80], and complete the leading-order radiated sector derived with an EFT worldline approch [1].
From eq. (4.1) one can compute the center-of-mass radiated energy, E (CoM) rad ≡ P rad · u CoM , obtaining [75] where we have defined the total mass m ≡ m 1 + m 2 , the symmetric mass ratio ν ≡ m 1 m 2 /m 2 and h(ν, γ) ≡ 1 + 2ν(γ − 1). This result has been used to check with the literature in different regimes. For instance, one can compare against post-Newtonian computations up to 2PN [63,126,127] by expanding it for small velocities. From eq. (4.57), one can also obtain the radiated energy in elliptic orbits in the high ellipticity limit via analytic continuation [57,58] and compare the small velocity expansion with known 3PN results [7]. Finally, the same radiated energy also appears as a tail effect in the 4PM Hamiltonian computed in [45]. We refer to [75,76,104] for a more thorough discussion.

Conclusion
We have used the worldline approach to directly derive, for the first time using purely classical methods, the four-momentum radiated during the encounter of two spinless bodies at leading-order in the post-Minkowskian expansion, i.e. at O(G 3 ). The calculation can be roughly split into two parts. The first (Sec. 2) involves the derivation of the classical amplitude of one graviton emission, which was the subject of our previous work with Mougiakakos [1]. Here we have reformulated the steps leading to the stress-energy tensor and the onshell radiation amplitude, providing details omitted in that reference (see e.g. Apps. A and B). In the worldline approach the two bodies are treated as non-propagating external classical sources, as far as gravitons are concerned. In particular, quantum contributions are represented by graviton loops that we neglect. The derived stress-energy tensor and the radiation amplitude are classical objects from the onset. In contrast with scatteringamplitude based methods, there is no need to take an exponential expansion that leads to the treatment of unphysical super-classical terms [77] in the derivation. The second part of the calculation, reported in Secs. 3 and 4, concerns the computation of the radiated observable. This involves a phase-space integration of a single graviton fourmomentum weighed by the probability density of emission, i.e. the modulo squared of the on-shell amplitude. Using reverse unitarity, we treated the phase-space delta function as a cut propagator, allowing us to reformulate the integration as a two-loop Feynman integral that we attacked using the paraphernalia developed in particle physics: integration-byparts identities and differential equations. We split the integration into four topologies that naturally follow from the worldline approach, each of which can be reduced by integration by parts to a set of only four master integrals. To compute these master integrals we solved the system of canonical differential equations that they satisfy. Their boundary conditions, given in the near-static limit (γ → 1), are discussed in App. C. The explicit reduction to master integrals is instead reported in the ancillary files accompanying the arXiv submission of this work. In would be interesting to compare the intermediate steps of our derivation with those of other amplitude-or eikonal-based methods such as those of Ref. [75,77,80].
An obvious future direction is the extension of the present computation to higher orders. This will require to include the effect of the interplay between radiation and potential gravitons. We also expect the Feynman rules and the number of diagrams to increase rapidly. The classical double copy [94][95][96]128] could be a promising approach to reduce these complications. Moreover, it would be useful to apply techniques analogous to those discussed in [76,129] to simplify the computation of the boundary conditions of the master integrals. Other natural extensions, of course, include tidal/dissipative [130,131] and spin effects [132,133] in the radiated four-momentum.
Our approach is promising also to derive other radiated observables, such as the angular momentum at orders higher than G 2 . More generally, radiation-reaction effects will be required to incorporate the dissipative dynamics in semi-analytic models such as the EOB (see e.g. [134][135][136][137][138]) and will be crucial to develop accurate templates to make the most of future gravitational wave astronomy.

Acknowledgements
We would like to thank Brando Bellazzini, Carlo Heissenberg, Gregor Kälin, Stavros Mougiakakos, Julio Parra-Martinez, Rafael Porto and Leong Khim Wong for insightful discussions. This work was partially supported by the CNES.

A Stress-energy tensor
Here we report the expression of the stress-energy tensor conserved off-shell. For generality, we have computed it in d spacetime dimensions. Since the cubic vertex in de Donder gauge is the same in any dimension, the only modifications of the Feynman rules are the graviton propagator (2.9), which becomes µν ρσ and the first-order deviation of the equation of motion (2.17) and (2.18), which are now For convenience, we define The stress-energy tensor at next-to-leading order in G is given by the sum of the following terms 6 It can be verified that it is transverse for on-shell and off-shell gravitons as well, i.e., By using momentum conservation, on-shell and harmonic gauge conditions, we recover the expressions in eqs. (2.25)-(2.27).

B Integrals involved in the amplitude
This appendix is devoted to solve the set of integrals defined in eqs. (2.30) and (2.31). In particular, we show how I (n) can be solved exactly, while J (n) can be at best rewritten as an integral over a one-dimensional Feynman parameter. Let us start with the scalar integral I (0) . Since the final result will be in terms of Lorentz invariants, we can solve this integral in a particular frame and it is convenient to pick the one defined in (2.36). Solving the two delta functions, we reduce I (0) to a two-dimensional integral over the components q ⊥ that lie on the plane perpendicular to the direction of the scattering bodies. Hence (B.1) 6 We thank Gregor Kälin and Rafael Porto for several consistency checks of this result.
where in the second step we introduced a Schwinger parameter t. Solving the Gaussian integral in q ⊥ eventually gives the final result, i.e., where K n are modified Bessel functions of the second kind and we defined Once the scalar integral is solved, the vectorial I µ (1) can be computed decomposing it on a complete basis, i.e.
where the dependence on the combination u µ 1 − γu µ 2 comes from the fact that u 2 · I (1) = 0. Contracting both sides with b µ and u µ 1 , one eventually obtains The second set of integrals, defined in eq. (2.31), is more involved due to the presence of a second massless propagator. Starting again with the scalar integral J (0) , we can use Feynman parametrization to rewrite it in terms of only one massless propagator, where for the second line we have performed the shift q → q + yk and we have imposed the k 2 = 0, because the amplitude must be eventually evaluated on-shell. At this point we can follow a procedure analogous to the one we used for I (0) . Choosing again the frame (2.36), we use the two delta functions to reduce the computation to a two-dimensional integral over q ⊥ , that we can solve using Schwinger parametrization. This yields where we have defined Notice that s(y) changes when computing the symmetric contribution 1 ↔ 2. At this point, the integral over q ⊥ in eq. (B.7) is Gaussian and can be easily solved, where we introduced the shorthand notation We can solve J µ (1) and J µν (2) analogously to what we did for I µ (1) with the difference that, before decomposing on a complete basis as in eq. (B.4), we find it convenient to use again Feynman parametrization. For instance, for J µ (1) we have where we performed again the shift q → q + yk. The contribution proportional to k µ can be computed using the result of J (0) , while the one proportional to q µ must be decomposed on a complete basis. The very same procedure can be carried out for J µν (2) , yielding where we omitted (irrelenvant) terms proportional to k µ . To find the coefficients B i in eq. (B.12) we must solve the system (B.14) Using eq. (B.9), the solutions are yk · u 2 − (y − 1)γk · u 1 s(y) K 1 (w(y)) .

C Boundary conditions
In this last appendix we show how to compute the master integrals defined in eqs. (4.13)-(4.16) in the near-static limit to obtain the boundary conditions that we wrote in eqs. (4.21) and (4.22). We are going to follow closely the appendices of Refs. [76] and [80].

C.1 Cutting Rules
Uncut master integrals are easier to solve than cut ones. Hence, to connect a cut master integral to its uncut version we can use Cutkosky's cutting rules [125], 7 as explained in App. C of [76]. To use the cutting rules it is helpful to depict the master integrals g 1 , g 2 , g 3 and g 4 as diagrams. To do this, for convenience we introduce a "propagator" also for the massive external source, which can be seen as the soft-expanded version of the propagator of a massive scalar field [76,80,81]. Note that this is only a convenient pictorial tool useful to solve the Feynman integral. (The compact bodies are external sources and do not propagate.) For example, for g 1 we have the following topology where a thick line denotes the "massive propagator" and a thin line denotes the massless one. Figure 1 shows all the topologies required to solve the four master integrals of eqs. (4.13)-(4.16). We then introduce the notion of scalar integrals, which are basically Feynman diagrams in which one isolates all the factors of i coming from the non-cut propagators and the factors of −i coming from the vertices. To make a concrete example, let us consider again g 1 , represented by the diagram in Fig. 1(a), In the above equation, I 1 is the scalar integral. At this point we can find a relation between g 1 and I 1 using Cutkosky's cutting rules [125,139,140]. These can be derived through Veltman's largest time equation [141] and then listed as follows: • The sum of all cuts in a given channel is zero.
• All uncut propagators and vertices on the left-hand side of the cut are unaltered, while the ones on the right-hand side are replaced by the complex conjugate of their usual expression.
• Cut propagators are replaced by on-shell delta functions.
Given this set of rules, one can find a relation connecting cut integrals with their non-cut version. Using again the diagram depicted in Fig. 1(a) as an example, we have In terms of scalar integrals, this relation becomes (i I 1 ) * + g 1 + (i I 1 ) = 0 → g 1 = 2Im (I 1 ) .

(C.4)
Thus, to find the solution of the cut master integral g 1 in the near static limit it is enough to compute the imaginary part of the corresponding uncut scalar integral.
We can now do the same for the other three topologies depicted in Fig. 1. We have Fig. 1 Following the same procedure as for eqs. (C.3) and (C.4), we eventually find We can use the near-static limit of eqs. (C.4), (C.8), (C.9) and (C.10) to find the boundary conditions of the master integrals g i . Let us focus now on solving the scalar uncut integrals I i in the near static limit.

C.2 Integrals in the near-static limit
In what follows, unless stated otherwise we will always consider an implicit Feynman prescription +i0 + for all the propagators.
Integral g 1 Let us start by I 1 defined in eq. (C.2). Sending ℓ 1 → ℓ 1 + q − ℓ 2 we can separate the integration in ℓ 1 and ℓ 2 (C.11) We can solve the integral in ℓ 1 using eq. (10.25) of [84], (C.12) For simplicity, we now perform a Wick rotation to Euclidean space, i.e., for each vector v µ = (v 0 , v) = (iv 0 E , v E ), and we use the metric η µν E = diag(+, +, +, +) to contract the indices. The above equation becomes (C.13) Notice that (C.14) Using Schwinger parametrization we can rewrite the integral over ℓ E of eq. (C.13) as a Gaussian integral, i.e., Finally, for a = 1, 2, we can split the integrations in t a and s a , by simply performing the shift t a → √ s 1 + s 2 t a , obtaining The integration over s 1 and s 2 can be performed using standard integration over Feynman parameters. Making the change of variables s = s 1 + s 2 ,s = s 1 /s one gets  2ε) . (C.17) The integration over t 1 and t 2 is a bit more delicate. Changing again variables as follows t 2 = t t 1 , one can solve the integration over t 1 , Note that the integrand in t is divergent for t = γ − γ 2 − 1 = x and t = γ + γ 2 − 1, so one must treat it with care. In the near static limit x → 1 we obtain cos(πε) cot(πε) − i + O(1 − x) .

(C.19)
Putting all together we arrive to our final result for I 1 in the near static limit, i.e., (C.20) Using eq. (C.4) we can finally find the boundary condition for the master integral g 1 , i.e., where C BC has been defined in (4.22) and we used that sin(πε) .

(C.22)
Integral g 2 Let us now analyse the second scalar integral I 2 defined in eq. (C.5). First of all, we perform the shift ℓ 1 → ℓ 1 + q and then go again to Euclidean space for simplicity, (C.23) Using Schwinger parametrization and then solving the two Gaussian integrals, one eventually arrives to where we have defined [80] t 13 ≡ t 1 + t 3 , t 23 ≡ t 2 + t 3 , T ≡ t 1 t 2 + t 1 t 3 + t 2 t 3 . (C.25) Now we shift t 4 → √ T t 4 and t 5 → √ T t 5 , splitting the computation in two integrals As shown by eq. (C.9), finding the boundary condition of g 3 requires the solution of I 3 in the near static limit and also the computation of another cut of figure 1(c). First of all, following a procedure analogous to what shown for I 1 and I 2 , we find (C.36) To compute g 3 we have to subtract the last term on the right-hand side of eq. (C.9). Following the rules described in section C.1, we find where we have defined (C.38) (C.39) The integral I L is a simple one-loop computation that can be carried out straightforwardly using Schwinger or Feynman parametrization, obtaining . (C.40) Inserting everything in eq. (C.37) and solving the two delta functions, one eventually arrives to (C.41) Using eqs. (C.36) and (C.41) into (C.9), we get Integral g 4 Finally, we discuss the boundary condition for g 4 . Because of the factor γ − 1 in front of the first term of eq. (C.10), I 4 does not contribute to the boundary condition of g 4 , and we do not need to compute it. Using the results computed before for I 1 and I 2 , one can take the near-static limit of (C.10), obtaining (C. 43)