Extremal black hole scattering at O(G^3): graviton dominance, eikonal exponentiation, and differential equations

We use $\mathcal N=8$ supergravity as a toy model for understanding the dynamics of black hole binary systems via the scattering amplitudes approach. We compute the conservative part of the classical scattering angle of two extremal (half-BPS) black holes with minimal charge misalignment at $\mathcal O(G^3)$ using the eikonal approximation and effective field theory, finding agreement between both methods. We construct the massive loop integrands by Kaluza-Klein reduction of the known $D$-dimensional massless integrands. To carry out integration we formulate a novel method for calculating the post-Minkowskian expansion with exact velocity dependence, by solving velocity differential equations for the Feynman integrals subject to modified boundary conditions that isolate conservative contributions from the potential region. Motivated by a recent result for universality in massless scattering, we compare the scattering angle to the result found by Bern et. al. in Einstein gravity and find that they coincide in the high-energy limit, suggesting graviton dominance at this order.


Introduction
The direct detection of gravitational waves at LIGO/VIRGO [1,2] has started an exciting new age of gravitational wave astronomy. Scattering amplitudes have emerged as the latest tool in computing the gravitational dynamics of binary systems in the perturbative regime. In contrast to the traditional post-Newtonian expansion, which is a simultaneous expansion in the Newton constant, G, and a relative velocity, v, relativistic scattering amplitudes can naturally lead to results up to a fixed order in G but to all orders in velocity, known as the post-Minkowskian (PM) expansion . A recent highlight is the result of Bern et al. [19,20] for the conservative dynamics of black hole binary systems at O(G 3 ), i.e. the thirdpost-Minkowskian (3PM) order. The result points to many interesting questions, some of which are explored in the present paper.
to two reasons. First, a direct calculation without a series expansion to high orders can be computationally more efficient. Second, Ref. [47] raised questions about the velocity resummation of Refs. [19,20] in the case of Einstein gravity. Since then, the correctness of the latter result has been verified at high orders in the velocity expansion [48,49], an alternative method for resummation of the velocity series has been used with identical results [50], and the unitarity cut construction of the loop integrand has been checked against direct Feynman diagram computations [51]. Still, a direct relativistic calculation that bypasses velocity resummation will be a valuable additional confirmation of the result, and will provide a way to streamline future calculations at O(G 3 ) and beyond.
The study of classical gravitational scattering in N = 8 supergravity was initiated in a beautiful paper by Caron-Huot and Zahraee [52], which we build upon. The large set of symmetries of this theory provides important simplifications, which make it the perfect theoretical laboratory to study various conceptual questions and test the technology to be applied at higher orders in perturbation theory. This is familiar to how precision QCD practitioners have often sharpened their axes with simpler calculations in N = 4 super-Yang-Mills theory before honing in on the beast. A lot is know about N = 8 supergravity, in particular the complete loop integrands for the quantum four-point amplitude are available through five loops [53][54][55][56][57][58][59][60][61]. These were constructed using the unitarity method [62][63][64][65][66] and the different incarnations of the double copy [67][68][69][70]. These results, being valid in D-dimensions, can be used to easily obtain massive integrands via Kaluza-Klein reduction, as we will do in this paper. 1 Moving on to integration, we will obtain the part of the amplitude relevant for classical conservative dynamics using the method of regions [71]. In particular, integration in the "soft region" produces the correct small momentum transfer expansion of the amplitude [23,72], up to contact terms that are irrelevant for long-range classical physics at any order in G. However, conservative classical dynamics actually arises from the "potential region" which is a sub-region contained in the soft region [16]. Strictly speaking, the potential region is defined in the near-static limit and produces an expansion of the Feynman integrals as a series in small velocity. But since the velocity series can be resummed to all orders, the resummed result will be also referred to as the amplitude evaluated in the potential region. In addition to isolating conservative effects, evaluating in the potential region also simplifies the integrals considerably.
Refs. [19,20] exploits the fact that infrared (IR) divergences cancel when matching the EFT against full theory, and circumvents the evaluation of IR divergent integrals. In this paper, all IR divergent integrals will be evaluated explicitly in dimensional regularization (which serves as both UV and IR regulators). This will allows us to check against the predictions from eikonal exponentiation, which expresses the divergent amplitude in an exponentiated form. Additionally, we will evaluate all integrals relativistically with full dependence on velocity, without constructing and resumming a velocity series. This is made possible by employing the method of differential equations for Feynman integrals [73][74][75][76], with a crucial new ingredient being the use of modified boundary conditions that isolate the contributions from the potential region. While Ref. [20] already presented a precursor of our differential equations method as an alternative to the "expansion-resummation method", it was only successfully applied to a subset of the needed integrals that do not involve infrared divergences due to "iteration" of graviton exchanges. This paper will use a finer control of boundary conditions to evaluate all integrals using differential equations. We also perform soft expansions prior to the construction of differential equations, resulting in dramatic speedups in computation. Another improvement is that we transform the differential equations into Henn's canonical form [77,78]. In this form, the differential equations have a simple analytic structure, and can be easily solved to higher orders in the expansion.
In the context of N = 8 supergravity, Ref. [52] put forward a tantalizing conjecture: that the energy levels of a pair of black holes in such theory retain hydrogen-like degeneracies to all orders in perturbation theory. This is tantamount to the classical black hole binary orbits being integrable and showing no precession. Two pieces of evidence were provided in support of this conjecture: first, the absence of precession for the full O(G 2 ) dynamics, which directly follows from an analog of the "no-triangle" hypothesis [79][80][81][82][83][84] for massive scattering; and second, various all-orders-in-G calculations in the probe limit for different charge configurations. It is known that O(G 3 ) (or any odd power of G) corrections to the conservative dynamics cannot yield precession [85,86]. Instead we will use the scattering angle at O(G 3 ) to test this conjecture, and see that it deviates from the integrable Newtonian result at this order. We will extract the scattering angle both from appropriate derivatives of the "eikonal phase" and via the EFT techniques of Refs. [19,20], finding agreement between both methods.
Although we perform our calculations in N = 8 supergravity, we expect the techniques here developed to be directly applicable to Einstein gravity as well. Such application is beyond the scope of the present paper and we leave it for future work.
This paper is organized as follows: In Section 2 we setup our conventions, we review some basic features of extremal black holes in N = 8 supergravity, and discuss the different limits that will be used in the paper. In Section 3 we construct the tree-level four-point amplitude, as well as the one-and two-loop massive integrands from the known massless integrands via Kaluza-Klein reduction and truncation to the appropriate sector. In Section 4 we briefly discuss the integration regions involved in our problem, and introduce our new integration method based on differential equations, which is applied to calculate the full oneand two-loop amplitudes in the potential region. In Section 5 we assemble the scattering amplitudes. In Section 6 we review the eikonal method, and use it to calculate the order G n≤3 eikonal phase, while checking exponentiation. Then we use the eikonal phase to calculate the gravitational scattering angle and we compare its high-energy limit with the result of Refs. [19,20] in Einstein gravity. In Section 7 we cross check our results via the EFT method of Ref. [16], and we comment on the advantages of this approach. We calculate the conservative Hamiltonian by matching, and find the scattering angle by solving the classical equations of motion. In Section 8 we present our conclusions. We include two appendices: Appendix A contains some technical details about the computation of integrals in the near-static limit and Appendix B collects the solution to our two-loop differential equations. The results are provided in computer-readable format in several ancillary files (see comments at the beginning of each file for detailed descriptions).

Kinematics and setup
We model the dynamics of two half-BPS black holes in N = 8 supergravity [87,88] by considering the scattering of two massive point particles in half-BPS multiplets, which interact via the massless supergravity multiplet. The incoming and outgoing momenta are p 1 , p 2 and p 3 , p 4 respectively, and the masses of the particles are (2.1) We will parametrize the scattering amplitudes in terms of the usual invariants s = (p 1 +p 2 ) 2 , t = (p 1 + p 4 ) 2 = q 2 and u = (p 1 + p 3 ) 2 , where we introduced the four-momentum transfer q = p 1 + p 4 for later convenience. As is common in the study of scattering amplitudes we will cross the incoming particles to the final state, so that all particles are outgoing. The physical scattering configuration corresponds to the region s > (m 1 + m 2 ) 2 , t = q 2 < 0 and u < 0. 2 The half-BPS multiplet in N = 8 supergravity contains massive states with spin 0 ≤ S ≤ 2. In this work we will focus on particular scalar components, φ andφ, with S = 0 and leave the study of spinning states for later work. The interactions between different half-BPS particles are mediated by the massless supergraviton multiplet. In addition to gravitons the N = 8 supergraviton multiplet contains 28 (vector) graviphotons, A IJ , and 70 scalars φ IJKL , as well as fermions which will not be important for our discussion. Black holes in N = 8 supergravity interact with the graviphotons and scalars with charges C IJ given by an 8 × 8 matrix. Here I, J, . . . are SU(8) R-symmetry indices and the vectors a scalars are in SU(8) representations of the appropriate dimension. We will not print the Lagrangian here, because it is lengthy. For our purposes, however, all scattering amplitudes could be built from the three-particle amplitudes: using factorization and unitarity, as done in Ref. [52]. Here ε are polarization vectors. In general the charges, C IJ are complex and the black holes are dyonic. The charges are also central charges of the supersymmetry algebra, and the BPS condition requires their magnitude to be equal to the mass When studying a pair of black holes we need only consider the relative phases in their BPS charges. These are parameterized by three angles along which the charges might be misaligned with Φ = diag(e iφ 1 , e iφ 2 , e iφ 3 , e iφ 4 ) and i φ i = 0. For the two-and one-angle cases, however, there always exist a duality frame where the magnetic charges are zero. We point the reader to Ref. [52] for a more detailed discussion of the charges. Although we will construct the full (quantum) loop integrands for the scattering amplitudes of these black holes, we are ultimately interested in their classical conservative dynamics. In the classical limit of hyperbolic scattering, the orbital angular momentum of the black hole binary system is much larger than . Thus, the classical limit of the scattering amplitudes simply corresponds to the large angular momentum limit J 1 (in natural units), which establishes a hierarchy of scales As a result, we are interested in calculating scattering amplitudes in the limit of small q, or more precisely as an expansion in small q. From a heuristic calculation in the Newtonian limit, the leading-order scattering angle θ is of the order Gm/(vr) ∼ Gmq/v, where m and r are the total mass and relative transverse distance of the system. So for generic values of v, the quantity Gmq is of order θ, and for each additional order of G, we need to expand the amplitude up to one additional power of q to obtain corrections to the scattering angle of order θ L , where L is the loop order. Terms that are more subleading in q at the same power of G are quantum corrections that vanish classically. In summary, at O(G n ), we only need to expand the scattering amplitude of massive particles up to O(|q| n−2 ) in the small-q expansion [72], in order to extract the classical dynamics. In practice this will imply, among other things, that when we calculate an amplitude some loop integrals can be discarded before any calculation, if they are beyond the classical order. Furthermore, we will only be interested in the conservative dynamics, so we will restrict the components of the momentum transfer q = (q 0 , q) to scale as so that the graviton multiplet mediates instantaneous long-range interactions. Note that the latter expansion involves an additional small parameter, a velocity |v| = q 0 /|q| 1, on top of the classical limit J 1. We will refer to this expansion as the near-static limit, and we delay a more detailed discussion to Section 4.
Finally, in comparing our results to Einstein gravity, it will be useful to take the highenergy or ultra-relativistic limit in which the black holes are highly boosted. This simply makes the hierarchy of scales in Eq. (2.7) more strict In this context, it will be useful to introduce the variable which is simply the relativistic factor of particle 1 in the rest-frame of particle 2 (or vice versa). In terms of this variable the high-energy limit simply corresponds to taking σ 1. Note that in our setup it is important that we take the classical limit first, before taking the high-energy limit, so that J σ. This is equivalent to having the hierarchy of scales in Eq. (2.9). The opposite limit, J σ, is closely connected to the regime of massless high-energy scattering considered in Ref. [24].
In summary, we will be interested in the three limits Generic classical limit:

Integrands from Kaluza-Klein reduction
In this section we construct the tree amplitude and loop integrands for the scattering of the two black holes via Kaluza-Klein (KK) reduction. Ref. [52] studied the case of three-angle misalignment in the BPS charges. While such case is the most rich and interesting, we will focus on the single-angle misalignment case, which is the one we can access via KK reduction from the existing integrands. Let us explain this in more detail: we consider Type IIA supergravity in D = 10 and perform KK reduction on a six-torus of radius R. When dimensionally reducing the massless integrand we will identify the massive black holes with KK gravitons, with ten-dimensional momenta k i and masses arising from the components of momenta outside of D = 4. The supersymmetry algebra in higher dimensions, upon reduction, identifies the extra-dimensional momenta as BPS charges (see e.g. Appendix B of Ref. [52]). There is only one relative angle between the extra dimensional momenta, so the dimensional reduction only provides the result for one-angle misalignment. Because of this, one might perform a rotation to set the momenta along all but two directions to zero and effectively reduce from D = 6. Henceforth, for simplicity, we shall then pretend we are reducing from six dimensions. Then we can write the momenta of the four particles as (3.1) The compactness of the extra dimensions requires the extra dimensional momenta to be discrete and of order ∼ R −1 . We will choose the masses m 1 and m 2 to correspond to the two lightest KK modes, φ 1 , φ 2 . Depending on the momentum in the extra dimensions the massless six-dimensional scalar, φ, will reduce to either of these.
We will see momentarily that the massless integrands for maximal supergravity, have two simplifying features which imply that we just need a few basic rules to perform the KK reduction. First, the loop integrands are proportional to the tree amplitude to all orders. This follows from the supersymmetry Ward identity [54], and implies that the polarization dependence is trivial and factors out of the integrand. Second, through two loops, the integrands are composed only of scalar loop integrals, so we only need to understand how to KK reduce propagators.
Let us first discuss the rules for reducing the massless loop integrand. The integration over loop momentum reduces as where the factors of (2πR) 2 simply relate the D = 6 and D = 4 Newton's constant G = G 6D (2πR) −2 , and the sum is over all possible ways to assign a KK momentum to each leg in a given diagram, subject to the constraint of momentum conservation in the extra dimensions. Since we are choosing our external legs to be two particular KK modes, this means in practice that we should sum over all the ways the external massive particles could route inside the diagram. We are interested in the diagrams that feature the exchange of massless particle in the graviton multiplet, so we will truncate the full massive integrand to this sector. We delay a discussion about the consistency of this truncation to the end of this section. The truncation to massless exchange, together with momentum conservation imposes an additional rule when routing the external particles through the diagram, namely that lines corresponding to different KK modes cannot cross at a three-point vertex.
As an example, consider the massless non-planar double-box integral in Fig. 1(a). It is easy to see that there are two alternative ways to route the mass/extra-dimensional momenta through the diagram, shown in Fig. 1(b) and (c). So this massless integral will yield two contributions to the massive integrand. In contrast, there are also examples in which there there is no way to route the masses. We will find several of these when constructing the two loop integrand.
Finally, using the identifications in Eq. (3.1) we find that the reduction of the external invariants is given by the following simple replacement rules and similarly for loop momenta which follows from the orthogonality of the four-dimensional loop momentum and the extradimensional components of the external momentum.

Tree level amplitude
As a warmup let's start with the tree level amplitude. We will write it as where K is the four point matrix element of the supersymmetric t 8 t 8 R 4 operator (see e.g. Ref. [89], Eq. (9.A.18)). In four dimensions K = [3 4] 4 / 1 2 4 δ (16) (Q), where Q is the on-shell super-momentum [90]. For simplicity we choose the incoming and outgoing states to be complex conjugate scalars φ andφ in the graviton multiplet. The corresponding component of K is simply s 4 and the D-dimensional scalar amplitude is Using our rules for dimensional reduction we find the result Although this is the full amplitude we want to restrict to the massless exchange sector. We can partial fraction (3.7) as which using s − |m 1 + m 2 e iφ | 2 = 2m 1 m 2 (cosh η − cos φ), where η is the relative rapidity, η = arcosh(σ), agrees with Eq. (3.18) of Ref. [52], restricted to the one-angle case.

One-loop integrand
The one-loop massless integrand was constructed long ago in Refs. [53,54] where all the integrals are one-loop boxes with the specified ordering of the external legs. Using the reduction rules described above and KK reduction maps the massless integrals to massive integrals as follows 1423 → I X .
(3.11) Where the integrals are shown in Fig. 2, and 0 indicates that there is no way to route the momenta so the reduction yields zero. Putting all together we find the one-loop integrand

Two-loop integrand
The massless two loop integrand was constructed in Ref. [54] using the unitarity method. It takes the remarkably simple form  where " + cyclic" means adding the two other cyclic permutations of (2, 3, 4) and the integrals, which are all scalar, are shown in Fig. 3. It is easy to find how the integrals map under the dimensional reduction. The planar integrals reduce as follows where I III is the ladder or double-box integral in Fig. 4(a), I H is the H integral in Fig. 5(a), I , I are the self-energy diagrams in Fig. 6(a-b), and the integrals with a bar denote their crossed versions, obtained by exchanging p 2 ↔ −p 3 , which are also shown in the same figures. It is interesting to note that the H and self-energy diagrams come from the dimensional reduction of the same massless diagrams. The non-planar integrals reduce as follows where we will refer to I IX and I XI as non-planar double-boxes, and the rest of the integrals are shown in Figs. 4 and 5. The KK reduced two-loop integrand is then given by     (d) Figure 6: Two loop integrals that include a self interaction.

Comments on the consistency of the integrands
Finally, let us make some brief comments about the consistency of the integrands we have constructed in this section. We have focused on the sector of the theory where KK modes with masses m 1 , m 2 of order R −1 exchange massless particles. This is not a consistent truncation, however, since there is no parametric separation between the masses of the KK modes which are all of order R −1 . 3 This manifests itself in various ways. For instance, the tree amplitude in Eq. 3.7, features the exchange of massive particle of mass ∼ m 1 − m 2 , which is required by crossing symmetry. Similarly, at loop-level, it is known that the sum of the box and crossed-box integrals contains a mass singularity (see e.g. [91]). Consequently, the amplitude in Eq. (3.12) has collinear divergences in violation of the theorem of Ref. [92] which precludes them in quantum gravity. 4 All of these issues are manifestations of the well known fact that there is no consistent quantum theory of a finite number of massive particles coupled to maximal supergravity. In attempting to fix these problems, one is bound to discover the tower of KK modes, which arise from a consistent massless theory in higher dimensions. In spite of these comments, our truncated theory has a well-defined classical Lagrangian and is a useful toy model to explore the questions we are interested in in this paper, so we will ignore all of these issues henceforth.

Integration via velocity differential equations
In the previous section we have constructed the full quantum integrand for the four-point amplitude through two loops. In this section we will calculate the integrals using the method of regions [71,93] to extract the contributions which are relevant for the conservative dynamics. After briefly reviewing the various regions involved the problem, we introduce a new method to calculate the integrals in the potential region, using single-scale fully relativistic differential equations with modified boundary conditions. We illustrate the method using several examples at one and two loops.

Regions and power counting
Following the discussion in Ref. [20], we consider an internal graviton line with fourmomentum = (ω, ), whose components can scale as where we take as reference scale m = m 1 + m 2 , and each scaling defines a region. Note that the four different regions are defined using two small parameters |q| (or J −1 ) and the velocity |v|, which define the classical and non-relativistic limit respectively. Of the four regions, only the potential region contains off-shell modes, which can be integrated out and yield the conservative part of the dynamics. Their contributions can be captured by a non-relativistic EFT which was introduced and put to use in Refs. [16,19,20], and we will utilize in Section 7.
The method of regions [71,93] instructs us to expand the integrand using the scaling corresponding to a given region, and then integrate over the whole space of loop momenta in dimensional regularization. Our goal is to calculate the contributions from the potential region.

Outline of the new method
Ref. [19] introduced a "non-relativistic integration" method by which one must first expand in velocity before expanding in |q|. This produces simple integrals akin to those appearing in NRGR [94,95] at the cost of breaking manifest relativistic invariance in the first step. As explained above the potential region is defined by a double expansion, and we might chose to expand in the opposite order, first in small |q|, and then in velocity. The expansion in small |q| is just the expansion in the soft region where all graviton momentum components are uniformly small (of order |q|). The result of this expansion is a power series in |q| truncated at an appropriate order, with each term in the expansion given by fully relativistic soft integrals with linearized matter propagators. To simplify the expressions, we will apply the well-known method of integration-by-parts reduction [96] to these soft integrals to rewrite them as a linear combination of master integrals. Then we will construct differential equations [73][74][75][76] in the canonical form [77,78] for these master integrals. The choice of a basis of the master integrals will be an important technical point to be discussed later. The selection of pure basis integrals is also facilitated by automated tools [97,98].
The upside of the soft expansion is that it keeps the integrals fully relativistic, but here we are only interested in the contributions from the potential region. Thus, in a second step we should re-expand the integrals in the potential region where graviton momenta are dominated by spatial components, since the potential region isolates conservative classical effects [16,19,20]. After the expansion in the potential region, each term in the previous small-q expansion would be rewritten as a Taylor series in the velocity (ratio of spatial to time components) of external momenta. Unlike the first step, which gives the expansion in small |q| to some finite order, in the second step the velocity expansion can be performed to all orders by using method of differential equations for the soft master integrals. A key observation is that we can construct differential equations for the soft integrals directly before re-expanding in the potential region, as the re-expansion does not change the differential equations, but changes the boundary conditions near the static limit. Thus, it suffices to expand the soft master integrals to leading order in velocity in the potential region, to obtain the boundary conditions that allow us to uniquely solve the differential equations and determine the integrals to all orders in velocity. 5 Let us now explain each of these steps in more detail.

Soft expansion using special variables
In order to carry out the procedure outline above, it will be useful to parametrize the external kinematics as 6 as displayed in Fig. 7. Note that s = (p 1 + p 2 ) 2 = (p 1 +p 2 ) 2 so the physical region is still given by s > (m 1 + m 2 ) 2 . By construction thep i 's are orthogonal to the momentum transfer q = (p 1 + p 4 ), We would like expand the full topologies in the soft region, which in these variables is characterized by the following hierarchy of scales where stands for any combination of graviton momenta ( 1 , 2 , 1 ± 2 , · · · ), or equivalently The massless graviton propagators typically take the form so they have uniform power counting |q| −2 in the small-|q| limit, without further expansion terms. Meanwhile, the momentum of each matter propagator is the sum of an external matter momentump i ± 1 2 q and the momentum injected by gravitons (here is generally 5 The true values of the soft integrals, which will be useful for future calculations beyond conservative classical dynamics, can be obtained by solving differential equations subject to the boundary conditions of the "full" soft integrals near the static limit or another suitable kinematic limit. 6 Notice that in our convention all external p µ i are outgoing, butpi can be either incoming or outgoing. some linear combination of one or more graviton momenta). We have to expand these matter propagators in the soft region, That is all massive propagators are replaced by "eikonal" propagators that are linear in loop momenta. We can further define normalized external momenta, We can then rewrite the denominators of Eq. (4.14) by following Eq. (4.15) and factoring out the scale associated top i from the propagators, where the relevant kinematic factor is again expanded in small |q|. This choice of variables, has the advantage that each order in the expansion is homogeneous in |q|, due to the absence of products between external and graviton momenta in the numerators. In summary, in the soft region the graviton propagators remain unexpanded, while the matter propagators have the form 1/(2u i · ), generally raised to higher powers when we look at terms beyond the leading order in the expansion. Thus, we can write down the following power counting rules applicable at any loop order, before we actually carry out the expansion in the soft region,  At successively higher orders in the expansion Eq. (4.14), we encounter integrals with propagators raised to higher powers as well as higher-degree polynomials in the numerators. Fortunately, all such integrals can be reduced to a finite number of master integrals via integration by parts [96] automated by the Laporta algorithm [99,100], and we use the FIRE6 software package [101] to perform the calculation. This allows the soft expansion result to be expressed in terms of a small number of master integrals, whose values will be calculated by the method of differential equations.

Velocity differential equations for soft integrals
Next we want to integrate the master integrals, which we will do by the method of differential equations. Importantly, by virtue of the normalization (4.15) we have Hence, after the soft expansion, the only dimensionful scale of the integrals is q 2 . The dependence on q 2 of each integral can be easily fixed by dimensional analysis, and the integrals only depend non-trivially on the following dimensionless parameter, (4.20) Hence our differential equations will depend on this single variable, y, which is related to the relativistic Lorentz factor in Eq. , We give this relation to the next-to-leading order in q 2 since it will be used later to convert amplitude results in y to results in σ.
We will construct the differential equations by taking derivatives of the master integrals. The choice of a basis master integrals is not unique; we choose a pure basis in which each master integral has an expansion where each term is a generalized polylogarithms [102][103][104] of uniform transcendentality. This is largely just a technical point, because at the order of expansion needed, the integrals in this paper do not contain any functions more complicated than logarithms (which are a special case of generalized polylogarithms). However, this will yield simple differential equations. A possible form of the differential operator ∂ y , rewritten as derivatives against normalized external momenta u µ i , is The original form d/dy is fine for differentiating the explicit y-dependent factors in the normalization of the master integrals, but the RHS of Eq. (4.22) is needed to differentiate the propagators and numerators expressed in terms of external and internal momenta. After differentiating any of the pure integrals with respect to y, the result can be IBP-reduced back to the basis of master integrals. We will rationalize the square root y 2 − 1 using the change of variable In terms of these variables, the physical region in our scattering processes is given by 1 < y < ∞, i.e. 0 < x < 1.
Our differential equations will take the canonical form [77,78] where A i are numerical matrices and each α i (x), called a symbol letter, is a rational functions in x, and = (4 − D)/2 is the dimensional regularization parameter. 7 The set of the α i is called the symbol alphabet, in the formalism of Ref. [103] which uses "symbols" to elucidate functional identities between generalized polylogarithms. These differential equations can be easily solved, given appropriate boundary conditions. While we could use them to calculate the full soft integrals, we will use them to directly extract the values of the integrals evaluated in the potential region. By expanding in the potential region and summing the expansion to all orders, we have localized the loop integration on the poles of matter propagators. We are essentially dealing with a version of cut integrals (see e.g. Refs. [107][108][109][110][111][112]), which satisfy the same IBP relation and differential equations as original uncut integrals. This is the reason why the only changes are in the boundary conditions, obtained in the near-static limit y → 1 by re-expanding the master integrals in the potential region.

Static boundary conditions from re-expansion in the potential region
We are ready to write down the power counting of momentum components in the potential region, in terms of a small velocity parameter v. Since we have first expanded in the soft region and transitioned to normalized external momenta in Eq. (4.17), we will write down the power counting for u µ i instead of p µ i , and for graviton momenta µ , The factor of |q| is unimportant in our two-step expansion procedure, where the integrals are already homogeneous in q 2 (i.e. proportional to a definite power of q 2 without further corrections) after the soft expansion is carried out. Now we can expand graviton and matter propagators. Recall that graviton propagators ∼ 1/ 2 are unchanged in the soft expansion. Their expansion in the potential region is On the other hand, matter propagators of the form (4.17) are homogeneous in v and the expansion consists of a single term, (4.29) 7 Henn's canonical form can also be used for finite integrals without a dimensional regulator, see [105,106].
The power counting rules for propagators and integration measure in the potential region are We will only need to expand to leading order in |v|, since we only wish to obtain the value of the integrals at one point, to supply a boundary condition. The expanded integrals can be evaluated by residues by performing contour integration over the graviton energies ω. Such energy integrals can be ambiguous until one applies a proper prescription [16,20]. Such a prescription is effectively part of the definition of the potential region which separates it from the larger soft region. Refs. [16,20] presented the prescription in the absence of double poles, i.e. squared matter propagators, but we will show in our examples that when the energy integral prescription is formulated in terms of residues, double poles can be treated in a natural manner and cause no difficulty. As explained in Ref. [20], this prescription generally implies that an integral in the potential region with less than one massive propagators per loop is necessarily zero. Finally, the resulting D − 1-dimensional integrals can be easily evaluated using traditional methods, and provide the desired boundary conditions to solve our soft integrals in the potential region.

One-loop integrals
Next we will illustrate the method above with some simple one-loop examples. We will evaluate all the box-type integrals, which appear in the one-loop N = 8 integrand in Eq. (3.12) with scalar numerator. Adopting the convention of Ref. [113], we remove from our integrals an overall factor of per loop, where µ is the dimensional regularization scale, which is to be restored at the end. In other words, we will write the integration measure for each loop as d D /(iπ D/2 ), where D ≡ 4 − 2 , and multiplying by a factor of Eq. (4.33) per loop in the end to recover results defined with the more common normalization d D /(2π) D .

Box integral
The box integral with two opposite masses has been evaluated in Ref. [91] in dimensional regularization up to order 0 . It has also been discussed in detail in Ref. [20]. As show in Fig. 8, a generic integral in the box topology is of the form (4.34) Where the propagator denominators are explicitlỹ The crossed box integral topologies are related to the box integral by the replacement

Integration-by-parts reduction of soft integrals
Using the soft power counting rules explained in the previous section we see that the box integrals are O(|q| −2 ). Thus, classical power counting requires expanding the integral to subleading powers. The box propagators reduce in the soft expansion to which upon expansion of the integral will generally appear raised to integer powers. The numerators appearing in the expansion are polynomials in ρ i , so each order in the soft expansion is a sum of integrals of the form with each such integral multiplied by a rational function of the external kinematic variables m 2 i , q 2 , and y. As we already mentioned, q 2 is the only dimensionful scale in such integrals. Whenever i 1 or i 2 is non-positive, the integral will become scaleless and vanish in dimensional regularization. 8 Using integration-by-parts reduction, all such integrals are rewritten as linear combinations of the following three master integrals 9 (4.38) 8 Physically speaking, this is because the soft expansion only captures the part of the amplitude that is non-analytic in q 2 and relevant for long-range classical physics. 9 In contrast the full system has 9 master integrals, see e.g. Ref. [20]. (c) f 3 Figure 9: Topologies for the box master integrals. whose corresponding topologies are depicted in Fig. 9. So all integrals given by Eq. (4.37) span not an infinite-dimensional, but a three-dimensional vector space. The above integrals are all proportional to (−q 2 ) − times a q-independent function of the dimensionless parameter y. The basis does not involve the other triangle integral G 0,1,1,1 with a (linearized) matter propagator on the bottom -this is because with the linearized propagator denominators 2u i · the two triangle integrals are identical and we may freely choose either one as part of the basis of master integrals.
Starting from the original box integral Eq. (4.34) with a k = 1, we expand the propagators as in Eqs. (4.14) and (4.17), and perform integration-by-parts reduction to obtain the small-q expansion of the box integral in terms of the three master integrals in Eq. (4.38), where the 1st, 2nd, and 3rd lines are of order |q| −2 , |q| −1 , and |q| 0 , respectively 10 . The bubble integral f 1 will be eventually set to zero because we will evaluate the integrals in the potential region.

Differential equations for soft integrals
Now we will construct differential equations for the three pure master integrals in Eq. (4.38). The original form of the differential operator d/dy is used for differentiating the explicit y-dependent factors in Eq. (4.38), such as y 2 − 1, and the RHS of Eq. (4.22) is used to differentiate the propagators in Eqs. (4.36) and (4.37). After differentiating any of the three pure integrals with respect to y, the result is a sum of integrals of the form Eq. (4.37), and can be IBP-reduced back to the basis Eq. (4.38). After IBP-reduction we use the change of variables from y to x in Eq. (4.23), to rationalize the square roots. The resulting differential equation is where the matrix A is explicitly given by This can be written in the form (4.25) so we recognize x as the only symbol letter for the integrals relevant at one loop.
Static boundary conditions from re-expansion in the potential region Finally, we need to obtain the appropriate boundary conditions to solve the differential equation (4.42) in the potential region. As explained above, we proceed by expanding the pure basis of master integrals Eq. (4.38) in the near-static limit |v| 1, using the rules in Section 4.2.3. After expanding in |v|, each order in the series consists of a sum of integrals of the form These integrals can be evaluated by performing integration over energy ω by residues. We work in a frame where the momentum transfer q µ has no energy component, so the energy of the two graviton lines are ω and −ω, respectively. For convenience, we can further boost our frame so that particle 1 is at rest 11 and u 2 moves in z-direction (4.44) The y variable defined in Eq. (4.20) is related to the above parametrization by v = y 2 − 1.
We symmetrize over the energy components of the two graviton lines, and rewrite Eq. (4.43) using the transformation Then we perform the ω integral by closing the contour either in the upper half plane or the lower half plane, and pick up contributions from poles at finite values of ω, discarding poles at infinity, i.e. neglecting possible non-zero contributions from the arc of a semi-circle contour whose radius tends to infinity. After the ω integral in Eq. (4.43) is carried out in this way, we are left with the spatial integral d D−1 , and the only denominators left are massless quadratic propagators in three dimensions and linear propagators The resulting spatial integrals only depend on a single scale q 2 , and are related to standard propagator integrals.
The bubble integral f 1 in Eq. (4.38) trivially vanishes in the potential region, because there are no poles at finite values of ω and poles at infinity are discarded in our integration prescription. Using the power counting rules in the potential region, Eqs. (4.30) to (4.32), we can see that f 2 and f 3 , i.e. triangle and box integrals with appropriate prefactors that ensure a canonical form of differential equations, both start at O(v 0 ) in the velocity expansion. For example, f 3 has a prefactor y 2 − 1 = v, two matter propagators giving O(1/v 2 ), and an integration measure of O(v), so overall f 3 is of O(v 0 ). This is not surprising, since it is well known that integrals of unit leading singularity can have at most logarithmic singularities in any kinematic limit. To obtain f 2 and f 3 evaluated in the potential region at the leading order in v, we keep only the leading term in Eq. (4.28) for each graviton propagator, and then use Eq. (4.45) to perform the energy integral, leaving spatial integrals . .
Looking forward to the next sections, we will find the solutions to differential equations to be more non-trivial for two-loop integrals.

Result
Substituting the results Eqs.

Crossed box integral
We end with a discussion of the crossed box integrals. As mentioned above, the unexpanded crossed integral is related to the box integral by the crossing replacement u 1 → u 1 , u 2 → −u 2 . Therefore, the same soft differential equations (4.42) are satisfied by these integrals, and one only needs to be careful about the boundary conditions. The specific choice of reference frame Eq. (4.44) is changed by crossing into In terms of Lorentz invariants, this is y → −y. However, our results for the box integral at y > 1 cannot be analytically continued to negative values of y, because the energy integration prescription produces non-analytic behavior in y. For example, when performing the energy integration for f 3 in Eq. (4.38) in the potential region, the two poles lie on the same side of the contour when y < 0, and the contour integration gives zero. The correct result for crossed integrals in the static limit (analogous to Eqs. (4.50)-(4.52) for the box) is Again, the above equations are derived from the static limit but are actually valid to all orders in velocity, because the velocity differential equations have trivial solutions at one loop.

Result
To obtain the small-|q| expansion of the crossed box, we also need to make the y → −y replacement in the coefficients of f i master integrals in Eq. (4.39). The end result for the small-|q| expansion of the crossed box integral is

Two-loop integrals
Next we will evaluate the two-loop integrals needed for the two-loop integrand in Eq. , only the H and H integrals at leading power survive the (q 2 ) 2 suppression of the numerator, and the rest do not contribute in the classical limit. 12 We will describe in detail the computation of the double-box (III) and H-type integrals, and present results for the rest of the integrals. As usual we will strip from our integrals a factor of (4.33) per loop in intermediate steps, to be restored at the end.

Double-box (III)
We first consider generic integrals of the form (4.60) Where the propagators arẽ (4.61) The double-box (III) topology can be embedded in this family of integrals, as shown in Fig. 10. Later we will see that the H topology can also be embedded in the same family. We note that the equal-mass double-box integral has been evaluated in Refs. [114,115] without expansion in the soft or potential region, but the case of generic masses has not been discussed in the literature.

Soft expansion and differential equations
In the soft region, we construct an expansion of the integrand around small | i | ∼ |q|. In the expansion, only the leading order parts ofρ i , denoted by ρ i and given by where negative indices represent numerators rather than denominators. There are a total of 10 master integrals for the III topology 13 as shown in Figs. 11, 12. A pure basis is given by where all the master integrals are normalized to be proportional to (−q 2 ) −2 . The corresponding topologies are depicted in Figs. 11 and 12, where we have separated the integrals which are even and odd in |q|. 13 For reference, in the full equal-mass problem there are 23 master integrals [115].  We perform soft expansion and use IBP-reduction to write the results in terms of the master integrals. The double-box integral is given as the following small-|q| expansion, where the first line is of order 1/q 2 , the second line is of order 1/|q|, and the remaining lines are of order |q| 0 , Since integration-by-parts will only produce analytic coefficients for master integrals, e.g. polynomials in q 2 but not −q 2 , the master integrals f III,1 to f III,7 appear in terms that are even in |q| in the small-|q| expansion of the amplitude, while f III, 8 to f III,10 appear in expansion terms that are odd in |q|. The differential equations for the master integrals are The even-and odd-|q| systems decouple and we can write where the matrices are given by We make a technical observation here. Previously we found that at one loop x, is the only symbol letter. As a consequence only powers of log x will appear in the solutions to the differential equations. In contrast, at two loops, there are multiple symbol letters appearing in the differential equations in Eq. (4.75): {x, 1 ± x}, so the symbol alphabet is larger. This generically results in the solution of the differential equations being (harmonic [116,117]) polylogarithms, but we will see that at leading order in all two-loop integrals only contain logarithms.
Re-expansion in the potential region As described in Section 4, we obtain boundary conditions for the pure basis of soft integrals by re-expanding the integrals in the potential region following Eqs. (4.28) and (4.29), and then integrate over energy components of loop momenta using an appropriate prescription [20]. The energy components of 1 and 2 are written as ω 1 and ω 2 , while the spatial components are written as 1 and 2 .
For the Roman III integral and non-planar variants with exactly three graviton propagators, we follow the prescription of Ref. [20], but with slight modifications to simplify the presentation. First, we symmetrize over 3! permutations of the energy components of the three gravitons, in a way that directly extends the one-loop prescription Eq. (4.45), with the definition ω 3 = −(ω 1 + ω 2 ), and then proceed as usual, i.e. perform the ω 1 and ω 2 contour integrals one by one, closing the contour either above or below the real axis and always neglecting poles at infinity. As an example, we calculate the static limit of G 1,0,1,0,1,1,1,0,0 , which appears in f III,6 in Eq. (4.69) and is shown in Fig. 11(d). In the y = 1 i.e. static limit, the graviton propagators are turned into (D − 1)-dimensional propagators, .
By the prescription Eq. (4.79), this divergent integral is turned into Now let us perform the ω 1 integral by picking up residues in the upper half plane. Only the 4th, 5th, and 6th terms in the bracket of Eq. (4.82) have ω 1 poles in the upper half plane, and in fact the 5th term contributes two poles whose residues add to zero. The result of ω 1 integration is Now we integrate over ω 2 by picking up residues in either the upper or lower half plane, obtaining the same result Putting it back into Eq. (4.80), we obtain (4.85) Now we check that the final result is also independent of the contour choice for ω 1 . If instead we perform the ω 1 integral in Eq. (4.82) by picking up residues in the lower half plane, we obtain a result identical to Eq. (4.83), so the subsequent ω 2 integration also gives the same result as Eq. (4.84). In conclusion, we have verified in this example that once the S 3 symmetrization over graviton energies are performed, the subsequent energy integration has no dependence on contour choice (in the sense of closing above or below the real axis).
Adopting the frame choice Eq. (4.44), and following this prescription, we find that in the static limit, the only non-vanishing master integrals are equal f  Here we just note that all functions have an overall factor of π 2 2 and therefore the transcendental weight of the solutions is effectively reduced by two. Consequently at the order considered, the only polylogarithmic function relevant is log(x) related to the arcsinh function characteristic of 3PM scattering [19,20] by the change of variable Eq. (4.23), (4.87) Going to O( 4 ) we find an additional weight-two function which has no singularity in the entire range 0 < x < 1, so has no singularity in either the static limit y → 1 or the high-energy limit y → ∞. Barring cancellations, it is natural to expect that this function will be relevant at O(G 4 ) (i.e. at the 4PM order). Finally, inserting in Eq.

H and crossed H (H)
Next we will consider the H integral, which can also be embedded in the family of indices in Eqs. (4.61) and (4.62) as shown in Fig. 13. We note that the case of equal masses has been evaluated in [118] without expansion in the soft or potential region. For this topology we only need the leading contribution in |q|, which is even in |q|. Therefore we only give the pure basis of ten master integrals needed to express the even-|q| terms 14 , The corresponding topologies are shown in Fig. 14. In terms of these the soft expansion of the H integral is simply given by The differential equations for these master integrals are where we have only kept the even-|q| sector and the matrices are given by (4.102) We also need to consider the crossed H, or H, integral, in Fig. 5(b), which is just a crossing of the H integral by p 2 ↔ −p 3 . We note, however, that the H and H integrals 14 For reference, in the full equal-mass problem there are 25 master integrals [118].  As mentioned above, we only need to perform the soft expansion of H and H to the leading order, due to the suppression by t 2 = q 4 factor in the numerator, and subleading corrections are not relevant classically. The leading soft expansion of Eq. (4.103) can be obtained from that of the H integral itself by the replacements followed by multiplying the resulting expression by 1/2. Effectively we have "cut" the matter propagators and turned them into delta functions. However, we still need to define how to "cut" matter propagators raised to higher powers, because integrals with squared matter propagators appear when we construct differential equations, and also appear in our choice of a pure basis of master integrals Eqs. (4.90)-(4.99). An appropriate prescription is with the definition HereĜ i 1 , i 2 ,...,i 9 vanishes whenever the integer i 1 or i 2 is non-positive, because the i0 prescription is of no relevance in the numerator, and the terms in one of the square brackets of Eq. (4.107) add to zero. The advantage of this prescription is that it preserves IBP relations and the differential equations Eq. (4.101). In particular, the pure basis of master integrals for the H topology, Eqs. (4.90)-(4.99) can be mapped to the "cut" version f H,n → f cH,n , 1 ≤ n ≤ 10, (4.108) using Eqs. (4.106) and (4.107), and the resulting integrals satisfy differential equations where the matrices, A i,cH , are identical to the ones in Eqs. (4.101) and (4.102) for the differential equations of the original uncut H topology. Hence the solution of the "cut H" differential equations will only differ from the full H in the boundary conditions. In order to obtain the boundary conditions for the "cut H" integrals, we follow a prescription for performing the energy integrals similar to that in the previous section. In this case, the prescription is simply to carry out the ω 1 integral by residues, and then performing the ω 2 integral by residues too. Each of the two integration steps is done by closing the contour either above or below the real axis, picking up residues from poles at finite values and discarding poles at infinity. We find that the only non-vanishing master integrals in the static limit are f (4.111) The part of Eq. (4.111), proportional to log(−q 2 ), which due to the q 4 suppression in the numerator is the only piece relevant for the classical dynamics, agrees with the result in Refs. [19,20].

Non-planar double-box (IX)
Next we discuss the non-planar ladder integral. We only consider the IX topology, noting that the integral XI is identical. The full integral has been discussed in the equal mass case in Ref. [119]. We first consider generic integrals of the form Where the propagators are, as depicted in Fig. 15 The small-|q| expansion consists of integrals of the form where the leading order parts of the propagators are  where the corresponding topologies are shown in Fig. 16 and Fig. 17.   The even-and odd-|q| systems decouple and we can write where the matrices are given by We proceed by computing the boundary condition in the static limit analogously to the planar ladder discussed above. As before the integrals in this limit are evaluated using the residue method, yielding 3d-integrals tabulated in Appendix A. Only

Crossed integrals
In order to evaluate the integrand (3.16) we also need the integrals that that are obtained from III and IX by p 2 → p 3 crossing (denoted III and IX). Since the energy integration step produces non-analytic behavior, these integrals cannot be directly obtained from analytic continuation, and we have to solve the differential equations again. From Eq. (4.23), we can see that x → −x corresponds to the change y → −y, y 2 − 1 → − y 2 − 1. The differential equations for the crossed integrals are thus obtained from the differential equations for the original integrals, Eqs. (4.75) and (4.132), when we change the LHS by where T ∈ III, IX denotes the topology, and change the RHS by The static boundary conditions for III and IX integrals are obtained by the same energy integration method covered before, and are explicitly given by and I (p)

Scattering amplitudes in the potential region
In the previous section we calculated the integrals necessary to evaluate the one-and twoloop conservative amplitudes in the potential region, which we will denote by M 4,(p) . In this section we will put together the integrals to construct such scattering amplitudes.

Tree-level amplitude
For completeness, let us start by considering the tree-level amplitude in Eq. (3.8). In this case the restriction to the potential region is trivial and we simply have which we have written in a form which will be convenient later.

One-loop amplitude
The one-loop integrand for the conservative black-hole amplitude in N = 8 supergravity is given in terms of the sum of the box and crossed box integrals in the potential region, From the results in Eqs. (4.54) and (4.59) we find Note that this formula is valid in arbitrary dimension. In particular it agrees with the soft-integrals in Eqs. (B.36) and (B.40) of Ref. [23]. This reference also calculated the contribution in the potential region at leading order in velocity, which, as expected, did not match the full soft integrals away from the static limit. It is well known that the contributions of soft and potential region coincide at one loop in D = 4, up to differences that are suppressed in the classical limit. Our result shows that this is also true in arbitrary dimensions. As a cross-check we have also calculated the result directly in the soft region, by solving the differential equations for the soft integrals subject to their full boundary conditions without restricting to the potential region, and found agreement to O( 0 ) for both the 1/(−q 2 ) coefficient and the 1/ −q 2 coefficient. Details will be given elsewhere.
With the sum of the boxes at hand we can evaluate the one-loop amplitude (3.12) with the result (5.4)

Two-loop amplitude
Next we use the integrals in Section 4.4 to assemble the two-loop amplitude. The two-loop amplitude in the potential region is given by III + I where the remaining integrals are suppressed in the classical limit. Naively, the ladders and crossed ladders appear with different prefactor in (3.16). We have where "analytic terms" stand for terms with polynomial (including constant) dependence on q 2 , with or without poles in . Such analytic terms give contact terms after Fourier transform to impact parameter space, and are irrelevant for long-range classical physics. Note that the classical log(−q 2 ) arises from the Taylor expansion of (−q 2 ) −2 . With these, together with the H-type integrals in Eq. (4.111), we can evaluate the conservative two-loop amplitude

Eikonal phase, scattering angle and graviton dominance
In this section we will study eikonal exponentiation of the conservative amplitudes directly in momentum space. We will check the exponentiation of the leading and subleading eikonal in the two-loop amplitude. Then we will use the eikonal phase to evaluate the scattering angle in N = 8 supergravity. Finally we will compare the high-energy limit of our result to that of Einstein gravity.

The eikonal phase in N = 8 supergravity
In traditional treatments of eikonal exponentiation, it is customary to Fourier transform the scattering amplitudes to impact parameter space in order to extract the eikonal phase. Here we will take a slightly different approach and study eikonal exponentiation directly in momentum space. There is a simple reason why we prefer this approach: First, in the presence of a Coulomb-like tree-level interaction, such as graviton exchange, the Fourier transform has the side effect of introducing an additional infrared divergence, which in dimensional regularization gives the appearance that one needs to carefully analyze the scattering amplitude at O( ) and keep track of / contributions to extract the eikonal phase at a fixed order. The momentum space approach has the advantage that the Coulomb-like singularities directly cancel, making clear that the O( ) pieces of the L-loop amplitude cannot contribute to the L-loop phase. 16 Working in momentum space comes at a cost nevertheless: simple products in impact parameter space become convolutions in momentum space. However, all convolutions can be easily evaluated as they are equivalent to iterated bubble integrals.
As usual in the eikonal approach, we will consider the amplitude as a function of a D−2dimensional vector, q ⊥ , transverse to the scattering plane, which has the same magnitude as the four-momentum exchange, i.e., q 2 ⊥ = −q 2 (see e.g. Ref. [32]). The conservative amplitude only depends on powers of q, so this poses no problem. The statement of eikonal exponentiation is that one can write the scattering amplitude in the potential region as a convolutional exponential of the eikonal phase where we defined the convolution as integral over the D − 2 dimensional transverse space with a normalization factor N = 4m 1 m 2 √ σ 2 − 1. Equivalently, one can write the inverse relation, We expand δ perturbatively . Then, from the discussion above we can write the phase in terms of the amplitudes δ (0) = M tree 4,(p) , (6.5) Looking at Eqs. (5.1)-(5.8), we see that to calculate the right-hand side of these equations we need the following convolutions , using Eq. (6.11) These expressions respectively cancel the O(|q| −2 ), the O(|q| −1 ) and the O(|q| 0 ) contributions to the two-loop amplitude Eq. (5.8), which arise from the ladder-type diagrams. Therefore, the ladder-type diagrams at two loops give exactly zero contribution to the eikonal exponent, up to the order of q relevant for classical dynamics at O(G 3 ). This cancellation is a check of the exponentiation of the leading and subleading eikonal phase in the two-loop amplitude. Henceforth we will assume exponentiation of the two-loop phase and leave a proof for further work. We note that this zero result relies on delicate cancellations between all six ladder diagrams which leave only the contributions of the H-type diagrams to the two-loop eikonal phase. In summary, putting together Eqs. (5.1)-(5.8) and (6.13)-(6.16) in (6.5)-(6.7) the result of calculation our of the eikonal phase is Note that δ (2) includes an O( 0 |q ⊥ |) which we have not calculated. This however goes beyond the classical power counting and so is a quantum correction to the phase. Finally, we can readily perform the Fourier transform to obtain the more familiar eikonal phase in impact parameter space with the result where we have dropped O( ) and quantum parts. As a cross-check we have verified that the same result is obtained by using the more common approach in which one directly transforms the amplitudes to impact parameter space.

Soft vs. potential and exponentiation
Let us stress that it was very important that we evaluated the amplitude in the potential region to extract the conservative piece. For the one-loop amplitude, the expansion in the soft region differs from that of the potential region at O( |q| 0 ), which, in addition to the suppression, is a quantum correction since the classical dynamics arises from O(1/|q|) terms. For the two-loop amplitude, however, the two expansions still differ from each other at O(|q| 0 ), which is at the same order as the terms responsible for the classical dynamics at two loops, and the difference is also no longer suppressed by . In fact, when we directly evaluate the integrals in the soft region at two loops we find non-exponentiating effects which cause infrared divergences that are not canceled by either matching to non-relativistic EFT or by extracting the eikonal exponent, signaling the appearance of contributions that cannot be interpreted as arising from a conservative potential. 18 The evaluation of the soft integrals using the differential equations above and a detailed discussion of this point will be presented elsewhere.

Scattering angle from eikonal phase
Let us now calculate the gravitational scattering angle from the eikonal phase. The formula relating the two can be derived from the stationary phase approximation of the Fourier transform of the exponentiated impact-parameter amplitude back to momentum space [30], which yields the relation The magnitude of q is related to the scattering angle χ and the magnitude of the threemomentum p in the center of mass by |q| = 2|p| sin χ 2 , (6.25) where in terms of the center of mass energy E = √ s and/or σ From Eqs. (6.24) and (6.25) we can derive the formula for the scattering angle Using this formula we find the following result for the scattering angle 18 This is reminiscent of the situation in the EFT formulation of the Regge limit of massless scattering [120], where contributions from the Glauber region exponentiate while the full soft regions contain nonexponentiating effects.
or separating the different orders Looking ahead, in order to more easily to compare with the results from EFT in the next section, we will write the formula in terms of the angular momentum, J. The angular momentum is defined as J = |b × p| = |b||p| , (6.32) where b is an impact parameter perpendicular the incoming center of mass momentum p. This is however not the impact parameter, b e , which arises naturally from the eikonal phase. Eq. For small angle scattering |b| ∼ |b e |, and the difference is unimportant at leading order. Our results, however, go beyond the leading order and the difference matters. Using the relation (6.34) we find the scattering angle in terms of the angular momentum For later convenience we can rewrite this in terms of the total mass, m, and symmetric mass ratio, ν, as follows Motivated by the universality in the massless case, we will study the high-energy limit of our result by taking σ → ∞ in our result for the scattering angle at order G 3 , which yields This can be compared with the high-energy limit of the Einstein gravity result in Ref. [20] Eq. (11.32) finding perfect agreement. This strongly suggests that the coefficient of the arcsinh term features graviton dominance, and universality also holds in the case of massive scattering. Note that this result does not trivially follow from the massless one since here we impose J σ (or |q| m) so the high energy limit of classical massive scattering is distinct from the Regge limit of massless particle. Admittedly, our calculation provides is only one point of comparison with Einstein gravity, so the question of graviton dominance merits further investigation, either by calculating the scattering angle in other supergravity theories or by directly proving universality. We leave this for future work.

Consistency check from effective field theory
In this section we will calculate the conservative amplitudes using the non-relativistic integration method of Refs. [16,19,20], which is optimized for EFT matching. This method avoids explicit computation of infrared divergent integrals in dimensional regularization, by canceling such integrals between the full theory and the effective field theory using a four-dimensional matching procedure. We will use the EFT Hamiltonian to calculate the scattering angle solving the classical dynamics. Finally, we will compare to our predictions for the amplitude and the angle from the previous section.
The EFT is defined in the center of mass frame where the magnitude of the three-momenta, |p| = |p |, is unchanged in the scattering and the energies are E i = m 2 i + p 2 . In this frame the momentum transfer is purely spatial and given by q = p − p , and the usual Mandelstam invariants are where χ is the scattering angle and we introduced the total center of mass energy, E, and the symmetric energy ratio, ξ, defined as (7.5) We will use these variables throughout this section.

Scattering amplitude with IR subtractions optimized for EFT matching
Here we will use the method of Refs. [19,20] which first expand in the small-velocity limit in the potential region to produce three-dimensional integrals, and then expand in the limit of small q. Divergent integrals will be kept unevaluated, to be canceled against EFT amplitudes in the matching procedure. First let us calculate the scattering amplitudes optimized for EFT matching. At tree level the relevant piece comes from the 1/t pole where we have divided by the non-relativistic normalization 4E 1 E 2 . We will use the notation in Refs. [19,20] and denote the conservative amplitudes in this section with calligraphic M to distinguish them from those evaluated in dimensional regularization in previous sections. The one-loop amplitude can be easily obtained from the one-loop integrand in Eq. (3.12).
As explained in Ref. [20], Sec. 7.2.2 and 7.3.3, the scalar crossed box gives a vanishing contribution in the potential region (in strictly four dimensions), and the box yields the following three dimensional integral Here evanescent terms refer to two classes of terms: (1) terms that are suppressed in or |q| after loop integration, (2) terms that arise from EFT diagrams with insertions of EFT operators suppressed by or |q| omitted from Eq. (7.17). Due to divergences associated with loop integration, terms of class (2) may be naively of the same order of and |q| as terms that directly correspond to four-dimensional classical dynamics, but nevertheless such evanescent terms cancel in the EFT matching procedure and do not contribute to the final results. The one-loop amplitude optimized for EFT matching is then Finally, we extract the two-loop conservative amplitude optimized for EFT matching from the two-loop integrand in Eq. (3.16). Let us first consider the integrals in the first line of such equation. As explained in Ref. [20], when using the non-relativistic integration method all the non-planar scalar ladders vanish. Intuitively this because the energy flow would require the propagation of an antiparticle, which is not allowed, so only the planar double-box contributes in the potential region as [20] We must note that the vanishing of the non-planar integrals is a consequence of the loopby-loop integration procedure used in Ref. [20], which at every stage drops evanescent contributions. In a two-loop integral these can hit at 1/ or 1/|q| pole coming from a different loop and generate finite contributions with classical power-counting such as those calculated in Section 4.4. These contributions arising from evanescent terms are scheme dependent, and, as mentioned above, their ultimate fate is to cancel in the EFT matching procedure. In particular, they will not affect any physical quantity. Thus, as long as the integration in full theory and EFT is done consistently one might drop such evanescent terms. This effectively gives us a four-dimensional regularization method which, in contrast to our eikonal calculation based on dimensional regularization, does not need quantum corrections of O(|q| 0 ), and O( ) contributions at one-loop in order to extract the classical dynamics at two loops.
Next we consider the integrals in the second line of Eq. (3.16). As explained in previous sections only I H and I H contribute with value given by Eq. (4.111), which we reprint here where we dropped 1/ pole terms that do not generate non-analytic dependence on q 2 .
Putting the pieces together we find the full two-loop amplitude optimized for EFT matching For later convenience we rewrite the conservative amplitudes in terms of the total energy, mass and cross ratios as (7.14)

EFT matching and classical Hamiltonian
Following Ref. [16], we want to match the amplitudes above to an EFT with an ordinary Hamiltonian with a potential, which we later will use to solve for the classical dynamics.
The EFT describes two massive scalars interacting with momentum space Lagrangian given by = 4πG q 2 c 1 p 2 + 2π 2 G 2 |q| c 2 p 2 − 2πG 3 log q 2 c 3 p 2 + · · · , (7.17) where for conciseness we have put the external legs on-shell. As in the full theory, here we have also dropped evanescent terms suppressed by or q 2 at each order in G, which can affect the scattering amplitudes but do not have physical effects. If we Fourier transform q back to position space this yields the more familiar potential with an expansion in G/|r|. The EFT amplitudes calculated with the Lagrangian above are very simple. Due to the absence of anti-particles they are given by iterated bubble diagrams. The results up to order G 3 are given by Ref. [20], where c i = c i (p 2 ) and the primes denote derivatives. The EFT matching is performed by requiring M n = M EFT n , which yields the following coefficients for the potential c 1 (p 2 ) = − m 4 ν 2 E 2 ξ 2(σ − cos φ) 2 , (7.19) so the scattering angle calculated from the EFT is χ 1PM = Gm 2 ν J 4(σ − cos φ) 2 √ σ 2 − 1 , (7.28) χ 2PM = 0 , (7.29) which precisely matches our results in Eqs. (6.39)-(6.41) from the eikonal analysis.
One might be tempted to use the Hamiltonian to also calculate the precession of the periastron, ∆Φ, but as explained in Refs. [85,86], there is a simple relation between this quantity and the scattering angle ∆Φ = χ(J) + χ(−J) , (7.31) which implies that odd orders in G (i.e. odd PM orders), which are also odd in J, do not produce a precession, which can be confirmed by explicit calculation using the Hamiltonian. This means that the absence of precession observed in Ref. [52] extends to O(G 3 ), although for trivial reasons, and a calculation at the next order will be needed to test their conjecture of no precession to all orders. The fact that there is a correction to the scattering angle at O(G 3 ), although suppressed in the probe limit, makes us less optimistic about the possibility of the orbits remaining integrable at higher orders.

Conclusion
In this paper we computed the conservative classical dynamics of scattering of two spinless extremal black holes in N = 8 supergravity at O(G 3 ). In Refs. [19,20] the O(G 3 ) (or 3rd-post-Minkowskian) conservative potential in Einstein gravity was calculated using an EFT matching procedure that avoids evaluation of infrared divergent integrals and provides a velocity expansion to high orders. Here, in contrast, we have directly calculated the IRdivergent scattering amplitude in dimensional regularization, and have directly obtained exact velocity dependence using differential equations, without the need to resum a series expansion. This has allowed us to probe the delicate IR structure of eikonal exponentiation, where terms that vanish in four dimensions or vanish in the classical limit have to be evaluated explicitly at one loop, in order to construct IR subtraction terms at two loops to isolate genuine classical contributions at O(G 3 ). Our novel integration method paves the way to a rigorous verification of the velocity resummation of Ref. [19,20], and to streamline further calculations. The ability to evaluate the divergent two-loop amplitudes in dimensional regularization also opens the door to applying the method of Ref. [17,18] which computes classical observables directly from appropriate phase space integrations of the S-matrix. Our differential equations method is highly flexible as the only difference between the soft region and the potential region is in the boundary conditions. The evaluation of the amplitude in the soft region at two loops and the emergence of non-exponentiating terms will be discussed elsewhere.
By computing the classical gravitational scattering angle in both the eikonal approximation and EFT formalism, we have explicitly established their equivalence at O(G 3 ) for the scattering of massive particles for the first time. While the EFT formalism gives a more direct connection to the classical Hamiltonian, the eikonal approximation provides a more direct relation between two gauge-invariant quantities, the scattering amplitude and the scattering angle. It would be interesting to prove the all-order eikonal exponentiation structure for massive scattering from first principles beyond the one-loop case [13], perhaps by generalizing the partially massive case studied at two loops in Ref. [92], and to prove the validity of the eikonal angle formula beyond two loops.
Remarkably, we found that the classical scattering angle of two extremal black holes in N = 8 supergravity coincides in the limit of high energy with that of two Schwarzschild black holes in Einstein gravity [19,20]. Since the classical limit satisfies |q| M and does not commute with the massless limit M → 0, our result is reminiscent of, but not a direct consequence of, the universality of massless gravitational scattering in the Regge limit recently unveiled in Ref. [24], and strongly suggests graviton dominance, whose mechanism still needs to be understood, is generic at this order.
Beyond universality, several aspects of the scattering of black holes in N = 8 supergravity deserve further study. For instance, it would be very interesting to re-analyze the two-loop calculation for dyonic black holes with generic charge misalignments. This might require an improved understanding of the structure of the S-matrix for mutually non-local particles. Furthermore, the integrability conjecture of Ref. [52] should be investigated at the next order, where precession can arise. Given the simplicity of loop integrands in N = 8 supergravity, we expect this highly symmetric theory to be an excellent theoretical laboratory for other aspects of black hole binary dynamics, such as spin-dependent scattering at O(G 3 ) and spinless scattering at O(G 4 ), both of which are unexplored frontiers in post-Minkowskian expansion of black hole binary dynamics, but are amenable to treatment by our techniques. 20 We hope to explore some of these questions in the near future.
In this appendix we present results for dimensionally regularized Feynman integrals in D − 1 = 3 − 2 spatial dimensions, needed for re-expanding the "soft integrals" in the potential region. All of these integrals are the result of evaluating the energy integrals using the residue prescriptions explained in the main text.
Following widely used conventions in the literature on Feynman integrals, the integrals are presented with the following normalization, In the frame chosen the external three-momentum transfer q is in the transverse (x, y) direction, while some integrals have linear propagators of the form 1/ z = 1/( · n z ), where n z is the unit vector in the z-direction. The final results are fully relativistic and functions of q 2 = −q 2 . Unless otherwise shown, we will consider the −i0 prescription to be implicitly present in every propagator.

A.1 One-loop integrals
At one loop we need to evaluate the linearized triangle and bubble integrals in Eqs. (4.49) and (4.48). These can evaluated using traditional methods. Concrete the general linearized triangle integral is given by [113] d D−1 π (D−1)/2 The usual bubble integrals with c = 0 can be recovered by Where the ⊥ i and q ⊥ are the components of i and q in the plane orthogonal to n z , respectively. Now we symmetrize over the two loop momenta, using