The Eikonal Approach to Gravitational Scattering and Radiation at $\mathcal O(G^3)$

Using $\mathcal N=8$ supergravity as a theoretical laboratory, we extract the 3PM gravitational eikonal for two colliding massive scalars from the classical limit of the corresponding elastic two-loop amplitude. We employ the eikonal phase to obtain the physical deflection angle and to show how its non-relativistic (NR) and ultra-relativistic (UR) regimes are smoothly connected. Such a smooth interpolation rests on keeping contributions to the loop integrals originating from the full soft region, rather than restricting it to its potential sub-region. This task is efficiently carried out by using the method of differential equations with complete near-static boundary conditions. In contrast to the potential-region result, the physical deflection angle includes radiation-reaction contributions that are essential for recovering the finite and universal UR limit implied by general analyticity and crossing arguments. We finally discuss the real emission of massless states, which accounts for the imaginary part of the 3PM eikonal and for the dissipation of energy-momentum. Adopting a direct approach based on unitarity and on the classical limit of the inelastic tree-level amplitude, we are able to treat $\mathcal N=8$ and General Relativity on the same footing, and to complete the conservative 3PM eikonal in Einstein's gravity by the addition of the radiation-reaction contribution. We also show how this approach can be used to compute waveforms, as well as the differential and integrated spectra, for the different radiated massless fields.

relations and colour-kinematics duality, render such calculations quite efficient as far as integrand construction is concerned, especially in the presence of supersymmetry, while systematic methods based on differential equations and integration-by-parts (IBP) identities automatise to a great extent the evaluation of the Feynman integrals.
The question then arises as to how classical gravity theories, in particular General Relativity (GR), emerge from such effective quantum theories in the classical limit, and in particular how classical observables can be recovered from perturbative scattering amplitudes, exporting as much as possible of the above technical tools. This question raises an immediate challenge: considering for simplicity four-dimensional 2 → 2 scattering in the centre-of mass-frame, where the total energy is given by √ s, the classical regime lies, for dimensional reasons, at trans-Planckian energies Gs ≫ , a manifestly non-perturbative kinematic region. The problem becomes tractable in the post-Minkowskian (PM) approximation, where one expands in powers of the dimensionless ratio G √ s/b ≪ 1, with b the impact parameter: the classical PM regime is captured by the Lorentz-invariant on-shell amplitudes when enforcing the hierarchy of length scales where √ s and G √ s can be regarded as the effective Compton wavelength and Schwarzschild radius, respectively. Notice that the first inequality in (1.1), which is equivalent to G √ s ≫ ℓ P in terms of the Planck length ℓ P = √ G , is incompatible with the perturbative expansion of quantum field theory.
Luckily enough, gravitational scattering at high energy and large impact parameter was extensively studied in the late eighties [1,2,3,4,5] and the fundamental progress achieved then indeed sheds light on the problem outlined above. The first process studied in the regime (1.1) was the elastic gravitational scattering of two massless states which, to lowest order in the perturbative expansion, is given by the one-graviton-exchange diagram. The graviton couples to the energy itself, and the tree-level scattering amplitude diverges strongly at high energies, violating partial-wave unitarity bounds. But it was soon understood that unitarity is restored by the observation that the lowest-order tree diagram, when looked at in impact-parameter space, is just the first term of an exponential series which corresponds to summing over the exchange of arbitrarily many gravitons. This non-perturbative resummation yields a phase factor exp(2iδ 0 ), called the leading eikonal, which can be used to compute classical observables such as the deflection angle or the Shapiro time delay. It turns out that, at this lowest order, the eikonal phase is nothing but the phase-shift obtained for the scattering of one massless particle in the Aichelburg-Sexl shock-wave metric [6] produced by the other particle [2]. Let us mention as an aside that the original motivations for studying trans-Planckian-energy collisions of light particles (or strings) were essentially theoretical, in the spirit of the old gedanken experiments of quantum mechanics. They had to do, in particular, with solving Hawking's celebrated information paradox [7] by trying to construct a unitary S-matrix even in the regime G √ s < b in which one expects the collision to lead to black-hole formation.
In spite of some interesting progress [8,9,10], this has remained, so far, an unfinished program.
For present-day purposes, it is important to stress that the gravitational eikonal exponentiation is not restricted to the case of massless scattering, as already pointed out in [11]. On the contrary, it has become clear over the last few years that it consistently provides the necessary non-perturbative resummation to extract the proper classical limit of the amplitude and to calculate relevant classical observables in many different setups, including the scattering of objects with masses m 1 and m 2 , in the whole range m 1 + m 2 ≤ √ s < ∞, i.e. for generic relative velocities between the two bodies [12,13,14,15,16].
For these reasons, since the first direct detection of gravitational waves, the emphasis in the study of the eikonal exponentiation has naturally shifted towards its application to the classical relativistic two body problem, for which various GR techniques, like the PM expansion mentioned above, the post-Newtonian (PN) expansion, and the effective-onebody (EOB) approach [17] had been developed in the past to describe the inspiral phase of two merging black holes. In particular, as emphasized by Damour [18], quantum scattering amplitudes can provide useful inputs to the EOB description of black-hole binaries by fixing some of the parameters appearing in its effective potential, and there has been considerable interest in the use of amplitudes-based techniques for extracting the conservative part of the EOB potential [19,20,21,22,23,24,25,26,27,28,29,30].
The techniques developed for the massless eikonal must be adapted to make room for the massive scattering problem at arbitrary relative velocities, while keeping suitable constraints to enforce the classical regime and ensure the validity of the post-Minkowskian approximation. To this effect, we shall incorporate the non-trivial mass scales m 1,2 in the hierarchy (1.1) according to √ s m 1,2 ≪ Gm 1,2 G √ s ≪ b , (1.2) or more precisely, since dimensional regularization is needed to deal with the standard infrared divergences that arise along the way, using its (4 − 2ǫ)-dimensional counterpart (2.20). In this setup, the Schwarzschild radii Gm 1,2 are taken to be much larger than the two types of quantum scales corresponding to the Planck length ℓ P ∼ √ G and to the Compton wavelengths m 1,2 , so that the leading eikonal 2δ 0 effectively scales as Gm 1 m 2 / ≫ 1 and captures the 1PM two-body dynamics. Higher-order PM corrections are instead parametrised by the quantity Gm 1,2 /b ≪ 1, which is indeed the typical size of Einstein's deflection angle, for arbitrary values of s m 1 m 2 . In this way, 2δ 1 is suppressed by an extra factor of Gm 1,2 /b compared to the leading eikonal 2δ 0 and captures the 2PM two-body dynamics. In the massive case, this 2PM eikonal 2δ 1 is non-trivial and its Ddimensional GR expression was discussed in [12,28]. In the massless case, the analogous 2PM contribution vanishes for kinematic reasons, and the same happens in the case of maximally supersymmetry gravity in D = 4 even when the external states are given a susy-preserving mass (for instance by assuming that they have a non-zero Kaluza-Klein momentum) [31]. Thus in the latter cases, the first non-trivial correction to the eikonal, 2δ 2 , comes at two loop (or 3PM) order.
In the massless, non-supersymmetric case, the 3PM eikonal was first derived more than thirty years ago [32]. Unlike its lower-PM counterparts, 2δ 2 has both a real and an imaginary part: the latter was computed in [32] from the three-particle cut involving the two massless energetic particles and a graviton, and Re 2δ 2 was determined from Im 2δ 2 through arguments based on analyticity and crossing. The analysis of [33] suggested that the same value of Re 2δ 2 should hold also for supersymmetric theories and this was confirmed in [34] for N = 8 supergravity. This universality was finally extended to massless theories with different amounts of supersymmetry including GR [35]. Universality of the real part depends crucially on a universal log(s)-enhanced term in the Im 2δ 2 , while the rest of the imaginary part is not universal and depends on the spectrum of massless states that contribute to the three-particle cut.
In the massive case, the first GR results at 3PM order [36,37] used the matching with an effective theory of interacting scalars to extract a conservative Hamiltonian relevant for the classical dynamics [22], which could be directly exported from the scattering setup to the case of bound trajectories. These results were later confirmed in [38,39], while the 3PM analysis in the maximally supersymmetric theory is more recent [13].
Before commenting further on these results, an important practical observation is in order. Of course, the program of extracting classical physics from scattering amplitude would hardly simplify matters compared to the PM expansion of classical gravity if one had to first evaluate the exact quantum amplitude and then take its classical limit/resummation. The goal is therefore to drop purely quantum contributions, such as short-range interactions localized near b = 0, as early as possible in the calculation of the amplitude. This can be achieved systematically by performing the near forward expansion of the amplitude, i.e. the expansion for small q ∼ b , via the method of regions [40,41]. This tool simplifies the integration by reducing the expansion of a given Feynman integral to the expansion of its integrand in suitable scaling limits, called regions, which leads, in turn, to the evaluation of much simpler integrals. In this language, neglecting analytic-in-q 2 terms in the amplitude simply amounts to dropping the so-called hard region, defined by the scaling ℓ ∼ O(m 1,2 ) for the momenta ℓ of the gravitons exchanged, and focus on the soft region, ℓ ∼ O(q), which contributes to the long-range effects.
All the investigations of massive 3PM scattering mentioned above focused instead on the contributions arising from the potential region, a sub-region of the soft region defined for small velocities v by the non-relativistic scaling ℓ 0 ∼ O(vq), ℓ ∼ O(q), hence breaking manifest Lorentz invariance. The rationale behind this further expansion was that the potential region, tailored to capturing conservative effects, would directly provide the conservative EOB data. Moreover, the restriction to the potential region allows one to consistently drop certain bubble-like contributions already at the integrand level, thus further simplifying the task at hand. However, the resulting conservative 3PM deflection angle exhibited a log(s)-enhancement in the high-energy limit √ s ≫ m 1,2 . This feature, which is tantamount to a log(s)-enhancement in Re 2δ 2 , was also claimed to be universal [13], in tension with the already mentioned and supposedly universal ultrarelativistic massless result [32], which did not show any sign of such an enhancement. This apparent discontinuity between the high-energy limit of massive scattering and massless scattering prompted a lot of discussion and some attempt, such as the one in [42], was made to solve it, but was eventually contradicted by subsequent checks at 6PN level [43,44]. This contradiction became even sharper with Ref. [15], where we first extended the analyticity-plus-crossing analysis of [32] to the ultrarelativistic massive case arriving at exactly the same conclusions. However in Ref. [15] we also presented for the first time the solution of the puzzle. Following [13], but performing the calculation of the various loop diagrams taking into account the full soft region -rather its potential part-we computed the 2 → 2 scattering amplitude and the associated eikonal Re 2δ 2 in massive N = 8 supergravity at 3PM level and we found that the resulting deflection angle is free of log(s)-divergences, thus restoring agreement with [32] and the universality in the ultra-relativistic limit. The physical interpretation of the new contributions emerging in this calculation was immediately suggested by Damour: the potential region provides only the conservative part of the deflection angle, while the soft region provides also the effects of dissipation due to the emitted gravitational waves, i.e. radiation reaction effects. While they are crucial to restore universality in the ultra-relativistic limit, these radiation-reaction terms can be more easily identified in the non-relativistic regime as those associated with odd powers of the velocity, or equivalently with half-integer PN order. In conclusion, while the unphysical "conservative deflection angle" has a spurious high-energy divergence, the physical (in principle measurable) deflection angle has a smooth limit.
By adding to the conservative result obtained in [36,37] a linear-response contribution [45] due to an O(G 2 ) radiative loss of angular momentum 1 , Damour [47] managed to extend to GR the result of the physical deflection angle presented in [15] for N = 8, obtaining the cancellation of the log(s)-divergence also in Einstein's gravity. His result has been confirmed by a completely different, amplitude-based approach in [49]. This alternative method consists in computing the infrared-divergent part of the three-particle unitarity cut due to the emission of a soft graviton. As in the massless case, this determines the infrared-divergent part of Im 2δ 2 (which can be equivalently derived from the exponentiation of infrared divergences in momentum space [50]) and also a specific logarithmic term in the finite part of Im 2δ 2 . This in its turn uniquely fixes the radiationreaction contributions to Re 2δ 2 via analyticity, and therefore the needed correction terms to retrieve the physical deflection angle from the conservative one.
In this paper we provide a fully self-contained discussion of the 3PM eikonal in massive N = 8 supergravity. In particular, we give the technical tools necessary to verify the N = 8 result for Re 2δ 2 presented in [15] and also derive the corresponding Im 2δ 2 . We extend the analysis of the imaginary part also to the non-supersymmetric case and, by building on previous results [36,37,47,49], present the full 3PM eikonal for the scattering of massive scalars in Einstein's gravity. As already mentioned, this requires evaluating the classical limit of the relevant two-loop Feynman integrals in the full soft region, and we achieve this goal by extending the method of ordinary differential equations (ODE) developed in [13], combining ODE and IBP-reduction techniques with the complete near-static boundary conditions 2 . This allows us to find from first principles the full 3PM eikonal in N = 8 and to elucidate the analytic structure of the classical result by writing it explicitly in a real-analytic and crossing symmetric way. We also extend the analysis of the 3-particle unitarity cut of [15,49] to generic centre-of-mass energies and beyond the Weinberg limit for the intermediate massless states, to include all classical contributions relevant to the scattering of massive states. This allows us to obtain the complete expression of Im 2δ 2 in GR, thus generalising the result of [32]. As a simple extension of this analysis, we are also able to calculate the total energy emitted at 3PM in the scattering process, both in N = 8 and in GR, finding perfect agreement with the results of [52]. We finally describe how our approach can also be applied to obtain waveforms, as well as full differential emission spectra, for the various massless fields as functions of the frequency/retarded time and the angles. For GR these results date back to [53], while the complete spectrum has been recently computed in [54,55] by using the world-line formalism and diagrammatic techniques (see also [56,57] for previous results in the massless case and [49] for the zero-frequency limit in the massive case).
The detailed outline of the paper is as follows. In Section 2 we describe the kinematics of the process under study and then give a quick reminder of the eikonal approach to connect quantum scattering amplitudes to classical gravitational observables. In Section 3 we give details on the tools used to perform loop calculations, both for the elastic 2 → 2 amplitude and for the 3-particle unitarity cut thereof. These are based on IBP techniques and ODE tailored to the evaluation of Feynman integrals in the soft region. In Section 4 we present the actual results for the massive N = 8 case up to two loops in momentum space. Integration covers the full soft region and is manifestly relativistic invariant. We then compute the corresponding eikonal phase and physical deflection angle. In Section 5 the same results are recast in an explicitly analytic and crossing symmetric form. This allows us to connect some contributions to the real part of the amplitude to corresponding imaginary parts endowed with a simple physical interpretation (e.g. elastic vs. inelastic contributions). One such connection is the shortcut used in [49] in order to obtain the radiation-reaction term of the deflection angle. In Section 6 we study the 3-particle unitarity cut of the elastic 2 → 2 amplitude, which we obtain directly from the appropriate "square" of the inelastic 2 → 3 amplitude at tree level. We show how the corresponding integrations can be carried out in momentum space, borrowing tools from the reverse unitarity method [58,59,60,61,52]. These techniques can be used to calculate Im 2δ 2 , which allows us to reconstruct the full GR eikonal up to 3PM, and the total radiated energy-momentum, for which we find agreement with Ref. [52]. We also indicate (and show in an example) how our method can be adapted to obtain the wave-forms in impact parameter space and to calculate the full differential spectrum of the emitted massless field waves from a completely amplitude-based approach. Section 7 summarizes our results and gives a brief outlook. The two appendices contain additional technical material. Appendix A shows how to obtain the complete near-static boundary conditions that are needed for the differential equations in Section 3 in order to encompass the contributions of the full soft region. In Appendix B we give the PN expansion of our results illustrating how the radiation-reaction terms correspond to odd powers of the velocity in the PN expansion of the real phase while containing logarithms of v in the imaginary parts.

Eikonal Approach for → Scattering
We are interested in the classical scattering of two spinless massive objects with masses m 1 and m 2 . Our approach is to start from the elastic 2 → 2 amplitude depicted in Fig. 1, which is evaluated in the usual perturbative expansion in powers of the Newton constant G, and then apply the eikonal resummation to extract its proper classical asymptotics. In this section, we first spell out our conventions concerning the process under scrutiny and then provide a short summary of the eikonal strategy for the calculation of the classical deflection angle.  Figure 1: Schematic representation of a 2 → 2 scattering process involving two objects with masses m 1 (blue) and m 2 (green).

The process under study
With reference to Fig. 1, we consider formally outgoing external momenta p µ 1 , p µ 2 , p µ 3 and p µ 4 subject to the conservation and mass-shell conditions The relative speed v of the asymptotic states is sized by the relativistic factor which tends to 1 in the static limit, v → 0, and to +∞ in the ultrarelativistic regime, v → 1. We work in the centre-of-mass frame and, following [13], we introduce the variablesp µ 1 ,p µ 2 and the perturbative transferred momentum q µ according to The usefulness of the new variables lies in the orthogonality relations which imply that q µ is constrained to a (D − 2)-plane, although we still regard it as a D-vector. We also introduce the parameter which, for small q, is related to σ of Eq. (2.2) by Here and in the following, we also employ q as a shorthand notation for q 2 .
In order to rationalize certain square roots that naturally emerge, it is also convenient to further define and similarly for log 1 x . Comparing with Eq. (2.2), we see in particular that log 1 z plays the role of a relative rapidity in this context, since Finally, we will also employ the shorthand definition especially in the discussion of differential equations and their boundary conditions. In terms of the variables introduced above, the Mandelstam invariants take the form Let us also record here some useful relations that link p, the momentum of either particle in the centre-of-mass frame, to the centre-of-mass energy E = √ s:

A reminder of the eikonal approach
As usual, we decompose the S-matrix according to so that the elastic 2 → 2 amplitude A(s, q 2 ) is defined by the matrix element Its perturbative expansion in powers of the coupling G takes the form with A j (s, q 2 ) ∼ O(G j+1 ) the j-loop contribution. For brevity the dependence of A on the masses m 1 , m 2 and on the dimensional regulator ǫ = 4−D 2 is left implicit. The classical weak-coupling regime obtains considering the following hierarchy of length scales, where the two quantities in the middle are constructed in analogy to the Schwarzschild radius. The first inequality ensures that quantum effects be negligible, i.e. suppressed by the ratio between the Compton wavelength and the size of the objects under consideration, while the last inequality allows us to focus on a regime of weak gravitational interactions, i.e. of near-forward scattering. In this near-forward limit, namely for small momentum transfer q 2 → 0, each of the above contributions to the amplitude in Eq. (2.19) can be schematically expanded as up to terms that are either proportional to non-negative integer powers of q 2 or subleading for small q 2 . The superscripts refer to the power counting, which marks the appear- terms. The most natural interpretation of the singular, super-classical terms is that such contributions arise from the formal small-G expansion of an exponential of the type e 2iδ that oscillates infinitely rapidly in the classical limit, where its phase Re 2δ becomes large. To see that this is indeed the case, however, one needs first to rewrite the amplitude in terms of the (D − 2)-dimensional impact-parameter b µ according to In impact-parameter space, the following exponentiation becomes manifest, where δ(s, b) is the classical eikonal, while ∆(s, b) is a quantum remainder. Employing term by term the formula

25)
where for instance and analogous relations link the f 's to the a's for higher-order terms. Expanding δ and ∆ appearing in Eq. (2.23) in powers of G, the eikonal resummation amounts to identifying and checking that The eikonal data is therefore determined by the classical and quantum terms of the amplitude, while the superclassical terms should be simply regarded as redundant objects.
Once the eikonal phase is properly retrieved, the inverse Fourier transform leads one to identify the total classical momentum transfer as Distinguishing between the real and imaginary part is crucial because, starting from two-loop order, the eikonal acquires an imaginary part connected to the suppression of the probability of the purely elastic process compared to the inelastic process including Bremsstrahlung. The deflection angle χ can be obtained by noting that Note that the impact parameter b µ used in all previous equations is by construction oriented along the direction of the momentum transfer Q µ . It is sometime more informative to recast the final result as a function of the impact parameter b µ J that is orthogonal to the asymptotic momentum in the centre-of-mass frame, so that the angular momentum J is given by pb J = J. To this effect, it is sufficient to note that so that tan This is the basic relation that we are going to employ in the following in order to retrieve the scattering angle χ from the eikonal.

Unitarity in b-space
In this paper we shall make use of (perturbative) unitarity constraints in impact parameter space up to the two-loop level. In momentum space the standard unitarity equation We will consider (2.38) between our initial and final two-particle states, using the loop expansion (2.19), but, of course, enforcing the full unitarity constraint involves elastic as well as inelastic amplitudes. Below we shall briefly review the implications of these constraints on the elastic impact-parameter amplitude as parametrized in (2.23). At tree level (2.23) simply defines the leading eikonal phase 2δ 0 which is real. If one could neglect all inelastic intermediate states, (2.23) with δ = δ 0 and ∆ = 0 would be the exact solution of the unitarity constraint, modulo a subtle point to be mentioned shortly.
Up to two loops, the possible intermediate states entering the right-hand side of the unitarity equation (2.38) consist of either two or three particles. Given that there is conservation for each heavy-particle species, the possible two-particle intermediate states are either the same as the initial/final ones or their (fermionic or KK) partners. The possible 3-particle intermediate states, occurring first at two-loop order, contain in addition a massless particle.
At one loop, the unitarity constraint (2.38) reads: where we indicated by ⊗ the standard on shell convolution with the appropriate phase space factors (see e.g. Section 6) and by A (in) 0 a generic 2 → 2 inelastic tree-level amplitude. On the other hand, going to impact parameter space and using (2.23) gives: where we indicated by F T [ · ] the Fourier transform according to (2.22) and the dots correspond to a possible lack of exact factorization when one converts the convolution in momentum space to a product in impact-parameter space. It turns out that in our frame (in which q is orthogonal to both p 2 − p 3 and As a consequence, 2 Im δ 1 = 0 while 2 Im ∆ 1 is just related to the inelastic two-body channels mentioned above 4 . In general, ReÃ 1 has a non-vanishing classical real part of O( −1 ), denoted by 2 Re δ 1 in (2.29), as well as a quantum contribution (O( 0 )) denoted there by 2 Re ∆ 1 . As discussed later in the paper, ∆ 1 should be kept up to O(ǫ) since it intervenes in the determination of δ 2 . It will be ignored here since it does not play an important role in the present discussion.
Consider now (2.38) at two-loop order. In momentum space it takes the form: where we have used the fact that the tree-level amplitudes A 0 and A (2→3) 0 are real. On the other hand, assuming eikonal exponentiation as in (2.23): that can be separated into an equation for its real and one for its imaginary part. For the latter we get: Using again that the two-particle convolution in momentum space factorizes in b space, we find: 44) and inserting this in the Fourier transform of (2.39) gives This is the sought-for connection between Im δ 2 and the three-particle cut of the elastic amplitude to be used later in this paper. We will show there how the r.h.s. of (2.45) can be factorized through a suitable definition of the Fourier transform of the 2 → 3 amplitude. We also point out that 2 → 3 processes in which the final state does not consist of the initial two particles and a massless quantum do not contribute to the classical eikonal δ 2 , but only to its quantum counterpart ∆ 2 in (2.29) which is not of interest at the 3PM level. The above arguments can be generalized to higher loops but we shall not need them in the context of this paper.

IBP+ODE Toolkit
In this section we apply the method of differential equations [13] to the calculation of the integrals that are needed to evaluate the 2 → 2 amplitude in the limit of small momentum transfer up to two-loop order. In general, the flowchart of the method can be summarized as follows: 1. Express the 2 → 2 amplitude in terms of elementary (possibly scalar) integrals I T , each corresponding to a certain topology T (which, by following the notation of [13], we dub II, III, IX, H and for their crossed counterparts II, III, IX, H); 2. Consider the limit of small momentum transfer q → 0 and expand each integrand in the soft region, defined by the scaling ℓ i ∼ O(q) of the loop momenta ℓ i ; 3. Perform Integration By Parts (IBP) reduction to identify a basis of master integrals and find the differential equations that link them to one another; 4. Solve these equations order by order in ǫ using as boundary conditions the values of the integrals for small velocities.
In step number 2 one neglects contributions arising from the hard region, ℓ i ∼ O(q 0 ), which amounts to disregarding terms proportional to non-negative integer powers of q 2 in the amplitude. This is fully justified, for the present purposes, because it amounts to neglecting contributions to the eikonal whose support is localized at b = 0 in impactparameter space.
Step number 3 is automated using LiteRed [62,63] and FIRE6 [64]. For N = 8 supergravity, steps 1-3 were already performed in [13] (except for the odd-q master integrals associated to the H family), where however step 4 was simplified by further restricting the evaluation of the integrals to the potential region, in the static limit.
In this section,we perform step 1-to-4 for the full soft region, evaluating the static limit of the master integrals without any extra ingredient (in particular, without restricting to the potential region). We first discuss ordinary (planar and non-planar) topologies arising at one and two loops, and then turn to the analysis of the integrals subject to the three-particle cut.

Box and crossed box
The relevant integrals at one loop are and i k are integers. The association between the labels i 1 , i 2 , i 3 , i 4 and the internal lines is depicted in Fig. 2. The master integrals can be taken as follows The scalar box integral 5 can be then expanded in the soft region and IBP-reduced to obtain where (3.6a) is superclassical, (3.6b) is classical, (3.6c) is the first quantum correction and we neglect higher-order quantum contributions. The differential equations for the master integrals are conveniently expressed in terms of the parameter x defined in Eq. (2.10) and read Their solutions are uniquely determined by fixing the integration constants c 1 , c 2 , c 3 using the boundary conditions at x = 1. For the first two master integrals, one obtains by standard methods and In this simple case, f 3 can be also evaluated by standard methods for any x, obtaining consistently with the general solution (3.8). Substituting back into the decomposition given by Eqs. (3.6a), (3.6b) and (3.6c), one then finds, up to subleading orders in q, The crossed box is given by replacing p 1 ↔ p 4 in the box. It can be obtained by first expressing I II in terms of z and then implementing crossing symmetry according to The result is up to subleading orders in q.

Planar double box
In this case the relevant soft integrals are . (3.15) Here the first four entries are associated to massive lines in the soft limit, while i 5 , i 6 , i 7 correspond to massless lines as depicted in Fig. 3. The last two entries are needed to complete the basis. The Feynman prescription has been left implicit but it can be restored by adding −i0 on the right-hand sides of (3.16), (3.17) and (3.18). We consider the master integrals which depend on even powers of q, and which depend on odd powers of q. The scalar double box integral can be decomposed as follows in the soft limit, and (3.24g) The differential equations for the master integrals decompose into two disjoint sectors, depending on even and odd powers of q, where the matrices are given by In order to solve the above differential equations, we need to first determine the static boundary conditions by studying x-dependence of the master integrals in the limit x → 1. The details of this calculation are presented in Appendix A. For the even sector, we have The solutions of the differential equations can be then determined in a perturbative fashion order by order in ǫ. Letting for j = 1, 2, . . . , 10, allows one to recast (3.25) in terms of the recursion relations d f which can be then solved for k = 0, 1, 2, . . . taking into account the ǫ-expansion of the boundary conditions at each step. We find, for the even sector, (3.34) In these expressions x is understood to have a small negative imaginary part, so that log(−x) = log x + iπ. For the odd sector, we have instead These explicit expressions for the master integrals in their turn determine the decomposition coefficients defined in (3.22), for which we find up to O(ǫ 0 ). In these equations, z is understood to have a small negative imaginary part, so that log(−z) = log z + iπ. The corresponding expressions for the crossed double box diagram, which is obtained from the double box by exchanging p 1 ↔ p 4 , can be retrieved via analytic continuation in the same manner as described for the crossed box. Applying therefore (3.13) we find

Non-planar double box
In this case the relevant soft integrals arẽ wherẽ The Feynman prescription −i0 is left implicit. Aside from some different sign choices compared to the planar case, an important novelty is thatρ 4 now depends on both ℓ 1 and ℓ 2 . A pure basis of master integrals is given by This family of master integrals embeds the nonplanar double box as in Fig. 4. The scalar non-planar double box integral can be decomposed as follows in the soft limit, and The master integrals satisfy the following differential equations The even-and odd-q systems once again decouple and we can write where the matrices are given by Most master integrals can be obtained from the ones discussed in the case of the III topology. Indeed, some integrals of the IX family are identical to their III counterparts, possibly up to the overall sign, while others are related to them by crossing symmetry, here denoted by "↔", Crossing symmetry can be effectively implemented by the replacement The only genuinely new integrals are therefore f IX,10 , for the even sector, and f IX,14 , for the odd sector, whose boundary conditions are given by More details about this point are provided in Appendix A. The first few terms of f IX,10 in the ǫ expansion, obtained by solving the differential equations perturbatively, are as follows: On the other hand f IX,14 vanishes up to O(ǫ 3 ) included. As before, in the above expressions x is implicitly understood to have a small negative imaginary part, so that log(−x) = log x + iπ.
The final result for the decomposition coefficients defined in (3.48) is The corresponding expressions for the crossed non-planar double box diagram can be retrieved applying (3.13), which yields We have checked, where applicable, that our expressions for planar and non-planar double box integrals agree with the soft limit of the exact results in [65,66] and [67].

H topology
Considering the same family of integrals defined in Eq. (3.15), a pure basis of master integrals for the H topology obtains letting, for the even sector, f H,10 = ǫ 4 q 2 G −1,1,1,−1,1,1,1,1,1 + 1 2 ǫ 2 (2ǫ − 1)G 0,0,0,0,1,1,0,1,1 For the odd sector, we used the software epsilon [68] to derive the following pure basis of master integrals, The scalar H diagram is given in particular as in fig. 5 by up to subleading corrections in q. The differential equations read in this case Let us note that some integrals coincide with those appearing in Section 3.2, We omit the boundary condition for f H,13 , which we shall not need in the following, while postponing the discussion of the static limit of f H,14 to its cut counterpart. For the remainder of this subsection, we shall focus on even master integrals only. Eqs. (3.78a), (3.78b), (3.78c) and (3.79a) hold for any x, and are independent of it.
(3.85) up to subleading orders in q. Applying Eq. (3.13) we also obtain the corresponding result for the crossed topology, (3.86) Results (3.85) and (3.86) agree in particular with the exact expressions obtained in Ref. [69] for the equal mass case without expanding in q. The real part of I H + I H also agrees with the result obtained via velocity resummation of the potential-region expansion [37,13].

Cut topologies
In the calculation of inclusive quantities associated to the three-particle cut, which will be discussed in Section 6.1, it is useful to consider master integrals where two massive propagators and one massless propagator are replaced with delta functions according to the so-called cut master integrals. The advantage of this step is to make differential equations and IBP reduction techniques directly available for the calculation of otherwise complicated phase space integrals [58,59,60,61,52]. Indeed, cut integrals can be handled just as the standard integrals discussed in the previous subsections above up to two novelties: • while performing the IBP reduction, integrals where any of the cut propagators cancel out should be dropped at all intermediate steps 6 , • once the corresponding differential equations are obtained, they should be solved imposing appropriate cut boundary conditions.
The needed cut boundary conditions can be directly determined via Cutkosky rules, taking twice the imaginary part of the standard boundary conditions, whenever the integral's topology only involves one cut. When more than one cut is present, we select imaginary parts arising from the singular static region discussed in Appendix A.
(3.89b) These obey the differential equations The new equations can be thus obtained from the old ones by erasing all integrals that do not posses the needed three-particle cut. The only exception to this rule is the differential equation for f cIII,7 , an integral with two distinct three-particle cuts, which acquires an extra factor of 1 2 , as can be seen from the last line of A cIII,0 . The nontrivial boundary conditions read Note that the imaginary part in the boundary condition of f III,10 in eq. (3.30c) arises from the two-particle cut (i.e. from the ordinary static region of Appendix A) and must be discarded. The solutions are, for the even sector, cIII,2 = 0 , f cIII,4 = 2π log x , f cIII,7 = π 2 log x , cIII,4 = 2πLi 2 x 2 + 2π log 2 (x) − and (3.98) For the odd sector, f cIII,9 = π 3 , f cIII,10 = 0 , f cIII,10 = π 3 log(x) . which satisfy the differential equations The boundary conditions can be determined via (3.59), which also hold for cut integrals.

Tree and one loop
We consider the elastic 2 → 2 scattering in N = 8 supergravity where the external states are massive thanks to a non-trivial Kaluza-Klein momentum in the compact directions [31,13]. For later convenience we focus on a s − u symmetric process where the states labelled by p 1 and p 4 in Fig. 1 are dilatons with KK momentum in one compact direction and those labelled by p 2 and p 3 are axions with KK momentum in another orthogonal compact direction. Then the tree-level amplitude reads in terms of the Mandelstam invariants (2.14). From the scattering amplitude in Eq. (4.1) we can get the leading eikonal by going to impact parameter space, following the general strategy outlined in Section 2.2 and using in particular Eq. (2.22). The tree-level amplitude A 0 in the classical limit and the eikonal phase 2δ 0 can be thus cast in the form The one-loop amplitude is instead equal to Inserting the results of Section 3.1 in Eq. (4.4) we get or equivalently in terms of σ, using the relations collected in Section 2.1, The term in the first line of the two previous equations, that we call A 1 , is the super-classical term. Its Fourier transform in impact parameter space reproduces the quadratic term of the expansion of the leading eikonal: The Fourier transform of the term in the second line of Eqs. (4.6) and (4.7), that we call A 1 , gives the sub-leading eikonal: .
Finally, going to impact parameter with the quantum term one gets the quantum part of the eikonal that has a real part given by and an imaginary part given by

Two loops
We combine here the results of Sections 3.2, 3.3 and 3.4 to give the expansion of the two-loop amplitude in powers of q 2 and ǫ. As already emphasized, we restrict to the terms that are non-analytic in q 2 , which are the only ones relevant for the long-range behaviour of the eikonal in b-space. Furthermore, we will neglect terms of order O(q 2 ) 1 2 −2ǫ or smaller, which are irrelevant for determining 2δ 2 and would only enter the calculation of the quantum remainder 2∆ 2 . To this effect, the relevant terms of the s − u symmetric two-loop scattering amplitude [52] are given by: with c(ǫ) as in Eq. (4.5). Inserting the previously computed diagrams we get the following result: The above result can be organized in the following convenient form: up to terms that are further suppressed in q 2 , i.e. O((q 2 ) 1 2 −2ǫ ) or smaller, or that are subleading in ǫ. Note in particular that the coefficient of ((q 2 ) 1+2ǫ ǫ) −1 vanishes identically. A term of O(ǫ 0 (q 2 ) −1−2ǫ ) is also needed in order to fully reproduce the O(δ 3 0 ) exponentiation and will not be discussed further. Similarly, O((q 2 ) − 1 2 −2ǫ ) terms start contributing to order O(ǫ 0 ), although they are omitted for simplicity. We checked that they reproduce the O(δ 0 δ 1 ) exponentiation to leading order. Explicitly, in terms of the variable z and σ defined in Section 2.1, the above coefficients read and Re (4.18)

Eikonal and deflection angle to 3PM order
In Subsection 4.1 we have computed 2δ 0 and 2∆ 1 that can now be used to verify that the eikonal exponentiation takes place as expected up to two loops and to extract 2δ 2 from the previous amplitude through the following relations: Re(2δ 2 ) = ReÃ 2 + 1 6 (2δ 0 ) 3 + 2δ 0 Im 2∆ 1 Im(2δ 2 ) = ImÃ 2 − 2δ 0 Re 2∆ 1 (4.19) where withÃ 2 (s, b) we mean the Fourier transform in impact parameter space of the two-loop amplitude divided by the factor 4pE as in Eq. (2.22). We get : and where we made the branch-cut singularity for z > 1 more explicit via and used Notice that in the ultrarelativistic limit σ ≫ 1 the terms proportional to (σ log(2σ)) 2 present in the second and third last lines separately cancel in agreement with the general pattern discussed in [15].
Let us finally turn to the calculation of the deflection angle χ as outlined in Section 2.2. To 3PM order, we have 2 sin where p is the momentum of either particle in the centre-of-mass frame as in Eq. (2.16).
Moving from b to the orthogonal impact parameter b J , we have 2 tan χ 2 = χ 1 + 1 12 from which we obtain: with (to leading non-vanishing order in ǫ at each PM order): The 1PM contribution corresponds to a tree diagram where both the graviton and the dilaton are exchanged, the 2PM contribution is absent for ǫ = 0 in agreement with the results of Ref. [31] at one loop. The first term in χ 3PM comes from the expansion of tan( χ 2 ) at small χ while the remaining terms are genuine new contributions from Re δ 2 .
Anticipating a more complete discussion in Subsection 5.3 we remark that in the second form of equation (4.20) (already presented together with (4.30) in [15]) the first and last term contain only even powers of p ∼ √ σ 2 − 1 (see (4.32) below), while the second term has only odd powers. This means that those two terms represent halfinteger-PN corrections which are known to be a consequence of dissipative processes (so-called radiation reaction). Instead, the second term in (4.20) corresponds to the more conventional integer-PN expansion due to conservative dynamics.
Looking now at the second form of (4.21) we note that its IR divergent term (the one appearing on the first line) is simply related to the above mentioned half-integer PN terms by a trivial factor − 1 ǫπ . Furthermore, the same half-integer PN terms also multiply, up to a factor 1 π , a log v 2 ∼ log(4(σ 2 − 1)) coefficient. More specifically one can immediately see that the divergent part and the part with log(σ 2 − 1) of Im 2δ 2 are related to the two half-integer PN terms that we call Re 2δ where Re 2δ In the next Section, using general analyticity and unitarity arguments, we will argue that these relations, rather than mere coincidences, should be regarded as generic relations between IR divergent terms in Im 2δ 2 and radiation-reaction corrections to Re 2δ 2 . Since analyticity and crossing are spoiled when we go from the amplitude to its partial waves (the Fourier/Legendre transform introduces spurious kinematical singularities and refers to one specific channel), we have to go back to the amplitude itself.

Real-Analytic, Crossing-Symmetric Formulation
In the previous subsection we have collected the contribution of the different diagrams both at the one and at the two-loop level, noticing that they combine into much simpler expressions than those of individual diagrams. Those expressions are typically in the form of an expansion in powers of q 2 and, for the two-loop calculation, of ǫ.
Furthermore, as a computer output, they simply collect independently real and imaginary terms and express the result in terms of a choice of the two independent Mandelstam variables, q 2 = −t and s (or of quantities like σ or z which are themselves functions of s and the masses).
Those properties are not at all apparent in the formulae of the previous subsection but they must be hidden somewhere, of course up to the order in q 2 at which we stop our expansion. Here we will explicitly show how to rewrite the tree, one-loop and two-loop amplitudes, given in the previous section in a real analytic crossing-symmetric form. Besides its usefulness for the interpretation of the result this also serves as a rather stringent test of the calculations themselves.
To this purpose it is useful to introduce variables which are related by s ↔ u exchange to those introduced earlier in Section 2.1:

Tree and one loop
The tree-level amplitude itself can be rewritten trivially as: At one-loop, equation (4.7) suggests trying the following analytic, crossing-symmetric ansatz for the super classical term: With the help of the following relation it is easy to check that, when expanded in q 2 using (5.1), (5.3) reproduces both the super-classical (first line in (4.7)) and the O(ǫ 0 m 1 m 2 ) terms in the square bracket of (4.7). The remaining terms can be described as follows: • The classical contribution proportional to (q 2 ) −1/2 (second line in (4.7)). It can be trivially symmetrized in s − u since the error in doing so is of order (q 2 ) 1/2 ; • A quantum contribution of O(ǫ) which can be written in the crossing-symmetric form: Note that, in (5.5) and (5.6), we can again neglect the difference betweenσ and −σ as well as the difference betweenz and −1/z. In order to facilitate the discussion of the various terms in Eqs. (5.11b)-(5.11f) we give here the Fourier transform in the impact parameter space (denoted as usual with a tilde) of the two previous quantum terms. We get

Two loops
Moving now to two loops we combine the different contributions in the known form of Eq. (4.13), which suggests to try the ansatz: In order to reproduce Eqs. (4.15) and (4.16) we find: and for Eqs. (4.17) and (4.18): The two previous equations have been derived using log(−z) = iπ+log z and the following identities and branch choices: Combining these with (4.22) we can also write the combination appearing in (5.11c) in a form: that is useful to discuss the non relativistic (σ, z → 1) limit. We can check that the expressions (5.10) and (5.11) satisfy crossing symmetry and real-analyticity. The first property follows simply from the fact that, at this order in q 2 , we can use the identificationz = − 1 z . Checking real analyticity is a bit more subtle. For (5.10) one needs to take into account the subleading (in q 2 ) term origination from the first line. This (purely real) term exactly cancels a similar term coming from the second line. Then one is left with a purely imaginary term from the second line: that can be written in real-analytic crossing-symmetric form as (log(−z) + log(−z))(log(z 2 ) − log(z 2 )) (5.15) At this point checking that the amplitude is real in the unphysical region σ 2 < 1 is straightforward once one realizes that, in that region, |z| 2 = 1, i.e. z is a pure phase and thus log z is purely imaginary. This implies that the term in Eq. (5.14) is real and the same is true for the last term in (5.10) because of an extra factor i coming from √ σ 2 − 1. Coming now to (5.11), real analyticity is easily checked along the same lines for (5.11a), (5.11b), (5.11d) and (5.11e). Concerning instead (5.11c) and (5.11f) one has to remember that, in the unphysical region, z −1 =z. Since both Li 2 (z 2 ) and log(1 − z 2 ) are real-analytic functions, the combinations involving them, appearing in (5.11c) and (5.11f) give a purely imaginary and purely real factor, respectively. Because of the different powers of (σ 2 − 1) appearing in the two contributions this is exactly as needed for real analyticity. In conclusion, we have shown that Eq. (5.10) and Eqs. (5.11a) . . . (5.11f) are real analytic functions. Note that this would not have been the case for the eikonal itself.
Before going further, in Subsection 5.3, into considerations of analyticity let us discuss, with the help of the Fourier transform of the pre-factor in Eq. (5.9) how the different terms in (5.11) have precise correspondences in the eikonal 7 : • Eq. (5.11a) is a purely real term. It is easy to check that its Fourier transform matches exactly the combination −4δ 0 Im ∆ 1 if we insert for Im ∆ 1 the corresponding term given by the Fourier transform of (5.6) in (5.8). Therefore (5.11a) does not contribute to Re δ 2 ; • Eq. (5.11b) matches only half of the corresponding contribution from 4i Re ∆ 1 δ 0 if we insert for Re ∆ 1 the second term in (5.7). This leaves an identical contribution in Im δ 2 .
• Eq. (5.11c) provides a new logs-enhanced contribution to be discussed below. It corresponds to the last term in Eq. (4.20); • Eq. (5.11d) gives a log(s)-enhanced contribution whose real part can be checked to match the result for Re δ 2 given in ref. [13]. It corresponds to the second term in Eq. (4.20); • Eq. (5.11e) This contribution has both a real and an imaginary part. The real part, like with (5.11a), matches exactly a corresponding term in −4δ 0 Im ∆ 1 (the imaginary term in (5.7)) leaving behind no contribution to Re δ 2 . Similarly, the imaginary part matches 4iδ 0 Re ∆ 1 (the term with log z in (5.7)) leaving no residual contribution to Im δ 2 ; 7 Comparison among the various quantities is performed up to terms that vanish when ǫ → 0 • Eq. (5.11f) is a new non log-enhanced contribution reproducing at high energy the ACV90 result. It corresponds to the first term in Eq. (4.20).
In summary, for what concerns Re δ 2 we are left with the result given in [13] accompanied by two new contributions given by Eq. (5.11c) and Eq. (5.11f). It is thanks to these new contributions, given already in Eq. (4.20), that the high-energy limit of Re δ 2 is finite, universal, and in agreement with [32]. The two above-mentioned new contributions to Re δ 2 have some distinctive features, in particular with respect to their PN expansion discussed in the Appendix B.

Connecting real and imaginary parts via analyticity
Before turning to the discussion of the analytic properties of our results (5.10) and (5.11) let us look at the imaginary parts they contain picking up, in particular, those containing, in (5.11), a log(1 − z) ∼ 1 2 log(σ 2 − 1) enhancement for σ → 1. The result is: After multiplying both equations by the appropriate power of ǫ we note that they fail to reproduce the combination we had noticed at the end of Subsection 4.3 for Im δ 2 . More specifically, the contribution of (5.11f) is accompanied by twice as large a contribution from (5.10). At the same time the contribution of (5.11c) has the correct counterpart in (5.10), but only for its −2σ piece in the last term of (5.17) while the σ 3 piece is missing. Fortunately, this apparent contradiction is completely resolved once we subtract from Im A 2 those pieces that correspond to the elastic cut of the amplitude itself (note that this subtraction only affects Eq. (5.10)) given by the Fourier transform from impact parameter to momentum space of the quantity 2δ 0 2 Re ∆ 1 4m 1 m 2 √ σ 2 − 1. In other words, the full analyticity argument has to be carried out on the (analytic and crossing-symmetric) piece of the two-loop amplitude whose imaginary part (entering suitably subtracted dispersion relations as discussed below) is restricted to inelastic contributions. Therefore, indicating byÂ (3pc) 2 that piece of the two-loop amplitude, we find: in full agreement with the structure (5.18).
The appearance of that combination is by no means an accident and can be understood [49] as following from unitarity and three-particle phase space considerations. In D = 4 − 2ǫ one finds that the imaginary part of the three-particle cut behaves, for σ → 1 as where ω is the energy of the massless particle and we have used the fact that the relevant upper limit in ω is O(v/b). We refer to [49] for further details.
Let us now look at (5.11) which corresponds to O(ǫ 0 ) contributions to δ 2 . Let us look, in particular, at the new contributions appearing in (5.11c) and (5.11f). Leaving the discussion of their PN expansion to Appendix B we note that, in both, the real parts we are interested in are accompanied by an imaginary piece that contains an overall log(σ 2 − 1) factor. More precisely, in the s-channel physical region real and imaginary parts appear in the combination: where we have used the relation: showing that the position of the branch points correspond to thresholds for the s and u-channels.
On the other hand we have seen that, combining real and imaginary parts, the combination (5.20) appears giving, in total, the structure: in full agreement with (and in addition to) what we have observed at the δ 2 level. The combination appearing in (5.23) looks to be forced upon by general considerations of unitarity, analyticity and crossing. Unitarity determines the behaviour of the imaginary parts of the amplitude near the branch points at σ = ±1. for the elastic and inelastic (i.e. three-particle) channels. The elastic channel leads to a square-root type branch point, while for the inelastic channel integration over the soft-particle momentum induces the logarithmic enhancement ∼ log(σ 2 − 1) of (5.20).
Concerning the implications of analyticity it is convenient to strip out of the amplitude a harmless (σ 4 +σ 4 ) as in (5.9) and meromorphic factors such as σ 2 /(σ 2 − 1). We can then write an unsubtracted dispersion relation for the resulting "reduced"Â (1) 2 amplitude. Also, because of s ↔ u symmetry, the left and right hand cut contributions combine to give the familiar form: 24) where P denotes the principal value of the integral. Finally, inserting in (5.24) each contribution to the imaginary part generates automatically the corresponding real part given in (5.11). In other words, the full real-analytic amplitude is uniquely reconstructed from its discontinuities across the physical cuts via a Cauchy integral going around the branch cuts and a circle at infinity. In particular a log(σ 2 − 1) factor in the imaginary part induces, via (5.24), a corresponding factor +π in the real part, in agreement with (5.23).
The above general considerations are neatly illustrated by considering, as an example, two different imaginary parts, one with and one without the log(σ 2 −1) factor, see (5.11b) and (5.11f). The above mentioned Cauchy integrals can be computed explicitly and give, at σ 2 > 1 and above the cut, showing the appearance of the extra π-term in the real part for the former case. The above considerations lead us to conjecture that the combination (5.23) is a generic feature for gravitational theories containing, besides a graviton, other massless degrees of freedom. This conjecture is indeed fully confirmed [49] by the GR example, for which we recover the recent result of [47] as well as by the N = 8 supergravity example at different values of the cos φ parameter. It is also confirmed by the explicit calculations of the three-particle cuts given in Section 6.

Direct-Unitarity Calculations
So far we focused on massive N = 8 supergravity where the 4-point amplitudes take a particularly simple form, see (4.2), (4.4) and (4.12). We then considered the softregion expansion, which is the one relevant for the classical limit, and used the results of Section 3 to evaluate the resulting integrals. This allowed us to extract the complete eikonal δ up to 3PM order by applying (2.23). When considering the contributions of the full soft region, one obtains both a real and an imaginary part for δ 2 , the 3PM term of the eikonal. The complete N = 8 results display an intriguing relation between the IR-divergent O(ǫ −1 ) term of Im 2δ 2 and the radiation-reaction terms Re 2δ (rr) 2 of the real part, see (4.31). In the previous section, we argued that this relation is in fact a general feature of gravitational field theories in the two derivative approximation, including GR. The main goal of this section is to provide further evidence for this by deriving both the O(ǫ −1 ) and O(ǫ 0 ) terms of Im 2δ 2 for GR and showing that they indeed display the structure appearing in (5.23) with a coefficient that agrees with the radiation-reaction result of [47,49]. Figure 6: Kinematic conventions for the three-particle cut. The lines in bold represent classical states. The drawing inside the dashed bubbles does not represent a specific topology, but just provides a visual help to recall the definition of the q i variables.
We will take a direct approach based on the 3-particle unitarity cut, which can be evaluated starting from the tree-level 2 → 3 inelastic amplitude A 5 and taking its "square" to obtain the corresponding contribution to the imaginary part of the elastic two-loop amplitude, schematically where k and k 1,2 are the momenta of the massless and the two massive particles involved in the cut, while p 1,2 denote the massive external momenta, as in Figure 6. We will specify more precisely how |A 5 | 2 must be understood in the next subsection. Starting from (6.1), after taking the classical limit, a Fourier transform to impact-parameter space directly yields Im 2δ 2 , which can be interpreted as the total number of gravitons emitted in collisions at given energy and impact parameter. Moreover, multiplying the integrand in (6.1) by k µ , the very same steps lead instead to the total emitted energy-momentum for such processes. The approach we adopt is thus complementary to the one of [52] as far as the integrand construction is concerned. In our approach we determine the integrand of (6.1) using the 5-point amplitude in (6.5), while, in the latter reference, the 3-particle cut is extracted from the 2 → 2 elastic amplitude. However, once the integrand is obtained, we then evaluate the phase-space integrals by trading the mass-shell delta functions for cut propagators and reducing the integrations to the cut master integrals discussed in Section 3.5, borrowing ideas and techniques from the reverse unitarity approach [58,59,60,61,52].
Our strategy allows us to treat both the N = 8 case and pure GR essentially on the same ground, bypassing the explicit form of the full 2-loop integrand in GR. 8 In the end, we can compare the result for Im 2δ 2 in N = 8 supergravity with (4.21) finding perfect agreement, while the GR result is new, provides a confirmation of the general considerations of Section 5.3, and can further serve as a crosscheck for a full derivation of the 3PM eikonal in that case. Performing the analogous calculations for the radiated energy and momentum, we also reproduce the results of [52].
Of course it is often interesting to consider differential (as opposed to inclusive) quantities such as the radiated energy per unit of frequency and/or per solid angle. In this case the direct approach taken here is probably the only one available since it allows one to obtain the full transition rate |A 5 | 2 without dropping any of its pieces at intermediate steps, not even scaleless terms that cancel out in the full integral (6.1). Armed with the complete integrand one can then integrate (numerically if needed) on any of the remaining independent variables, obtaining for instance the graviton emission spectrum or the emitted power spectrum. Also, the remaining quantities can be expressed by using different variables depending on what is most convenient. For instance, one can take a (D − 2)-dimensional Fourier transform to impact parameter space of A 5 and obtain the waveforms in terms of the orbital angular momentum and the frequency. This is what was done in [49], albeit only for the leading soft contribution, i.e. the leading term of A 5 as k → 0. If one starts from the full classical amplitude A 5 , as for instance discussed in [53] for the non-relativistic limit and recently done exactly in [55], then instead of the simple rational dependence on b, the waveforms involve Bessel functions, which provides an exponential suppression for very large frequencies. Alternatively, it is possible to perform a (D − 1)-dimensional Fourier transform and write the waveforms in terms of the retarded time and recently refs. [29,55] showed how to derive the results of [53] starting directly from a Feynman diagram approach.

Direct unitarity in momentum space
In this subsection we focus on inclusive quantities, such as the imaginary part of the 3PM eikonal or the total radiated energy-momentum to 3PM order, as discussed above.
As a starting point to calculate the five-point amplitude A µν 5 (p 1 , p 2 ; k 1 , k 2 , k) on the left in Figure 6, we use the 2 → 3 tree-level amplitude in Eq. (3.1) of [15], which had been obtained from the low energy limit of a string amplitude in a toroidal compactification. Thus, in particular, the theory naturally includes the contribution of the dilaton, together with vectors and scalars arising in the Kaluza-Klein compactification. We then enforce the classical limit: using the conventions summarised in Fig. 6, so that we take the momenta q 1,2 and k to be simultaneously small, of the order of the elastic momentum transfer q (see (2.20)), Equivalently, reinstating momentarily , this can be regarded as a formal → 0 limit in which the wavenumbers q 1,2 / and k/ are held fixed. The particles corresponding to the bold lines are classical and so their energy and momentum are instead independent of q, to leading order, To include contributions stemming from the graviton and the dilaton, together with Kaluza-Klein vectors and scalars, simultaneously, we find it convenient to once again promote all vectors p µ 1,2 , k µ 1,2 , q µ 1,2 , and k µ to 10-dimensional vectors P M 1,2 , K M 1,2 , q M 1,2 , and k M . Their precise relation to the effective four-dimensional kinematics depends on the theory under consideration and will be specified below.
The leading contribution to A 5 in the classical limit is of order O(q −2 ) and reproduces the result of [20,21,29]. It can be written in the following convenient form , while the quantity β is defined in (6.8) below depending on the theory under consideration. The main feature of (6.5) is that it sat- = 0 for arbitrary values of the free index, which makes the calculations in general dimensions easier [12,70].
Let us now describe how the N = 8 supergravity and GR amplitudes are "squared" in order to construct the integrand in (6.1). In both cases, the kinematics is given as in Figure 6 by and q 1 + q 4 = P 1 + P 4 = q , q 2 + q 3 = P 2 + P 3 = −q . (6.10) In the maximally supersymmetric case, we simply work in 10D and saturate the Lorentz indices with the Minkowski metric so as to get where we do not need to subtract the contribution of the longitudinally polarised states as the amplitude A 5 is transverse. As already stressed, the advantage of this strategy is that, doing so, one effectively encompasses all the contributions that were taken into account separately in [49], i.e. besides the graviton and the dilaton, also the vectors and scalars arising from the toroidal compactification. For GR all indices become 4D and on the top of this we need to subtract the contribution of the dilaton. Thus we have where the structure within square brackets arises from the sum over the physical polarisations and is the one appearing also in the de Donder propagator, simplified by using the symmetry of the A 5 in its two Lorentz indices. Once the integrand |A 5 | 2 of (6.1) is constructed according to the above strategy, one then needs to perform the phase space integrals. This task can be greatly simplified by resorting to the reverse-unitarity strategy [58,59,60,61,52]: trading the mass-shell delta functions appearing in the second line of (6.1) with the corresponding cut propagators, one can recast the integral as an effectively more tractable (cut) two-loop integral. After performing these steps, we need to decompose the resulting integrand into structures that match the cut integrals evaluated in Section 3.5, that is, to isolate terms having each the appropriate denominator structure for a given topology once a suitable identification of the integrated momenta is applied. This decomposition can be performed in several equivalent ways associated to different possible polynomial divisions.
A simple-minded approach is to just expand out all terms in the product |A 5 | 2 and use their denominators written in terms of p 1 k, p 2 k and q 2 1 , q 2 2 , q 2 3 , q 2 4 to group them into the following classes: where the massive propagators are in general raised to non-negative integer powers. These structures can be associated to the topologies of Figure 7. With proper identifications of q 1 and q 2 in terms of ℓ 1 , ℓ 2 and q, all such structures can be embedded in the cH or cIX families of Section 3.5, possibly after use of crossing symmetry and of the u 1 ↔ u 2 symmetry. For example, letting the first structure of Eq. (6.13) can be recognized as belonging to the cH family: inserting the cut propagators, denoted with an underscore, Moreover, letting we find, for the first structure in the second line of (6.13), 2 , (6.17) which we recognize as belonging to the cIX family up to an (immaterial) interchange of u 1 and u 2 and the exploitation of crossing symmetry. After the appropriate family is identified, the integral under scrutiny can be reduced to the corresponding cut master integrals discussed in Section 3.5 using LiteRed. When expanding the mushroom 7b and "IY" 7d topologies in terms of the cIX family, however, one should note that, again via crossing symmetry,G 0,0,k 1 ,k 2 ,1,1,0,0,1 ↔ G 0,k 2 ,k 1 ,0,0,1,1,0,1 (6.18) which relates integrals of the typeG 0,0,k 1 ,k 2 ,1,1,0,0,1 written in terms of the cIX family to integrals G 0,k 2 ,k 1 ,0,0,1,1,0,1 appearing in the definition of cH masters, f cH,8 for the even sector and f cH,12 for the odd sector.
A slightly more systematic approach is to start out by fixing a particular choice of independent momenta, for instance thus rewriting the integrand in the form N(u 1 ℓ 2 , u 2 ℓ 2 , ℓ 1 q, ℓ 2 q, u 1 ℓ 1 , u 2 ℓ 1 , ℓ 2 1 , ℓ 1 ℓ 2 , ℓ 2 2 ) , (6.20) where N is a Lorentz-invariant numerator and the denominators are In particular, the powers α 1 and α 2 in (6.20) take the value 2 in N = 8 supergravity and 4 in GR. One can then apply simple integrand-reduction techniques using basic tools from algebraic geometry such as Groebner basis and syzygy relations (see e.g. [71]). To simplify matters, one can first trade N with the remainder R of the polynomial division between N itself and the Groebner basis of D 7 , D 8 and D 9 , thus eliminating many terms that vanish on-shell. Denoting byg 7 ,g 8 ,g 9 the Groebner basis of D 7 , D 8 and D 9 , this means for some quotient polynomialsq 7 ,q 8 ,q 9 and a suitable 3 × 3 conversion matrixã ij . Then, one can perform the polynomial division between R and the Groebner basis of the full set of denominators, thus generating a number of integrands with simpler sets of denominators, determined by the corresponding quotients. One can also check that the remainder of this division is zero. Explicitly, letting g i with i = 1, . . . , 9 denote the Groebner basis of D 1 , . . . , D 9 , with conversion matrix a ij such that one then finds quotient polynomials q i such that and therefore 25) up to terms where the cut propagators would again cancel out. One can automatise this procedure and repeat it for each of the newly generated ratios, collecting along the way those terms for which the remainder of the polynomial division is nonzero, which constitute the irreducible integrands. Moreover, the form of the numerators generated at each step can be further simplified by using the syzygy identities obeyed by the relevant denominators. To implement these steps we used the very convenient software Singular [72] in combination with its Mathematica interface. Proceeding this way, we were able to regroup the integrand into three contributions associated to the cH and cIX topologies of figures 7a and 7d, in the N = 8 case, and into seven contributions associated to the cH, "IY" and cIX of figures 7a, 7c and 7d in the GR case. All integrals can be then mapped to the cIII and cIX families via appropriate changes of variables and judicious use of the u 1 ↔ u 2 symmetry. For instance, for the denominator we can perform the change of variables ℓ 1 → −ℓ 1 − ℓ 2 − q, followed by ℓ 1 ↔ ℓ 2 and by the harmless replacement u 1 ↔ u 2 to obtain the denominator 1 (u 2 ℓ 1 )ℓ 2 2 (ℓ 1 + ℓ 2 + q) 2 (ℓ 2 + q) 2 ℓ 1 2 u 1 ℓ 2 u 2 (ℓ 1 + ℓ 2 ) , (6.27) which is appropriate for the cIX family, up to crossing symmetry. The observation (6.18) again comes in handy when performing the reduction of integrals associated with figure 7c. We include in an ancillary file the result of this decomposition of the N = 8 and the GR integrands.
In the calculation of the three-particle cut contribution to the imaginary part of A 2 , Eq. (6.1), only the even sectors of master integrals are actually involved. As already mentioned, the 3-particle cut is the only and full contribution to Im 2δ 2 , so, in the case of N = 8, evaluating (6.1) and taking the (D − 2)-dimensional Fourier transform to impact parameter space according to (2.22), we reproduce Eq. (4.21). For the case of GR we obtain instead the new result Im 2δ (6.28) The first line agrees with the result obtained in [49] by using the soft approximation which is sufficient to capture the IR-divergent term. The second line matches the analytic behaviour as σ → 1 derived on general grounds in Section 5.3. The ultrarelativistic limit σ → ∞ reproduces the universal behaviour derived in [15] by using an argument based on analyticity/crossing and an explicit calculation in the double Regge limit. In particular the terms proportional to (σ log(2σ)) 2 in the second and third line cancel and the leading term is actually (2σ) 2 log(2σ) in agreement with [15]. Combining the conservative eikonal phase that can be extracted from the scattering angle given in [36,37] with the radiationreaction contribution derived in [49] and confirmed by (6.28), let us complement (6.28) with the corresponding (IR-finite) real part: Re 2δ (6.29) in order to match as much as possible the two results. For massive N = 8 supergravity we reproduce precisely f 1 and f 3 of Eq. (6.31) but we get f 2 = 0. In the case of GR we get: matching the results for f 1 and f 3 in Eqs. (6.32) only for σ → 1 + . The vanishing of f 2 follows, in both cases, from the fact that this contribution comes from diagrams where the massless particle is emitted from an internal leg and thus do not contribute to Weinberg's soft limit. Let us also mention that, while the approach developed here directly captures the classical result (6.30), the analogous calculation based on taking suitable cuts of the elastic two-loop amplitude relies on the cancellation of spurious superclassical contributions O(b −1 ), O(b −2 ) which appear at intermediate steps.

Direct unitarity in impact-parameter space
We now consider an alternative approach to the one of the previous subsection. This approach was adopted in [49], even if in that paper the analysis was restricted to the Weinberg limit [73] of the inelastic amplitude A 5 , i.e. to the leading term for k ≪ q i . The basic idea is to introduce the Fourier transform to impact parameter space at an earlier stage in order to rewrite Eq. (6.1) as the product of two inelastic amplitudesÃ 5 expressed in terms of the impact parameter b, which is then integrated over the momentum k of the massless particle in the cut. Most importantly, the integrand obtained by this procedure represents the number of produced quanta per unit of phase space for each massless particle species and polarisation.
As a first step we write Eq. (6.1) for the three-particle cut in the explicit form where k 0 1 , k 0 2 , k 0 are the energies of the three intermediate particles and we are using the kinematics of Eqs. (6.6) or (6.7) depending on whether we consider N = 8 or GR. The sum over physical polarisations in the square brackets reproduces the structures in (6.11) or (6.12) again depending on the case under consideration.
To fix ideas, it is useful to further specialize the kinematics (2.3) of the elastic process, separating transverse and longitudinal components with respect to the classical direction of propagation: Here and in the following, where needed, transverse (D − 2)-vectors shall be denoted by boldface letters. Conversely, components along the (D − 1)th direction, referred to as longitudinal, shall be labelled by a superscript L. On the other hand, spatial (D − 1)vectors shall be denoted with an arrow. For instance, for the q 1,2 defined in (6.2) and the massless momentum k, In Eq. (6.35), the integral over the longitudinal components k L 1,2 of k 1,2 can be performed using the longitudinal and energy components of the δ-function of momentum conservation, explicitly The resulting Jacobian cancels the factors of k 0 1,2 in the measure and produces an extra factor of |k 0 1 k L 2 − k 0 2 k L 1 | −1 . We are then left with only the integrals over the remaining D − 2 transverse directions: where we changed variables in the integration to the q 1,2 introduced in (6.2) (see also Fig. 6 for a visual help). In the classical limit we can safely approximate k L 1 ≃p ≃ p and k L 2 ≃ −p ≃ −p, where p is the absolute value of the momentum in the centre of mass frame. Similarly k 0 1 ≃ E 1 and k 0 2 ≃ E 2 , so that, in terms of the total energy E = √ s = E 1 + E 2 in the centre of mass frame, we have 4|k 0 1 k L 2 − k 0 2 k L 1 | ≃ 4Ep. This is precisely the factor appearing in the Fourier transform to impact parameter space (2.22).
In order to treat the two 5-point amplitudes in Eq. (6.39) more symmetrically, let us impose the relations in (6.10) by introducing integrals over q 3 and q 4 and two extra δ-functions. This can be done inserting the factor in Eq. (6.39). Going to impact parameter space then gives where we used q = 1 2 (q 1 − q 2 + q 4 − q 3 ) in the exponent. As reviewed in Section 2, Re 2δ 2 is directly related to the classical deflection angle, while, as will become clear shortly, 2 Im 2δ 2 is equal to the total number of massless states that are produced in a collision of the two massive particles, for fixed energy and impact parameter, to order G 3 . The form obtained in Eq. (6.41) suggests also to focus on each one of the square parenthesis separately, which represents, for each i, the corresponding wave-form 10 in terms of its frequency and angular direction. Defining M N we are thus lead to consider: where ∆ ≡ 1 2 (q 1 − q 2 ), and we used the delta function to carry out the integration over q 1 + q 2 , which sets q 1 + q 2 = −k . (6.43) We display only the dependence on b and k after the Fourier transform, even if, of course,Ã depends also on the masses and the centre of mass energy. At leading order in the Weinberg limit k ≪ q 1,2 only the last line of (6.5) survives and, since in this limit q 2 2 ≃ q 2 1 ≃ q 2 1 , the Fourier transform can be easily performed and one obtains the results summarised in Section 4 of [49]. As discussed in that reference, this limit is sufficient to obtain the IR divergent part of Im 2δ 2 , the zero-frequency limit of the radiated energy, and the radiation reaction on Re 2δ 2 .
One can go beyond this approximation by keeping k to be generically of the same order as q 1 , q 2 as in (6.3) and then obtain the full spectrum for each massless particle in the classical limit. In the rest of this section we will illustrate this procedure for the case of the dilaton spectrum in N = 8 supergravity.
Let us start from the classical amplitude A 5,dil for dilaton emission that is obtained from Eq. (6.5) by taking β = 4(p 1 p 2 ) 2 and saturating it with the dilaton projector (6.44) 10 In particular, for i = +, ×, it will give the two linear polarizations of the gravitational wave.
Using in particular the exact relations which follow from the mass-shell constraints, and neglecting irrelevant terms in the classical limit, we obtain (6.46) Let us now introduce rapidity variables y 1 , y 2 and y according to where y 1 > 0 , y 2 < 0 , (6.48) andm 1,2 were introduced in (2.6). Using the conservation conditions and solving first two constraints (6.45) to leading order in the classical limit, so that effectively − p 0 we obtain q L 1 = q 0 1 coth y 1 = |k| sinh y tanh y 2 − cosh y tanh y 1 − tanh y 2 , q L 2 = q 0 2 coth y 2 = |k| cosh y − sinh y tanh y 1 tanh y 1 − tanh y 2 , (6.51) so that q 2 1 = (q 1 ) 2 + c 2 1 k 2 , c 1 = cosh(y − y 2 ) sinh(y 1 − y 2 ) , q 2 2 = (q 2 ) 2 + c 2 2 k 2 , c 2 = cosh(y 1 − y) sinh(y 1 − y 2 ) . (6.52) Note, however, that the expressions in terms of rapidities hold in a general frame boosted in the longitudinal direction (corresponding to an overall shift of all the rapidities). Consider now, for fixed k, the Fourier transform of (6.46) rewritten in the form: where each B (a) corresponds to a different scaling with k at fixed q i . The basic objects to be computed are then:B Let us remark once again that, at the classical level, one can neglect corrections of order q i , k approximating for instance −p 1 ≃p. The simplest case isB (2) , which reduces to the standard Bessel integral using equations (6.43) and (6.52) to rewrite q 2 1 = 1 2 (q 1 − q 2 ) + 1 2 (q 1 + q 2 ) 2 + c 2 1 k 2 = ∆ − 1 2 k 2 + c 2 1 k 2 (6.63) (and similarly for q 2 2 ), and shifting ∆ appropriately. We thus find, restricting for simplicity to D = 4, 11 The case ofB (1) is slightly more involved because of the kq i appearing in the numerators. Using (6.53) we can separate the longitudinal part, for which we can again apply (6.62) while for the transverse part we need to replace q i by the appropriate derivative with respect to b. Using K ′ 0 = −K 1 , we get: Let us finally turn toB (3) . The most convenient way to compute it is by introducing two Schwinger parameters t 1 , t 2 in order to write the two propagators in exponential form. One then performs the Gaussian integral over ∆ followed by the integral over T = t 1 + t 2 that can be done in terms of a K 1 Bessel function. One is left with an integral over x ≡ t 1 /T . The final result is: We may further simplify (6.64), (6.65), (6.66) by replacing the remaining scalar products in terms of masses and rapidities. The outcome is that eachB depends only upon the products m 1 m 2 and |k|b, the rapidities y, y 1 , y 2 , and the azimuthal angle φ between b and k.
Assembling together the various terms and prefactors we arrive at the following result: Here, c i , d i and f are defined, respectively, in (6.52), (6.53) and (6.66). The above result uses the transverse momentum k and rapidity y to represent the dilation's momentum. Note that M dil is dimensionless and, when written in these variables, it is also independent of the masses. Depending on circumstances it can be more convenient to express(6.68) in terms of the centre of mass energy ω and angles θ, φ of the produced dilaton. One finds: where one should insert the relations given in equations (6.58), (6.59), (6.57) and write explicitly the different scalar products. The amplitude in impact parameter space (6.69), appropriately squared as in (6.41), determines the differential spectrum of the number of emitted dilatons according to Equivalently, using the relation d 3 k = d 2 kd(|k| sinh y) = d 2 k|k| cosh ydy and the fact that ω = |k| cosh y, one rewrite dN dil in the form with M dil given in (6.68). Other useful relations that allow one to recast dN dil in different forms are (6.57), (6.58), (6.59), together with It is easy to check that, if we interpret k and ω as a classical wave-vector and a frequency, respectively, (6.70) must have an extra factor −1 on its right-hand side. Indeed the number of emitted dilatons per unit of phase space diverges in the → 0 limit. On the other hand the dilaton's energy emission spectrum is simply obtained multiplying dN dil by ω, and has therefore a smooth classical limit. One can also compute less differential quantities, such as the frequency spectrum dN dil dω obtained after angular integration or the angular distribution dN dil d 2 Ω integrated over the frequency. It should be clear how the above considerations can be applied to the other massless particles produced in the collision, in particular to the graviton in GR, where we can also disentangle the two polarizations of the gravitational wave. Furthermore, in principle, evaluating the integrals of our spectra corresponding to the total number or total energy of the emitted particle one ought to reproduce the results derived in the previous section. So far we have been discussing waveforms, as well as number and energy differential rates in the frequency (i.e. ω) domain. Performing, at the level of the amplitude, one more Fourier transform from ω to the retarded time u, one can obtain the corresponding waveforms and rates in the time domain and make contact with other results in the literature [53,54,55].

Discussion and Outlook
In this paper we provided a detailed derivation of the full 3PM eikonal in massive N = 8 supergravity from the two-loop 2 → 2 elastic amplitude. The real part of the eikonal captures the classical deflection angle while its imaginary part is related to the total number of massless particles emitted in the scattering process. This analysis explicitly shows the necessity of taking into account contributions from the full soft region of the loop integrals in order to obtain the correct results for the various physical observables, with a smooth high-energy limit. This is done be using the appropriate boundary conditions of the key scalar integrals entering the relevant differential equations for the classical limit.
We also showed how to rewrite the N = 8 result so as to make the analyticity properties and crossing symmetry manifest. These properties can be used to obtain a direct link between certain terms of the imaginary part and the radiation-reaction contribution to the real part which, in turn, determines the corresponding terms in the deflection angle. The dispersion relation ensuring this relation makes it manifest that the radiation-reaction contributions have a characteristic PN expansion which is shifted by a single power of the velocity v with respect to the contribution arising from the potential region. However, in the ultrarelativistic limit, the two structures scale in the same way and cannot be separated. In fact, as already discussed in [15], there are important cancellations of terms enhanced an by extra-factor of log(s) and this ensures that the deflection angle for the massless case [32] is smoothly reproduced. This 3PM result is a universal property of all gravitational field theories in the two derivative approximation.
In Section 6, the imaginary part of the 3PM eikonal was derived by following an alternative approach where the 2 → 3 tree-level classical amplitude is used to calculate explicitly the 3-particle cut. One can again use the same technical tools developed to evaluate the 2-loop integrals relevant to the elastic two-loop process. In the N = 8 case we reproduce the same result for Im 2δ 2 as the one obtained from the evaluation of the full two-loop amplitude, but this "direct unitarity" approach can be easily extended to obtain Im 2δ 2 , waveforms and fully differential spectra even for the case of pure GR. By following this route we derived the radiation-reaction part of the 3PM eikonal for GR which, when combined with the previous results [36,37] for the potential region contribution, provides the full 3PM eikonal for massive scattering in GR given in (6.28) and (6.29). This fully amplitude-based approach reproduces the physical deflection angle derived in [47] by using the linear response formula of [45].
Of course, a first natural generalisation of our analysis is to use again N = 8 massive supergravity as a laboratory for studying the full soft-region contribution at 4PM. Very recently, the potential contribution in GR was derived [74] and it would be certainly interesting to complete that result by including the radiation-reaction and the tail effects so as to extract a physical, IR-finite result for the deflection angle both in N = 8 and in GR. More in general, the lesson of our 3PM analysis is that analyticity properties and crossing symmetry can be very useful to organise the result and find connections between different terms. On the practical side, this may simplify the derivation of (some parts of) the full 4PM eikonal. On the conceptual side it would be important to extend systematically the discussion of the eikonal exponentiation to the inelastic terms: this presumably requires to promote the eikonal phase to an operator containing creation and annihilation modes for the massless particles that can be emitted in the process, see e.g. [57]. Another interesting development would be to use the direct unitarity approach to continue the analysis of gravitational waveforms, energy spectra and angular momentum emission at 3PM as a step towards an extension to higher orders.
The generalisations mentioned above will also help to provide a more clear physical interpretation of the radiation-reaction contributions to the deflection angle also beyond the 3PM order. Again, working in a maximally supersymmetric setup will provide the perfect arena to uncover the physical principles and technical shortcuts that can then be applied in more general situations. Of course, an ambitious goal would be to provide a complete post-Minkowskian, amplitude-based approach to the analysis of the inspiral phase of gravitational binaries. In the spirit of [26,27], it would be useful not to rely on gauge-dependent intermediate results and obtain the relevant physical observables directly from on-shell scattering amplitudes. This approach has been successful when restricted to the conservative part of the problem, but, as this work hopefully illustrates, radiation effects should be included in order to complete the programme. and with a standard integral over Feynman parameters The integral over t yields √ π 2 while the remaining integration over s 1 and s 2 can be handled in a manner similar to the one in Eq. (A.7), eventually retrieving (3.30a).
The remaining boundary condition for f III,7 can be instead discussed as follows. In view of the prefactor τ 2 , nontrivial contributions to the static limit can only emerge from the leading singularities of the integral G 7 ≡ G 1,1,1,1,1,1,1,0,0 as y → 1. These can be evaluated by first writing the integral as and change variables to 55) In this fashion, β, b parametrize the directions along the two independent zero modes q, k of M(1) and α, a parametrize the orthogonal directions q ⊥ , k ⊥ . Therefore, to leading order small positive τ , we find (A.56) Note the appearance of a τ −2 singularity as τ → 0. Substituting into Eq. (A.51), the resulting integral over t 1 , t 2 and t 3 can be tackled in a manner analogous to f III,6 discussed above. The final result indeed reproduces (3.29g) as the overall factor of τ 2 and the τ −2 singularity cancel out.
Let us now briefly comment on the boundary conditions for the IX family. In view of eqs. (3.59), all master integrals for the IX family can be deduced from the ones discussed for the III family except f IX,10 and f IX, 14 . The derivation of the corresponding boundary conditions (3.61) proceeds in a manner similar to the ones for f III,7 and f III, 10 . The crucial point is that, for these topologies, the singularities arising as τ → 0 are not sufficiently strong in order to compensate the overall factors of τ 2 or τ appearing in their definitions, so that the static limit eventually gives zero.
Moving to the H family, let us first recall that some integrals coincide with those appearing in the III family by Eq. (3.76). The master integrals f H,2 , f H,3 , f H,4 can be directly evaluated via Schwinger parameters. Direct evaluation is also viable for the last three integrals appearing on the right-hand side of the definition (3.68j) of f H, 10 . To discuss the boundary condition for f H, 8 , it is again necessary to distinguish between ordinary and singular static contributions to the small-τ limit of G 0,1,1,0,0,1,1,0,1 . One can verify that the ordinary static contributions cancel between the two terms on the righthand side of the definition (3.68h) of f H, 8 . The corresponding boundary condition (3.78d) is therefore entirely due to the singular static contributions to G 0,1,1,0,0,1,1,0,1 . Finally, the (ordinary) static limit of the integral . (A.58) The discussion of the odd integrals proceeds in a similar fashion. The integral f H,11 can be evaluated directly since it is independent of x, while the boundary condition for f H,12 given in eq. (3.79b) emerges both from the ordinary static region (first line) and from the singular static region (second line).
Note that imaginary parts are all (Z+1/2)PN but the new contributions are distinguished by exhibiting an extra log v enhancement. This makes the corresponding real parts to be also (Z + 1/2)PN. For the usual (Z + 1/2)PN imaginary parts the corresponding real parts are ZPN. We should also mention that in (5.11c) and (5.11f) the leading 0.5PN term cancels both in the real part and in the logarithmically-enhanced imaginary part so that their sum is 1.5PN (it becomes 2.5PN in pure Einstein gravity as discussed in [47] and [49], and after equation (6.29)).