Radiative Classical Gravitational Observables at $\mathcal O(G^3)$ from Scattering Amplitudes

We compute classical gravitational observables for the scattering of two spinless black holes in general relativity and $\mathcal N {=} 8$ supergravity in the formalism of Kosower, Maybee, and O'Connell (KMOC). We focus on the gravitational impulse with radiation reaction and the radiated momentum in black hole scattering at $\mathcal O(G^3)$ to all orders in the velocity. These classical observables require the construction and evaluation of certain loop-level quantities which are greatly simplified by harnessing recent advances from scattering amplitudes and collider physics. In particular, we make use of generalized unitarity to construct the relevant loop integrands, employ reverse unitarity, the method of regions, integration-by-parts (IBP), and (canonical) differential equations to simplify and evaluate all loop and phase-space integrals to obtain the classical gravitational observables of interest to two-loop order. The KMOC formalism naturally incorporates radiation effects which enables us to explore these classical quantities beyond the conservative two-body dynamics. From the impulse and the radiated momentum, we extract the scattering angle and the radiated energy. Finally, we discuss universality of the impulse in the high-energy limit and the relation to the eikonal phase.


Introduction
The increasing experimental success of current gravitational wave astronomy [1,2] combined with the design specifications of future detectors [3][4][5] require theoretical predictions for the classical general relativistic two-body problem to keep up with the experimental accuracy [6].
An important tool for the generation of waveform templates, used for detection and parameter estimation, are fast and reliable semi-analytic models of the binary merger. Prominent examples of are the ones provided by the effective one-body (EOB) formalism [7], which take as input numerical simulations combined with analytic results of the late ringdown and early inspiralling phase of the merger. In the latter phase, the two bodies are still widely separated and in a weak-field, slow motion regime-amenable to perturbation theory. Traditionally, this phase is analyzed in the Post-Newtonian (PN) approximation, which is an expansion both in Newton's constant, G, as well as the relative velocity v of the constituents. Both are linked by the virial theorem for a bound system.
Besides the gravitational inspiral, one can also consider scattering (or hyperbolic) events in which compact objects fly past one another while interacting gravitationally. In such events, the deflection of the objects' trajectories is accompanied by the emission of gravitational radiation (dubbed gravitational Bremsstrahlung in analogy to the electromagnetic case). Even though hyperbolic motion events currently appear to be out of experimental reach [8,9], it has been suggested [10] that scattering observables, nonetheless, can be used as input to determine parameters in EOB models which are subsequently applied to the bound state problem. Furthermore, in certain favorable circumstances, it is even possible to directly link bound and unbound observables via analytic continuation [11][12][13]. For scattering kinematics, as opposed to virialized bound state systems, the velocity and G are not necessarily linked expansion parameters and it is therefore possible to explore perturbation theory in G only, to all orders in v, the so-called Post-Minkowskian (PM) expansion. The leading order waveform in this regime has been constructed in papers by Peters, Kovacs, Thorne, and Crowley in the 1970s [14][15][16][17][18], and was recently revisited in Refs. [19][20][21][22][23][24].
In a beautiful paper, Kosower, Maybee, and O'Connell (KMOC) [81] pointed out how to directly relate classical observables to scattering amplitudes using a version of the in-in formalism. In their original work (and in its extension to spinning objects [82]) the formalism was verified at leading order (LO) and next-to-leading order (NLO) by comparing their formulae for the electromagnetic impulse to expressions obtained by solving classical equations of motion. Beyond leading order, however, these checks were performed at the level of unintegrated expressions, and the full evaluation of the corresponding loop and phase-space integrals was left as a problem for the future. Naively, the solution to this problem poses considerable technical challenges. In our previous letter [83], we have emphasized the similarities between the KMOC setup and cross-section calculations in traditional particle-physics settings. We also emphasized how collider physics ideas, like reverse unitarity [84][85][86][87], integration-by-parts reduction [88,89], and (canonical) differential equations [90][91][92][93][94] are ideally suited to render the KMOC formalism a practical computational tool to derive state-of-the-art results for the relativistic two-body dynamics. As an example we computed the radiated momentum in a binary black hole encounter [83].
It is the aim of this work to elaborate on the technical details of our recent letter, and to obtain additional gravitational observables. More concretely, we calculate the O(G 3 ) gravitational impulse, i.e. the momentum change between the initial and final state of one of the scattering black holes. This verifies the classical calculations of Portilla [95,96] and Westpfahl [97,98] done several decades ago, and extends them to one higher order using very different methods. From the impulse, it is possible to extract the radiative scattering angle at O(G 3 ). Inspired by a computation in maximal supergravity [99], the radiative GR angle has been obtained earlier by Damour [100] from a linear response computation and was later confirmed (subject to certain assumptions) by DiVecchia et al. [101], using eikonal methods. This scattering angle, which includes radiation reaction corrections, also cleared up some of the confusions arising in the high-energy limit of the conservative result of Refs. [41,43].
Since the technical bottleneck involves the evaluation of loop integrals, we give a detailed account of the computation of all relevant master integrals that appear in the KMOC setup at O(G 3 ), i.e. at two loops. As will be reviewed in the main text, the KMOC formalism includes both virtual loop amplitudes as well as phase-space integrals over products of lower order amplitudes. Since we are interested in classical physics, these amplitudes are expanded in the → 0 limit, which is equivalent to the soft region in the language of expansion by regions [102]. We therefore lay out how to expand, reduce, and evaluate all relevant two-loop soft integrals via differential equations that were previously adapted to the Post-Minkowskian expansion in classical gravity [79]. One of the key advantages of the KMOC formalism, together with the integration tools described here, is that we can treat inclusive observables such as the total radiated momentum or the impulse, that only depend on a small number of kinematic scales in an efficient and streamlined fashion without having to go through the multi-scale gravitational waveform where subsequent integrations are technically challenging.
The remainder of this work is organized as follows. In section 2, we briefly recall some of the features of the KMOC formalism and summarize how to represent the gravitational impulse (subsec. 2.1) as well as the radiated momentum (subsec. 2.2) in terms of scattering amplitudes and their unitarity cuts. In section 3, we give a broad review of the methods that have already been successfully applied to the conservative two-body dynamics. We discuss the scattering kinematics, the relevant classical regions together with a brief reminder of generalized unitarity that allows us to efficiently derive the relevant loop integrands. Starting from these integrands, we recall the classical expansion (i.e. the soft expansion in the method of regions [102]) in subsection 3.3, before reducing all integrals to a set of independent masters with the help of integration-by-parts relations. In order to obtain analytic expressions for the resulting master integrals, we also review the differential equations of Ref. [79]. Section 4 comments on the novelties of the soft region and introduces reverse unitarity in subsection 4.3 to treat phase-space integrals (involving on-shell delta functions) on the same footing as virtual ones. Section 5 is devoted to the evaluation of the virtual and cut master integrals in the soft region. Evaluation of these integrals require the knowledge of differential equations and their boundary conditions at certain convenient kinematic points. We furthermore discuss the analytic continuation of these solutions to the relevant physical scattering regions. Next, we comment on several general properties of the KMOC observables and simplified properties in terms of scattering amplitudes in section 6. In section 7 we present the final results for the gravitational impulse in N = 8 supergravity and general relativity and we investigate universality in the high-energy limit and the relation to the eikonal calculation in Refs. [101,103]. We close with our conclusions and an outlook to future directions. Appendices A, B, C, and D respectively include details on relevant Fourier transformation identities, the relation between the impulse and the scattering angle, unitarity relations and cutting rules to determine certain phase-space integrals from the imaginary part of virtual diagrams, as well as our conventions for the soft master integrals. The results of all soft master integrals, as well as our conventions are attached to this arXiv submission as computer readable files.
Note added: In the course of this work, we learned from an independent computation of several of the two-loop soft master integrals by Di Vecchia, Heissenberg, Russo, and Veneziano [103] in the context of the eikonal approach to classical gravitational scattering. We are grateful for discussions and comparisons as well as for coordinating publication.

Gravitational observables via the KMOC formalism
To begin, let us briefly review the key features of the KMOC formalism. For a detailed discussion, the reader is referred to the original article [81]. The basic idea of KMOC is to set up a thought experiment for the scattering of two wavepackets evolving from the asymptotic past to the asymptotic future and to measure the change of some observable, ∆O. In the asymptotic past, the wavepackets are represented by |in , an in quantum state constructed from two-particle momentum eigenstates |p 1 , p 2 in with wavefunctions φ i (p i ), which are well separated by an impact parameter b µ |in = dΦ 2 (p 1 , p 2 ) φ 1 (p 1 )φ 2 (p 2 )e i b·p 1 / |p 1 , p 2 in . (2.1) From now on we will drop the in subscript in the momentum eigenstates and leave it implicit. For convenience, we have introduced the Lorentz-invariant multi-particle on-shell phase-space measure 1 dΦ n (p 1 , · · ·, p n ) = We work in an all-outgoing convention for the momenta p i , in which physical incoming (outgoing) states have negative (positive) energy components. The sign in the Heavyside function, Θ, is chosen accordingly. The incoming state |in will evolve to an out state in the asymptotic future, |out , which might contain additional particles produced during the interaction. The change in an observable, O, can be simply obtained by evaluating the difference of the expectation value of the corresponding Hermitean operator, O, between in and out states In quantum mechanics, the out states are related to the in states by the time evolution operator, i.e. the S-matrix: |out = S|in and we can write The stripped kernel I O is related to the matrix elements via where, to arrive at this expression, we have used the unitarity of the S-matrix, S † S = 1.
Following [81], I O can be related to scattering amplitudes by writing S = 1 + iT such that where ∆O is a measurement function acting on the scattering amplitude, M, defined as On the other hand, to evaluate the real contribution, I O, r , we insert a complete set of states which includes a sum over an arbitrary (possibly empty) set of intermediate "messenger" particles, X. Here dΦ 2+|X| (r 1 , r 2 , X) is the multiparticle on-shell phase-space measure of the massive scalars with momenta r 1 and r 2 together with the massless messengers in the set X. The real kernel turns into 10) and the measuring function ∆O only acts on the amplitude on the left of the unitarity cut which was introduced by the intermediate sum over states. Henceforth, we will use the abbreviated graphical notation in the last line of Eq. (2.10), in which the phase space integral with measure d Φ 2+X is understood to be computed over the legs crossing the dashed blue line, and includes the momentum-conserving delta function. We remind the reader that we work in an all-outgoing convention for the particle momenta and note that all momenta crossing the cut flow from the left to the right in Eq. (2.10) and the following. While this formalism can be applied fully quantum mechanically, in this work we are interested in classical observables. This corresponds to the regime where the Compton wavelength of the external particles is the smallest length scale in the problem. Tracking the KMOC argument carefully, this feature implies that most computations boil down to simple plane-wave calculations in the classical limit. Once the dust settles, in the classical limit, the wavepackets sharply peak about their classical values of the momenta which leads to the appearance of on-shell delta functions and one arrives at a compact expression for the classical change of the observable O in terms of the impact parameter b µ , conjugate to the small momentum transfer q µ = p µ The KMOC analysis suggests that we ought to focus on kinematic regions where the massive particle momenta p i are large and scale like O(1) in the classical counting and the fourmomentum transfer q, as well as graviton loop variables that we will denote by i below, scale like O( ). In the effective field theory context, employing terminology from the "method of regions" [102], the classical expansion is therefore equivalent to the so-called soft expansion. 2 We drop an O(q 2 ) quantum piece in the argument of the delta functionsδ(x) and furthermore do not explicitly write the positive energy theta functions Θ(−p 0 1 + q 0 ) Θ(−p 0 2 − q 0 ) which can be set to 1 in the classical limit. The quantum terms originated from parametrizing the final state momenta as p4 = −p1 + q and p3 = −p2 − q. We note that we compute all kernels in terms of fully on-shell objects in terms of specialized variables introduced in section 3.3. From here on out, we work in natural units = 1, but the counting can be restored from the q-expansion.
Furthermore, we will also expand scattering amplitudes in G where the L-loop amplitude is O(G L+1 ). The impulse (kernels) have analogous expansions which we are going to describe in detail in section 6. However, it is already clear from the current statements that loop amplitudes and their unitarity cuts are essential ingredients.

Gravitational Impulse
In this work, we discuss two observables relevant to classical gravitational scattering. The first is the gravitational impulse, ∆p µ i , which is defined as the total change in momentum of one of the particles during the collision. In the KMOC setup this is encoded by the appropriate quantum momentum operator P i , which is measured asymptotically far from the collision region as follows As summarized above, in the classical limit, this is simply a Fourier transform of the impulse kernel I µ p 1 from momentum transfer q to impact parameter space b which is separated into virtual and real contributions, given in terms of the amplitude as where the numerator insertions q µ and 1 arise from the measurement function ∆P µ 1 acting on the respective amplitudes, which extracts the momentum change of the particle 1 line. Note that relative to Eq. (2.10), we have changed variables in the real contribution by shifting the massive intermediate momenta r i = −p i + i , so that all i are small, O( ), in the classical expansion. The impulse on particle 2 can be obtained by simple relabelling.

Radiated momentum
Another observable of interest is the total radiated momentum ∆R µ (which has been first computed in Ref. [83]) carried away in the form of gravitational waves during the scattering of two black holes. This observable is defined by measuring the momentum operator R µ of the emitted messenger particles. As explained in Ref. [81], this observable only receives real contributions and its kernel is given by unitarity cuts

Setup and review
In this section, we review the relevant technology to calculate the scattering amplitudes that serve as building blocks in the KMOC formalism. Most tools were introduced in Refs. [41,43,79], in the context of classical conservative gravitational scattering, i.e. in the potential region (c.f. Eq. (3.7). Readers familiar with Refs. [41,43,79] can skip this review. The novelties of the soft region, which includes dissipative effects, are the subject of Sec. 4.

Kinematics
Before detailing the integrand construction via generalized unitarity, we briefly review the relevant kinematics for the two-to-two scattering of massive black holes in the classical limit. This allows us to link the classical expansion to the soft (small |q|) expansion familiar from the method of regions [102] and further motivate certain truncations in the integrand construction of the classical, conservative sector. The four-particle scattering of massive scalars in an all-outgoing momentum convention is characterized by the following kinematic invariants, We work with a mostly minus metric η µν = diag(1, −1, −1, −1) and the Mandelstam invariants s, t, and u are subject to the usual constraint It will be useful to introduce the total mass and symmetric mass ratio as well the combinations where σ is the relativistic factor of particle 1 in the rest-frame of particle 2 (or vice versa) and h is the total energy-mass ratio. Physical particle scattering in the s-channel corresponds to the region s > (m 1 + m 2 ) 2 (or σ > 1), t = q 2 < 0, and u < 0. In contrast to the massless case, for massive 2 → 2 scattering one can define a Euclidean region where all invariants are negative and the amplitude is real. This region also plays an important role in the evaluation of Feynman integrals in section 5.

The method of regions and the classical limit
We are ultimately interested in classical dynamics and we would like to only retain the minimal amount of information necessary to describe classical black holes. In this limit, the orbital angular momentum of the scattering black hole binary system is much larger than and 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 as an expansion in small |q|. From a heuristic analysis of scales, we perform an expansion in orders of r s /|b|, where the Schwarzschild radius is r s ∼ Gm for some common mass scale m ∼ m 1 + m 2 , and |b| is the relative transverse distance (conjugate to the momentum transfer q) of the system. This implies that the relevant dimensionless expansion parameter is r s /|b| ∼ Gm/|b| ∼ Gm|q|.
For each additional order of Newton's constant G, we need to expand the amplitude up to one additional power of |q|. The relevant term in the q-expansion at tree-, one-loop, and two-loop level are 1/|q| 2 , 1/|q|, and log |q| respectively. The fact that we are interested in the non-analytic terms in the |q| expansion is related to the long-distance effects in impact parameter space (Analytic terms in the |q|-expansion transform to δ-function contributions in b-space). At a given loop order (order in G), terms that are more subleading in |q| are quantum corrections. In summary, at O(G n ), we only need to expand the scattering amplitude of massive particles up to O(|q| n−3 ) in the small-q expansion [104], in order to extract the classical dynamics. In practice, this implies that some loop integrals can be discarded in the amplitude construction, if they are beyond the classical order. Before moving to the actual integrand construction, it is helpful to recall the relevant kinematic scaling of external and loop momenta. We separate temporal and spatial momentum components k = (ω, k) to define the relevant momentum regions [43] hard: ∼ m soft: ∼ |q| . (3.6) The classical limit is equivalent to an expansion in the dimensionless power counting variable λ = |q|/m. For convenience, instead of counting power of λ we can count powers of |q| relative to the scaling dimension of the amplitude. In terms of this scaling, all matter lines of the heavy black holes have hard momenta of p i ∼ m ∼ O(|q| 0 ). As mentioned above, we are interested in classical, long-distance physics mediated by graviton exchange related to the soft region, which further splits into .
In comparison to previous work that primarily focused on conservative dynamics [41,43,79], associated to potential modes, in the KMOC setup, we will directly work in the full soft region. This is owed to the fact that the classical expansion in the KMOC formalism is intimately tied to the q-expansion and we do not have to make further assumptions about small velocities v or instantaneous interactions. However, even in the KMOC setup, we can impose this additional velocity restriction to reproduce conservative results which serve as nontrivial cross-checks of our assembly.

Loop integrands and generalized unitarity
We have seen in section 2 that the extraction of classical gravitational observables within the KMOC formalism involves virtual higher-loop scattering amplitudes as well as their unitarity cuts. Crucially, all ingredients are on shell. This allows to employ a number of simplifying features developed in the context of the scattering amplitudes program over the last several decades. Making use of these novel tools leads to a highly simplified and streamlined construction of the relevant loop-integrands, i.e. the rational functions before the loop or phase-space integrations are performed. For example, the key technology at work is generalized unitarity [25][26][27] and color-kinematics duality [28][29][30][31][32][33] which effectively reduces all gravitational calculations to the computation of sums of products of Yang-Mills tree-level amplitudes. Most details of the loop-integrand construction have already been described elsewhere [41,43], so that we can limit our review to only the most relevant points and remain telegraphic otherwise. Here we give a general review of the conservative sector result, before pointing out certain additions relevant in the full soft region in section 4. It is well known that loop integrands in quantum field theories can be reconstructed from their singularity structure which in turn is entirely dictated by factorization. This idea is formalized in the generalized unitarity framework which allows to construct amplitudes from their unitarity cuts. In the following, we limit our discussion to the case of two-loop contributions, which is of main interest in this work. The generalization to other loop orders is (at least conceptually) straight forward. We are therefore interested in two-loop scattering processes of two massive scalar particles that are minimally coupled to gravity.
To extract classical physics, it is not necessary to obtain the full quantum amplitude, but instead we can directly focus on the relevant pieces that contribute to long-distance interactions between the two black holes mediated by the exchange of gravitons. In the case of conservative dynamics, the gravitons mediate instantaneous interactions causing further simplifications due to the vanishing of certain potential region integrals [43]. In particular, one requires at least one matter line per loop, together with the absence of certain mushroom type graphs (see e.g. Fig. 7) that will become important in the full soft region shortly. Since we are interested in long-distance effects, we can neglect all contact interactions between the two black holes. In the field theory language, this corresponds to setting to zero the four-scalar interactions. From a practical perspective, this limits the number of relevant terms that are interesting for the amplitude construction. For further discussions, see Ref. [43].
In our calculation, we use generalized unitarity in the following way. At two-loops, we write down an ansatz of cubic, local Feynman-like diagrams for all graph topologies that pass our (conservative,) long-distance, non-scaleless, classical power counting constraint. The resulting list of graphs is summarized in Figure 1 [41,43]. For each of the diagrams, we write a yet undetermined numerator. Power-counting dictates that the two-loop numerators have to have mass-dimension 12 as is easily seen by counting derivatives in simple Feynman diagrams. Our numerator ansatz for a given cubic graph Γ is then written in terms of Lorentz dot products between the three independent external momenta p 1 , p 2 , p 3 and the two independent loop momenta 1 , 2 (3.8) Note that this representation of the numerator ansatz secretly includes possible contact terms of the cubic graph, where a contact term denotes any graph Γ with one of the propagators pinched due to a numerator factor. Here, we chose not to write the numerator monomials in such a way that the contact-term stratification is manifest. (For advantages of a contact term basis, see e.g. Refs. [105,106].) In the generalized unitarity setup, it is then the goal to fix the yet undetermined coefficients a i,j by matching the integrand ansatz in terms of cubic graphs against field theory cuts. This is done by taking unitarity cuts of the amplitude integrand, given in terms of products of tree-level amplitudes (on-shell functions), and equating these cuts to the cuts of the cubic diagrams. This procedure yields linear relations for the free coefficients a Γ,j that can be solved in a straightforward manner. Once we match a spanning set of cuts, we are guaranteed to have a correct amplitude integrand. For the conservative classical impulse and radiated momentum at O(G 3 ), such a spanning set of cuts is given in Figure 2. As mentioned before, these unitarity cuts are nothing but products of on-shell tree-level amplitudes summed over the on-shell states. These tree-level amplitudes can in principle be obtained in whichever way one can imagine, even via Feynman diagrams if necessary. For us, we make further use of recent developments in scattering amplitudes, where we can express tree-level gravity amplitudes as the square of Yang-Mills amplitudes via the celebrated BCJ double-copy procedure [28][29][30][31][32][33]. How this is done in practice has already been summarized in the work of Ref. [43] so we will not reiterate these steps here.   Figure 3: Depiction of the parameterization of external momenta useful for the soft expansion.
Ultimately, we write the unitarity-based L-loop amplitude as a sum over cubic L-loop graphs with classical counting where the numerators N Γ ( j , p i ) have been determined by matching unitarity cuts and the denominator is the product of all Feynman propagators P Γ ( j , p i ) of graph Γ. Concretely, for the conservative result at two-loops, this involves the diagrams in Fig. 1 whose numerators have been constraint by matching the unitarity cuts in Fig. 2. We note that all propagators retain the full kinematic dependence and have not yet been expanded in the soft region, i.e. matter propagators are schematically of the form 1/(( + p i ) 2 − m 2 i ).

Soft expansion, integration-by-parts and differential equations
In this section, we briefly review the relevant tools that allow us to go from the integrands derived above, closer to the final integrated result. We first introduce special kinematic variables to facilitate the classical or equivalently soft (small |q|, or small λ) expansion in the context of the method of regions [102]. All relevant definitions have already appeared in Ref. [79], so we are going to be brief. In order to facilitate integration, we perform an integration-by-parts reduction to a minimal set of master integrals that will be solved by differential equations. In the following, we introduce specialized kinematics, depicted in Eq. (3.10), tailored towards the discussion of the soft expansion of the relevant integrals. These variables have the advantage that the new vectors p i are orthogonal to the momentum transfer q, p i · q = 0 , which directly follows from the on-shell conditions p 2 1 = p 2 4 = m 2 1 and p 2 2 = p 2 3 = m 2 2 . Furthermore, s = (p 1 + p 2 ) 2 = (p 1 + p 2 ) 2 so that the physical region s>(m 1 + m 2 ) 2 , q 2 <0 is unaltered. We also define the soft four-velocities of the two black holes u µ i = p µ i /|p i |, such that u 2 i = 1, and where x rationalizes various naturally appearing square-roots later on. Note that these soft velocities coincide with the classical four velocities of the black holes up to irrelevant corrections of O(q) that do not affect the classical observables. As mentioned above, we are interested in the soft expansion, with the following hierarchy of scales | | ∼ |q| |p i |, m, √ s , or equivalently λ 1. Here, schematically represents arbitrary combinations of graviton momenta of the form ( 1 , 2 , 1 ± 2 , . . .) and typical graviton propagators take the form 1 2 , 1 ( −q) 2 , so that they have a homogeneous |q|-scaling and therefore do not require any non-trivial expansion. Note that we can choose a momentum routing so that graviton lines do not involve the individual momenta p i of the external massive particles. On the other hand, matter propagators do have a non-trivial |q| expansion which we express in terms of the dimensionless velocity variables u i such that each order in the expansion is homogeneous in |q| and the mass dependence factorizes. The matter propagators effectively "eikonalize" and the soft expansion to higher orders in |q| can lead to raised propagator powers. Graphically, we denote these eikonalized (or linearized) matter propagators by a double-line notation, see e.g. the diagram in Fig. 7, to distinguish them form unexpanded propagators e.g. in Fig. 1. In order to reduce tensor integrals and integrals with doubled propagators that appear unavoidably in the soft expansion of matter propagators (3.12), we make use of the standard practice in collider physics and use IBP identities [88,89]. These are due to the fact that total derivatives identically vanish in dimensional regularization (see e.g. [107]). By writing sufficiently many total derivatives, one obtains a set of linear relations in the space of Feynman integrals with a given set of propagators. A key insight is that such a space is in general finite dimensional [108] and solving IBPs reduces the task of computing a general integral to the evaluation of a basis of master integrals. The most common approach to solving the system of IBP identities is Laporta's algorithm, [109,110], implemented in a variety of different packages.
In the present work we use FIRE6 [111]. The soft expansion not only implements the classical J limit by truncating at appropriate orders in |q|, but also leads to an enormous simplification of the resulting integrals. Indeed, consistent with the spirit of effective field theory, the appropriate separation of scales allows us to focus on one scale at a time, here |q|, which essentially reduces classical gravitational scattering to a single-scale problem.
The advantage of the new soft variables u µ i and q µ lies in fact that the mass dependence of a general soft (linearized) integral completely factorizes (due to the properties of the expansion of matter propagators in Eq. (3.12)) and the only remaining dimensionful scale is q 2 which can be extracted by simple dimensional analysis. Therefore, in the soft region, the only nontrivial kinematic variable is y = u 1 · u 2 and we are going to find the values of all soft integrals by deriving and solving differential equations in y here and in section 5, respectively. In order to take derivatives with respect to y at the integrand level, we can express the y derivative in terms of the vectors u 1 and u 2 Acting with (3.13) on the set of master integrals g will produce another set Feynman integrals with the same set of propagators, which can subsequently be reduced to the master basis g to yield a differential equation where A(y, ) is a matrix with rational dependence on y and = (4 − D)/2. We can use the freedom in choosing the basis of integrals and the parametrization of the kinematics to simplify the differential equation (3.14). In all cases discussed in this article we are able to choose an appropriate set of canonical master integrals f [93], 4 in terms of which

Full soft integrands and reverse unitarity
As alluded to before, most of the tools describing classical conservative dynamics in terms of scattering amplitudes have been successfully applied before, see e.g. Refs. [41,43,79]. Taking radiation effects into account leads to a few novelties that we are going to discuss in this section. First of all, there are additional contributions to the integrand, which we summarize in subsection 4.1. One additional feature of the KMOC framework is the presence of on-shell phase-space integrals. We employ reverse unitarity [84][85][86][87], well-known from collider physics computations, in subsection 4.3, to efficiently handle such integrals.

New contributions from the soft region
In our general review of the integrand construction in subsection 3.2, we mostly discussed the conservative sector, previously presented in Refs. [41,43]. Here, we are interested in going beyond conservative dynamics and taking radiation effects into account. This requires us to slightly augment the known integrand and include additional terms. The generalized unitarity strategy to determine these additional contributions, though, is the same as for the conservative result. Compared to the conservative dynamics considered in Refs. [41,43], we have additional diagrams depicted on the second line of Figure 5. This is owed to the fact that there are additional terms that survive in the full soft region but are zero in the potential region. In particular, we are interested in classical physics (including radiation contributions) which, as explained in section 2, is encoded in the soft expansion where the momenta of the black holes scale like O(m) and the momentum transfer and the momentum of internal graviton lines scale like O(|q|). This leads to the following summary of |q|-counting rules that facilitate the classical counting (i.e. soft counting) of linearized integrals 5 graviton propagator: ∼ |q| −2 , matter propagator: ∼ |q| −1 , As was pointed out in Ref. [79], graviton propagators scale homogeneously like |q| −2 , whereas matter propagators have a |q|-expansion starting at order |q| −1 . The leading order k in |q|, O(|q| k ) for a given two-loop graph is where the factor 4L comes from the loop measure L i d 4 i ∼ |q| 4L , n v 3 is the number of three-graviton vertices of the graph, n pg the number of graviton propagators, and finally n pm the number of matter propagators. We call two-loop diagrams superclassical (or classically singular), when their leading order term in the |q| expansion starts with k < 0, classical when k = 0, and quantum, if k > 0. With these simple counting-rules, we see that the graph on Figure 4a has k = 8 + 2 × 3 − 2 × 6 − 1 = 1 > 0 which is, as advertised, quantum and we can therefore neglect such contributions.
Another simplification, related to the soft-expansion of Feynman integrals comes from the knowledge that certain integrals become scaleless and therefore integrate to zero. There are simple rules to identify such topologies even before integration which allows us to neglect such terms in the integrand from the outset (see e.g. Figure 4c).
Taking into account the above rules, we find the list of cubic graphs relevant for radiative classical dynamics at O(G 3 ); depicted in Figure 5. Compared to the conservative diagrams shown in Figure 1, the second line is new. These additional diagrams require us to enlarge our spanning set of cuts, compared to the conservative ones depicted in Figure 2 in order to get constraints on the new numerator ansaetze. The spanning set of cuts that allows us to fix the classical integrand in the soft region is given in Figure 6.
We note that some of the unitarity cuts also involve quantum terms and in order to match the full integrand before further classical truncation would require additional cubic graphs not listed in Figure 5. Since these additional terms are purely quantum, we can in principle drop them from our discussion and write the amplitude analogous to Eq. (3.9) where the sum over cubic graphs now contains the additional contributions of the diagrams on the second line of Figure 5 with the numerators fixed by matching the ansatz against the spanning set of cuts in Figure 6. Consistently ignoring such quantum terms in the cut matching procedure, however, is rather subtle.
In the full soft region, completely fixing the unitarity based ansatz requires the matching of rather "deep" unitarity cuts with very few lines put on shell (such as the three-graviton cut on the l.h.s. of Fig. 6). Fully matching these cuts becomes increasingly cumbersome due to the addition of a multitude of quantum terms that are irrelevant for classical physics. One way to circumvent this situation is to step back from the unitarity setup and instead employ simplified gravitational Feynman rules [113,114], closely following the implementations in Ref. [45,115], to target the set of classical Feynman diagrams directly. The simplification of the gravitational Feynman rules is possible due to a judicious choice of gauge-fixing functions. The Feynman diagrams are generated by QGRAF [116], ignoring ghost particles which have no contributions to classical physics. 6 The Lorentz index contractions are carried out with an in-house code to produce numerators in terms of dot products for each diagram. We have checked that the integrand constructed via the simplified Feynman rules matches all relevant classical parts of the unitarity cuts of Fig. 6, so that we are confident in our implementation.
Note that our integrand contains the box-bubble graph on the very right hand side of the bottom-row of Fig. 5 which naively looks like a quantum contribution due to the internal graviton loop. However, by our |q|-counting arguments, this graph is of classical order by virtue of a 1/|q| iteration hitting a O(|q|) quantum contribution. As we will show explicitly in subsection 6.2, these contributions cancel in the classical observables within the KMOC formalism in a way that is similar to eikonal subtractions, see e.g. [99,101].
With the relevant classical virtual two-loop integrand at hand, we also have all relevant terms required for the real contributions in the KMOC setup in Eqs. (2.10) and (2.18). In fact there is no need to construct these cut contributions separately. Instead, we can take our virtual integrand and perform the required unitarity cut. The relevant measurement function in the form of the appropriate loop-momentum insertion for either the impulse or radiated momentum is then simply linked to the labeling of the cut. As such, we have now constructed all integrands that make an appearance in the KMOC formalism and we can turn our attention to the novelties of integrating in the full soft region as well as tools that handle these cut, or phase-space integrals. This is what we turn to next.

Soft expansion and partial fractioning
There is one aspect of special soft integrals that warrants discussion. A particular feature of mushroom-type integrals (see e.g. Fig. 7) pertains to the fact that once the matter propagators are linearized via the soft expansion of Eq. (3.12), some of the propagators in a given diagram might become linearly dependent. However, this can be addressed by partial fractioning. In the example given by Figure 7, the three matter propagators on the top are These expressions are linearly dependent and partial fractioning allows us to split diagrams with all three propagators into terms with at most two of the propagators at a time 1 Pictorially, this identity is expressed as a relation between mushroom-type integrals Figure 7: Linear mushroom integral.
where the dot represents a doubled propagator. Due to the soft expansion to higher orders, we also need to treat raised propagator powers with partial fractioning analogously to Eq. (4.3). The propagators on the right-hand side of (4.4) do not satisfy any linear relations and they can be embedded into different top-level families where the four-point vertex is blown up. For example Similarly, we can proceed for all other mushroom topologies of Figure 5. In the soft region, these additional diagram topologies were required to match the classical part of the amplitude. However, upon soft expansion, their propagator structures overlap with existing topologies so that no new integral families are required. The first integral on the right-hand-side of Eq. (4.4) actually vanishes in dimensional regularization, because it factorize into a box integral times a matter self energy diagram which is scaleless in the soft region. More generally, non-factorizing integrals, where the loop momenta can be routed such that the integral is independent of the momentum transfer q are zero in dimensional regularization. One such example was presented in Figure 4c. (4.6)

Reverse unitarity
We have seen in section 2, that classical gravitational observables, like the impulse kernel (2.17), or the radiated momentum kernel (2.18) involve not only virtual amplitudes, but also certain unitarity cuts. These cut contributions include an integral over the onshell phase space of the exchanged states. In order to efficiently evaluate such phase-space integrals, we follow our earlier letter [83], where we took inspiration from the enormous progress in cross-section calculations and the computation of collider physics observables where similar real contributions appear. For some time, it has proven advantageous to handled phase-space integrals on the same footing as virtual integrals. This idea has formally been implemented via reverse unitarity [84][85][86][87], where one replaces on-shell delta functions and their n-th derivatives by the difference of (appropriate powers of) propagators with varying iε prescription Trading all delta functions by differences of propagators allows us to employ standard tools for loop integrals such as dimensional regularization, IBP reduction [88,89], and (canonical) differential equations [90][91][92][93][94] to evaluate a minimal set of master integrals. From a practical perspective, we can treat any on-shell delta function as a regular propagator. This is owed to the fact that integration-by-parts identities, crucial in the derivation of the differential equations, are insensitive to the Feynman iε. The same is true for the partial fractioning discussed in subsection 4.2. This significantly simplifies our computations and circumvents the difficulties in having to evaluate integrals containing derivatives of delta functions that would otherwise appear.
As will be explained in section 5, the differences between the cut integrals compared to the virtual ones arise from the following facts: • Certain cuts can break diagram symmetries of the virtual diagram.
• Various terms in the differential equations can be omitted because the appropriate master integrals do not have the desired unitarity cut.
• The boundary conditions for the differential equations for the cut master integrals change, relative to the virtual integrals.
Of course, all the described properties of cut integrals are well known from collider physics applications and we adapt them to the gravitational setting here. To reiterate, the huge advantage of the reverse unitarity setup arises from the fact that we can directly treat the phase-space integrations for the inclusive classical observables (such as the gravitational impulse or the radiated momentum) in one go without having to perform sequential integrations over the gravitational waveform. Even for more exclusive observables, like the radiated energy spectrum, we can add one new variable at a time, which still leads to simplified integrals. A detailed discussion of such quantities is left to the future.

Evaluation of soft master integrals
In this section, we explain the evaluation of the soft master integrals relevant for the computation of radiative observables up to O(G 3 ). We first review the relevant kinematic domain and show how to compute the soft integrals at one-loop level as a warm-up exercise. Subsequently, we explain how to evaluate virtual integrals as well as two and three particle cuts at two-loop level. This completes the set of relevant master integrals for classical observables at O(G 3 ). Ultimately, our set of master integrals can be recycled to obtain analogous classical observables for different theories (such as quantum electrodynamics), or for spinning black holes once the relevant integrands are available.

Soft one-loop integrals: Euclidean region and analytic continuation
Before elaborating on the evaluation of soft master integrals via differential equations, it is illustrative to recall the kinematic dependence of soft integrals. As mentioned before, upon soft expansion, all soft master integrals have their external mass and momentum transfer (−q 2 ) dependence determined by simple dimensional analysis. The only nontrivial Figure 8: Kinematics in x and y space. The x plane is a double cover of the y plane, i.e. the two points x and 1/x map to the same point y, we therefore focus on points inside the unit disc |x| ≤ 1. Physical s-channel scattering (green region) corresponds to y > 1 and Im(y) = 0 + , i.e. 0 < x < 1 and Im(x) = 0 − . Physical u-channel scattering (red region) corresponds to y < −1 and Im(y) = 0 − , i.e. −1 < x < 0 and Im(x) = 0 + . There is also a Euclidean region connecting the two for real −1 < y < +1 which corresponds to points on the unit circle |x| = 1. Due to the double-cover property, x = +i and x = −i map to y = 0.
kinematic dependence of the integral is through the dimensionless variable y = u 1 ·u 2 (or equivalently in terms of x, s.t. y= 1+x 2 2x ). The =(4−D)/2 dependence can sometimes be computed exactly, otherwise we work in an expansion around =0.
The soft expansion, reviewed in section 3.3, is defined in a manifestly relativistic way and therefore soft integrals are genuine D = 4−2 dimensional Feynman integrals, although involving linearized propagators. The manifest covariance has the benefit that the integrals are analytic functions of y and we can use analytic continuation to relate integrals in different kinematic regions. In contrast to the massless case, for massive 2 → 2 scattering one can define a Euclidean region where all Lorentz invariants are below production threshold and the amplitude is real. As will become apparent from the explicit examples below, it is advantageous to compute the soft master integrals in the Euclidean region and then analytically continue to the desired scattering kinematics. These regions, together with the analytic continuation are summarized in Figure 8. As an illustrative example, consider the one-loop box family of the form (see also Ref. [79] In the following, we adopt the normalization conventions of Ref. [117] and remove an overall factor of where the linearized propagators are explicitly There are 3 master integrals The differential equation is The system only has a single letter and can thus be integrated to all orders, using x = −1 as a boundary condition For the boundary conditions, we can directly evaluate the bubble and triangle integrals, using the master formula for the linearized triangle with arbitrary powers of the propagators (see e.g. Ref. [117]) Applied to the integrals of interest, we find For the boundary condition of the box integral, we resort to the method of regions. For this analysis it is convenient to choose a frame which coincides with the rest frame of u 1 up to q-corrections In this frame we have y = √ v 2 + 1. By crossing the limit x → −1 corresponds to the static limit of the crossed box integral. The leading contribution in this region comes from the potential region, where the integral vanishes, the subleading contribution from the "quantum soft" region scales as O(v), so we find f 3 (−1) = 0. With this we can evaluate the crossed box integral log(x) By analytic continuation, we obtain the box integral Now as we discussed before, we can use the same differential equation (5.4) for the twoparticle cuts. In this case the triangle and bubble functions are trivially zero because they do not have a relevant cut. This directly implies that the cut box integral is constant and given by its value at y = 1. The integral can be directly evaluated and reduced to a (D − 2)-dimensional Euclidean bubble integral = 1 The triangle integral has no x dependence and is therefore completely specified by the boundary conditions at x = −1.

Virtual two-loop integrals
The most involved part in our KMOC computation of the classical gravitational observables at O(G 3 ) is the evaluation of the virtual two-loop soft integrals. The differential equation matrices have been constructed in Ref. [79], with the exception of the odd-in-|q| integrals for the H family which we add in this work. 9 A complete list of all master integrals, together with our conventions, is given in Appendix D.
Similarly to the one-loop discussion, it is advantageous to first evaluate all integrals in the Euclidean region and then analytically continue to the desired scattering kinematics. In the Euclidean region, the integrals are real-valued which serves as a valuable cross check on the calculation and also facilitates numerical verification against e.g. PySecDec [118,119]. 10 In the most general (nonplanar) case, which the IX integral in Figure 11 is an example of, the Euclidean region is −1 < y < 1. The scattering regions are (1) s-channel: y > 1, Im(y) = 0 + , (2) u-channel: y < −1, Im(y) = 0 − . For planar integrals, the Euclidean region is larger, given by y < 1, and includes the u-channel scattering region, so a nontrivial analytic continuation is needed only for the s-channel scattering region.
The boundary conditions can be fixed by various methods, we discuss a set of sufficient conditions given by known single scale integrals, regularity and the method of regions. 8 The additional factor of i is due to our conventions of the integral measure. 9 Since IBP reduction of an integral gives a sum of master integral with analytic-in-q 2 coefficients, integrals that scale like odd powers of |q| form a decoupled system under IBP relations and differential equations. 10 Some pure integrals have a prefactor y 2 − 1, and become purely imaginary without a real part.
(a) Figure 9: Single-scale integrals appearing as boundary values for two-loop soft integrals.
Single scale integrals First, there are a handful of single-scale integrals independent of y, for example the sunrise integral in Figure 9b. These integrals are either factorizing into one-loop integrals (Figure 9e-9g), or can be be performed loop-by-loop (Figure 9b-9d), eventually reducing them to one-loop integrals which can be evaluated using the master formula (5.6). Some integrals, like the double triangle in Figure 9a, can be computed using the trick of symmetrizing over the graviton momenta. This turns the matter propagators into delta functions and the integral becomes three-dimensional (cf. Appendix A of Ref. [79]).
Regularity Another input is the regularity of integrals in the Euclidean region. In our kinematic parametrization, this translates to the statement that the s-channel planar integrals have to be regular at y = −1 or x = −1. At two loops, this only provides nontrivial constraints for planar integrals in the u-channel. Finally, integrals odd under parity y 2 − 1 → − y 2 − 1 have to vanish at y = −1. Using all these conditions, we obtain all-order boundary conditions for all integrals in the H family.
Analysis of regions For the remaining integrals we obtain boundary conditions by the method of regions, splitting the soft region into subregions defined in Eq. (3.7). We again adapt the frame defined in Eq. (5.9). As an illustrative example of how we can use the method of regions to obtain boundary values at two-loops, we consider the scalar III integral. By naive velocity power-counting the leading contribution in the small velocity limit comes from the region where all gravitons are in the potential region. In this region the integral scales as 1/v 2 . All other regions are suppressed in velocity -the next-toleading contribution arises from the region where one of the gravitons is in the radiation region and scales as v −1−2 by naive power-counting. Therefore the value of the canonically normalized scalar III integral at v = 0 equals the potential-region boundary value, which has been computed in Ref. [79], Combining these different methods allows us to determine a complete set of boundary conditions to all orders in for the virtual soft integrals relevant at O(G 3 ). The complete velocity-dependent analytic values of the soft master integrals can then be obtained by integrating the differential equations from Ref. [79] in conjunction with the soft boundary conditions. The explicit values of all master integrals are given in the ancillary files accompanying this work.

Two-particle cut integrals from sub-loop integration
A general two-particle cut integral can be evaluated by sub-loop integration. This is true for any loop order, however we focus on the two-loop case here. This can be seen as follows: all two-particle cut integrals that we encounter in the KMOC-formalism, including those with a numerator insertion, can be cast into the form where I L and I R are the sub loop integrals to the left and right of the two-particle cut, respectively. Since the momentum transfer is the only dimensionful quantity, it can be fixed by dimensional analysis This leads to the following formula relating the sub-loop integrations to the cut The δ-functions localize the integral to an (euclidean) integral over transverse space 11 As a concrete example we consider the two-particle-cut of the scalar III, where we have log x+iπ Likewise, we can evaluate all other two-particle cut integrals by our loop-by-loop integration technique, using the known one-loop building blocks and the general master formula in Eq. (5.18). Notably, the imaginary parts of the one-loop building blocks change, depending on whether they are inserted to the left or to the right of the cut legs.

Triple-cut integrals from differential equations and Cutkosky rules
We first start by considering triple cuts of integrals inside the H family. It will turn out that all other triple-cut integrals can be obtained from these via differential equations. Using the cutting rules, reviewed in appendix C, we can relate the triple cut to the imaginary part for a H-type integral we find Since we already computed the virtual integrals in subsection 5.2, this allows to obtain all triple-cuts in the H family. To give a concrete example, we have We notice that there are only four integrals which cannot be embedded into the H topology or its crossing, namely the III and the IX and the planar and non-planar box triangles. For these integrals we can make use of a particular feature of the differential equation, which is that derivatives of these integrals are expressible in terms of integrals inside the H topology. Therefore these integrals can be computed by direct integration. As a concrete example we consider a triple cut of the IX integral. The derivative of the canonically normalized IX integral is proportional to a N-type integral, which is known from the H family For the boundary condition we make use of the method of regions. The power-counting for a cut integral is identical to the corresponding virtual integral. For the triple-cut contribution, one of the gravitons is on-shell and therefore this integral receives no contribution from the potential region. The leading behavior as v → 0 is therefore dictated by the potentialradiation region which scales as O(v −1−2 ) by the power-counting of Eq. (3.6). As the integral appearing in the canonical differential equation is normalized by a factor of y 2 −1 = v 2 (see appendix D for details) it vanishes in the static limit and using Eq. (5.22) we find Similar ideas also apply to the three-particle cut integrals of III and will not be displayed explicitly. This concludes our discussion of all relevant virtual and cut master integrals required for the determination of O(G 3 ) classical observables in the KMOC formalism.

Integral checks
We have performed several consistency checks on our virtual and phase-space integrals. As mentioned previously, for the virtual two-loop integrals, we have performed extensive numerical checks to high precision against PySecDec [118,119] in the Euclidean region where the integrals are real-valued. As a cross-check of our analytic continuation, we have furthermore performed numerical comparisons for scattering kinematics in s and u-channel regions, where our analytic results agree with numerical values within numerical errors.
We also compared our results to available analytic expressions for the full integrals in the equal mass m 1 =m 2 case, finding agreement for the non-analytic-in-q parts for the Htype integrals [120] and ladder-integrals [121][122][123], to the orders of = (4 − D)/2 available in the literature.
Furthermore, we have checked the results of our master integrals against cutting rules. For example, we have checked that Eqs. (C. 16) and (C.17) hold as relations for the softexpanded master integrals, i.e. with all quadratic matter propagators replaced by their linearized expressions at the leading order expansion in |q| given by the first term on the r.h.s. of Eq. (3.12).
Lastly, our virtual master integrals have been checked against an independent calculation [103] which we learned from private communications. Each integral has been checked to the maximum order of that has been computed in both papers.
We note that all velocity-dependent functions satisfy a first-entry condition [124], where only x is allowed as first symbol [125][126][127] entry. This is obvious for the one-loop integrals (which only contain log x), but becomes nontrivial at two-loop order and suggests potentially further simplifications by eliminating more explicit boundary value evaluations due to this analyticity property.

Simplifications in the KMOC setup
We have reviewed the KMOC formalism in section 2, together with general formulae for the gravitational impulse kernel (2.17), and the radiated momentum kernel (2.18). Here, we would like to discuss a convenient organization of these quantities as well as aspects of their perturbative expansions, before presenting their explicit results in maximal supergravity and general relativity up to O(G 3 ) in section 7. More concretely, we use unitarity and the cutting rules to obtain simplified KMOC formulae where certain properties of the impulse kernel, such as its reality properties or the absence of superclassical term are more manifest.
To begin the discussion, it is convenient to decompose the total impulse into its transverse, ∆p ⊥ , and longitudinal, ∆p u , components such that u i ·∆p ⊥ =0 and q·∆p u =0. (For the relevant kinematic definitions, c.f. the beginning of section 3.3.) Correspondingly, the impulse kernel can be written as where we have defined dual four-velocitieš which satisfy u i ·ǔ j = δ ij and are still orthogonal to the momentum transfer q. Decomposing the loop momentum dependent impulse numerator in a similar fashion 12 reveals that only the transverse part of the impulse has a virtual contribution whereas the longitudinal part is purely real, i.e. it only receives contributions from the unitarity cut terms Note that due to the difference in the |a| scaling of q µ andǔ µ i , we have to expand the longitudinal impulse kernels to one higher order in |q| compared to the transverse ones. Finally, all classical observables are real 13 (not complex), and the various factors of i in the KMOC setup serve this purpose. In particular, the transverse KMOC kernels need to be purely real to yield a real result after the final Fourier transform (Eq. (2.11)), whereas the longitudinal kernels are purely imaginary. Indeed, it will serve as a nontrivial check of our computation, that all imaginary contributions to the classical observables cancel.

Leading and next-to-leading order impulse
At leading order, O(G), the impulse kernel is given by the tree level scattering amplitude (6.7) 12 In principle, there is an orthogonal direction ε(·, q,ǔ1,ǔ2) that, however, does not play a role. 13 The reality properties of (possibly) complex quantities should not to be confused with our nomenclature of real, i.e. cut, contributions to various classical observables.
There is only the virtual contribution at this order, since the scattering amplitude starts at O(G) and the real contribution is quadratic in the amplitude. At next-to-leading order, O(G 2 ), the impulse kernel receives both virtual and real contributions. The transverse component is 8) and the longitudinal component is In Ref. [81] it was shown that the superclassical part of the one-loop virtual amplitude cancels at the integrand level when expanded in the classical limit. Here, we offer an alternative argument that will streamline the calculation of the kernel. The basic observation is that the cut has a horizontal flip symmetry which does not change the sign of the integral. Thus one might average over the two different labellings of loop momenta, related by 1 ↔ q− 1 1 2 This means that the transverse impulse numerator insertion is in fact constant and the cut contribution can be related to the imaginary part of the amplitude via the unitarity relation in Eq. (C.5). Thus, by virtue of Eqs. (6.10) and (C.5), all imaginary parts cancel and we find that the transverse classical impulse kernel is given by the real part of the one-loop amplitude, We will explicitly show in section 7.2, that at this order the superclassical pieces are contained in the imaginary part of the amplitude. Roughly, this can be understood from the fact that the imaginary part is related by unitarity (C.5) to a cut that corresponds to the iteration of lower orders and therefore cancel in the impulse kernel. The computation of the longitudinal kernels I u i can also be simplified. We begin by noting that the four velocities u i in the numerator insertions in Eq. (6.9) satisfy u i ·q=0. Thus, we can express them in terms of the momenta of the scattering amplitude (see the definition of the soft variables in section 3.1) as u 1 =(−p 1 +q/2)/m 1 + O(q 2 ) and u 2 =(−p 2 −q/2)/m 2 + O(q 2 ), and we can write (6.12) We have used the on-shell conditions ( 1 −p 1 ) 2 −m 2 1 = 0, ( 2 −p 2 ) 2 −m 2 2 =( 1 +p 2 ) 2 −m 2 2 =0 and the · · · indicates contributions that pinch graviton propagators as well as quantum suppressed terms, which can be ignored. At one-loop, the former yield short distance matter contact terms which are irrelevant for widely separated black holes. The resulting numerator is again loop-momentum independent. Using unitarity (C.5), we find that the longitudinal impulse kernels are directly proportional to the imaginary part of the one-loop amplitude. Collecting all the ingredients we obtain a more direct relation between the impulse kernel and the scattering amplitude at this order from which we learn that the impulse can be directly extracted from the virtual amplitude without the need of evaluating any phase space integrals.

Next-to-next-to-leading order impulse
Next we discuss the simplifications at next-to-next-to-leading order, O(G 3 ). We will focus on the transverse part of the impulse and offer some comments about the longitudinal part.

Transverse part
At this order the transverse impulse kernel is given by 14) The three-particle cut in the second line always produces a real contribution. This cut enjoys the same horizontal flip symmetry as the one-loop two-particle cut. Thus one might once again average over the two different labellings of loop momenta as in Eq. (6.10) The same considerations are valid for the real part of two-particle cuts. Any contribution and its horizontally flipped version combine to give a trivial impulse numerator upon averaging as in Eq. (6.10) such that where we have dropped the restriction to the real part in the second line since the imaginary parts cancel in the sum in the absence of a nontrivial numerator. Therefore, the three-, and (real part of the) two-particle cut real contributions combine to cancel the imaginary part of the virtual amplitude in the impulse kernel, by virtue of unitarity (C.6). The remaining pieces are the real part of the virtual amplitude together with the imaginary part of the two-particle cuts. The latter arises from the imaginary parts of the one-loop amplitudes on either side of the cut, which have opposite sign. Once again, we can simplify these by using one-loop unitarity (C.5) on the one-loop amplitude on the left of the cut where the additional phase-space integration on the r.h.s is over the newly cut legs, denoted by 3 −p 1 and 4 −p 2 . The contribution with the one-loop amplitude on the right of the cut has the opposite sign. Combining both into a single term by choosing uniform labels we find our final expression for the transverse impulse Eq. (6.18) provides a simplified prescription for the calculation of the impulse kernel, which reveals among other things, that the purpose of the (real part of) two-particle and threeparticle cuts at this order is to cancel part of the imaginary part of the virtual amplitude. The only non-trivial contribution of the cuts (i.e. the real part of the impulse kernel), take the form of the cubed of the tree amplitude, which is reminiscent of the calculation of the eikonal phase. Eq. (6.18) also manifests the fact that the transverse impulse kernel is real, as required by the fact that its Fourier transform is the transverse impulse which is a physical observable and hence also real. It also facilitates exposing the cancellation of certain quantum contributions to the transverse impulse. For example, diagrams involving self-energy corrections to virtual graviton propagators are quantum corrections (i.e. not relevant for classical physics). Individual diagrams of this class can still be of order q 0 , i.e. naively of classical order, according to the power counting rules in Eqs. we learn that such contributions are imaginary and indeed absent in Eq. (6.18). We note that this integral only receives contributions from the "quantum soft" velocity region (see Eq. (3.7)). It would be interesting to explore whether the contributions from such region can be consistently dropped at an earlier stage in the calculation without spoiling some of the advantages of the full soft region computations.
In addition, Eq. (6.18) allows us to show the integrand level cancellation of superclassical terms. We first notice that all tree-amplitudes entering the iterated two-particle cut are on the same footing and are functions M(σ, q 2 i ) of the respective momentum transfer. Only the on-shell delta functions break the invariance of the cut under the permutation of the q i 's. However at leading order in the classical expansion we have and likewise for the delta functions involving p 2 . Therefore, realizing that 1 − 3 = − q 2 and summing over the cyclic relabellings, we find The leading superclassical part of the virtual amplitude is purely contained in the planar double-box diagram. The same symmetrization relation can be applied to this diagram. Carefully keeping track of the i in the denominators yields which generates on-shell delta functions from the four matter propagators and cancels against (6.21). Naively, there still remains a superclassical contribution in the virtual amplitude at O(|q| −1 ), but this can be shown to be purely imaginary so it is manifestly absent in Eq. (6.18). An alternative, and perhaps more explicit, derivation of our simplified formulae proceeds by considering the contribution from each diagram in our integrand to Eq. (3.9), together with the cutting rules in appendix C. We now give an explicit example of such proof focusing on the contribution of the III diagram to the impulse on the bottom matter line. The impulse numerator combines into a scalar in the sum over the two two-particle cuts as follows, where we used 1 + 2 + 3 = q in the last line. In the second term of the second line we used the fact that the horizontal flip symmetry of the diagram is unaffected when considering only the real part of the involved diagrams, despite the fact that the r.h.s of the cut in each diagram represents a complex conjugated amplitude. Similarly, the impulse numerator combines into a scalar in the sum over the two three-particle cuts. Here we do not need to take the real part because all these expressions are real by themselves, with tree-level expressions on both l.h.s and r.h.s of the cuts.
where we have again used the horizontal flip-symmetry of the cut diagram between the second term on the first and second line. Having canceled the nontrivial impulse numerators for both the three-particle cut and two-particle cut contributions, it can be seen that the imaginary parts cancel between all contributions to Eq. (2.17) originating from the III diagram, using the three-term relation from Cutkosky rules given by Eq. (C.16).

Longitudinal part
Finally, we discuss briefly the simplification in the computation of the longitudinal part of the impulse. Recall that neither part of the longitudinal impulse kernel receives virtual contributions, I u i ,v = 0 and at this order the real contributions are The subsequent Fourier transform does not flip reality properties so that we find the following is true for the longitudinal part of the impulse: Three-particle cuts give real contributions to the impulse kernel. Every diagram and its horizontally flipped version contribute equally, due to the identity The same considerations are valid for the real part of two-particle-cut contributions, so every diagram and its horizontally flipped version contribute equally. For the imaginary parts of double-cut contributions, an extra sign difference causes cancellation between each diagram and its horizontally flipped version. Therefore, to calculate the longitudinal impulse, we only need three-particle-cut contributions and the real part of two-particle-cut contributions. Since the double-cut integrand contains an overall factor of i, and only odd-in-|q| master integrals are needed for the longitudinal impulse, we only need the box-triangle master integral Eq. (D.14) with a cut on the box side as well as its horizontally flipped version.

Results
In this section, we present the results of our computation of the two classical gravitational observables studied in this work: the impulse and the radiated momentum for the scattering of two black holes both in N = 8 supergravity and in general relativity through O(G 3 ).
For the maximally supersymmetric case, the study of this scattering process was initiated in Ref. [128], and revisited in Ref. [79] in the conservative sector. The loop integrands up to two-loops (or next-to-next-to-leading order) were constructed in Ref. [79] by dimensionally reducing the known massless loop integrands [129,130], and are reproduced here.

LO impulse
Before discussing the impulse computation at loop level, for completeness, we give a lightening discussion of the leading order impulse, which is purely transverse (see Eq. (6.7)) and basically the Fourier transform of the classical long-distance tree-level scattering amplitude An analogous leading order analysis has already been performed in the original work of KMOC [81]. The impulse kernel is given by In maximal supergravity, we discuss the scattering of non-identical scalars and include the BPS angle φ [79,128]. The factor s − |m 1 + m 2 e iφ | 2 = 2m 1 m 2 (cosh η − cos φ) expresses the prefactor in terms of the relative rapidity η = arccosh(σ) between the two massive states. Upon Fourier transforming to impact parameter space, we find the leading order impulse in both theories ∆p µ,(0) Our result in general relativity agrees with the earlier expressions derived in Refs. [95,96].

NLO impulse
At next-to-leading order, the structure of the one-loop classical amplitude is as follows M (1) (p 1 , p 2 , p 3 , p 4 ) = c II I II + c X I X + c tri,1 I tri,1 + c tri,2 I tri,2 where c i are the rational coefficients of the loop integrals. As we will see, at the classical level, the structure of the amplitude reveals that there is no difference between the conservative and radiative impulse at this order. The reason is that in general c II = c X so the box and crossed box integral appear in the combination I II + I X . Using the values of these soft integrals computed in Section. 5 and comparing to the values in the potential region given e.g. in Eqs. (4.54) and (4.59) of Ref. [79] we see that the difference between soft and potential region cancels in the sum. In addition the value of the triangle integral is identical in both regions.

N = 8 supergravity
Let us begin our one-loop discussion with the appropriate loop integrand for the scattering of non-identical scalars [79,128] where we use the same notation as in the tree-level analysis.
Equipped with this integrand, we can soft expand both the I II and I X integrals and plug the resulting expressions into the impulse kernels Eqs. (6.8) and (6.9) to obtain the transverse and longitudinal components. Let us discuss the transverse part first. From the real term, we only get a contribution from the cut of the box integral, where we have expanded the impulse numerator in our preferred basis (6.4) and truncated at the classical order. On the other hand, summing the box and cross-box is equivalent to a symmetrization of the graviton loop momentum. This effectively cancels the real parts of the box and crossbox integrals and sets on-shell the two matter propagators in the box and yields a purely imaginary term in the virtual part of the impulse kernel which exactly cancels the cut contribution so that we are left with consistent with the general expectation of Eq. (6.11). The fact that the real parts of the oneloop amplitude is zero shows that there is no transverse deflection of the black hole orbits at one loop in maximal supergravity. Ultimately, this is due to the no-triangle property [131][132][133][134] which was also linked to the non-precession of black hole orbits in Ref. [128]. A similar cancellation of the imaginary parts in the impulse was shown in Ref. [81] for the electromagnetic impulse. On the other hand, there is a contribution for the longitudinal impulse kernel. Naively, one could have guessed that the longitudinal impulse numerators 1 · u i , together with the cut conditions of the matter lines do not yield a contribution. However, in accord with the general expectation of Eq. (6.13), due to subleading terms in the soft expansion there remains a contribution. Alternative to the general expectation from Eq. (6.13), we can directly compute all diagrams in Eqs. (6.8) and (6.9) using the explicit results for all soft integrals from the previous section. This is quite instructive and will be used at higher loops as well. Upon soft expansion and IBP reduction of the real longitudinal impulse contributions, we find where the double line, again, denotes linearized soft propagators. The factor i/(16π 2 ) originates from the difference in normalization conventions between our soft integrals and standard Feynman diagrams. The value of the cut soft box is given in Eq. (5.12). Next, We perform the Fourier transform (2.11) to impact parameter space, using the results collected in appendix A to arrive at the final (purely longitudinal) result for the impulse Note that we can replace the soft velocities inǔ by the usual ones for free, since all superclassical pieces have canceled and we only need the leading in q terms.
We would like to mention that at one loop order, the impulse is the same in the soft and potential region and receives no radiative contributions at the classical order. This is owed to the fact that soft bubble integrals (that vanish in the potential region) only contribute at higher orders in the expansion and are therefore irrelevant. This also allows us to compare the one-loop impulse to the conservative result obtained from the scattering angle [79], finding full agreement. The extraction of the conservative impulse from the scattering angle is reviewed in appendix B, where it becomes clear that the (conservative) longitudinal impulse is completely dictated by lower-order information due to on-shell conditions.

General Relativity
Next, we consider general relativity. In principle, the same computational tools that led to all results in maximal supergravity are also applicable here. The only complication originates from a more complex loop integrand and more contributing soft master integrals.
Just like in maximal supergravity, we begin our discussion of the one-loop impulse with the relevant integrand, which is known from e.g. [128,129]. Notably, we find that at one loop in D = 4 there is no distinction between the conservative result and the full soft region and it suffices to consider the following covariant diagrams M (1) (p 1 , p 2 , p 3 , p 4 ) = c II I II + c X I X + c tri,1 I tri,1 + c tri,2 I tri,2 = c II where I tri,i is the triangle with matter propagator of mass m i , and the coefficients are We could proceed with the general relation in Eq. (6.13), however, here we explicitly check its validity by soft expanding Eqs. (6.8) and (6.9). Subsequently, we insert the explicit results for all soft integrals from the previous section. Let us begin by discussing the transverse part of the impulse, that has in principle two contributions, one from the virtual amplitude and one from the real piece. In the impulse formula, only the box diagram contributes to the cut and by the same reasoning as in N = 8, it just cancels the virtual boxes and we are left with only triangle contributions where the linearized triangle is given in Eq. (5.8) and we have again taken into account the standard normalization of Feynman integrals leading to the additional factor of i/(16π 2 ). The transverse impulse kernel therefore reads From the kernel we can easily calculate the transverse impulse via (2.11) ∆p µ,(1) 14) The remaining longitudinal impulse computation is essentially identical to the one in maximal supergravity, as only the box integral has the two-particle cut. Consequently, we simply have to replace the box coefficient c II in Eq. (7.8) by its pure gravity counterpart (7.11) ∆p µ,(1) so that the leading high-energy limit (σ 1) of the longitudinal impulse coincides between GR and maximal supergravity.
This concludes our one-loop calculation of the gravitational impulse within the KMOC formalism. We agree with previous results [45,97,98,135] that can also be obtained from the scattering angle only (see appendix B), since conservative and soft region results are identical in D = 4 up to classical order.

NNLO conservative impulse
Before deriving novel results in the full soft region, it turns out that we can test our two-loop setup by reproducing known results for the conservative dynamics from the KMOC point of view. This can be done by performing the calculation in the potential region, defined in Eq. (3.7). If we separate the impulse into conservative and radiative pieces the potential region only captures the conservative contribution, ∆p µ 1 ,cons. . In the potential region, the gravitons are off-shell and therefore there cannot be real (onshell) graviton emission. Hence, only the virtual integrals and the contribution from twoparticle cuts to the real impulse kernel in Eq. (6.14) survive. Furthermore the "mushroom" integrals are identically zero in the potential region, and hence do not contribute to any conservative quantity. All the remaining integrals can be evaluated using the differential equations of sections 3.3 and 5, although with modified boundary conditions appropriate for the potential region as originally described in Ref. [79]. Reproducing the conservative impulse [45] and the scattering angle [41,43] constitutes a highly nontrivial check of the most complicated parts of our assembly.

N = 8 supergravity
The two-loop integrand of maximally supersymmetric gravity is obtained by dimensional reduction of the massless integrand [130] with the following result [79] where (2 ↔ 3) instructs to add terms with p 2 and p 3 interchanged. This integrand is the complete quantum integrand for the supergravity amplitude, and hence is valid both for conservative dynamics, as well as in the full soft region discussed below. Note that the final four scalar diagrams are quantum suppressed in N = 8 because they are accompanied by the q 4 prefactor and only the H-diagram survives. Upon classically expanding the integrand (in |q| or ), subsequent IBP reduction to a set of soft master integrals, and inserting the appropriate values for both the virtual and cut pieces evaluated in the potential region, we find the impulse kernels u 1 ,cons = I (2) u 2 ,cons = 0 . (7.18) Note that three particle cuts are zero in the potential region, as the internal graviton lines are never on-shell. As expected, the superclassical terms in the transverse impulse kernel cancel between the virtual diagrams and two-particle cuts. Furthermore, the longitudinal impulse kernel does not receive a virtual contribution and vanishes due to a cancellation between the two-particle cuts of double boxes and crossed double boxes. Fourier transforming to impact parameter space yields ∆p µ,(2)

General Relativity
We can repeat a similar conservative two-loop analysis for pure gravity. The integrand is more complicated than (7.17) and contains additional diagrams but has already been constructed previously in Refs. [41,43], and reproduced by the present authors using generalized unitarity. It only involves the cubic graphs in the first row of Figure 5. Taking said integrand, expanding it in the classical limit, reducing to a minimal set of master integrals and inserting the appropriate potential region values of the master integrals [79] allows us to obtain the impulse kernels Fourier transforming to impact parameter space, we obtain The conservative impulse has a logarithmic divergence at high energies, corresponding to that in the scattering angle of Refs. [41,43]. By comparing to the maximal supergravity result we find that the coefficient is universal in agreement with Ref. [79]. We will come back to this point when considering the full impulse including radiation reaction. The longitudinal impulse is ∆p µ, (2) 1,u,cons = We note that the longitudinal part of the conservative impulse does not contain new information. Its purpose at a given order in G is to ensure that the energy transfer between the two particles is such that that they remain on-shell after transverse deflection at previous orders. In other words, the longitudinal impulse is the solution to the equation which must be satisfied at each order in G. This was used in [46] to obtain the O(G 3 ) longitudinal impulse in General relativity. In contrast our result in Eq. (7.24) follows from direct calculation and Eq. (7.25) serves as a check on our methodology. For the radiative impulse in maximal supergravity, we start from the full integrand in Eq. (7.17) obtained via dimensional reduction. Upon soft expansion of the integrand, subsequent IBP reduction to a set of soft master integrals, and inserting the appropriate values for both the virtual and cut pieces, we find the following impulse kernels After taking the Fourier transform in Eq. (2.11) and subtracting the conservative result in Eqs. (7.23) and (7.24) we obtain the following result for the radiative impulse in general relativity ∆p µ,(2) The structure is very similar to the result in maximal supergravity. However, unlike in the case of supergravity, where the longitudinal and transverse impulse was controlled by the same algebraic functions f i , in general relativity, the structure is different. Interestingly, the radiative transverse impulse is captured by algebraic functions f LS i that purely encode leading soft (LS) dynamics of gravitons [101]. In hindsight, the fact that in N = 8 supergravity the leading soft theorem also controls most of the longitudinal impulse can be understood as a consequence of the no-triangle property of theories with maximal supersymmetry [131][132][133][134]. In the absence of triangles the Weinberg soft factor exactly captures the contributions from all diagrams where the radiated gravitons are emitted from a matter leg. Thus the only new contribution can arise from the H diagram, which indeed produces the term with log σ+1 2 , as pointed out above. Using Eq. (B.14) we obtain the radiative contribution to the scattering angle in general relativity This result agrees with the computation by Damour in Ref. [100] via a linear response formula derived in Ref. [136]. This was later reproduced in [101] using a beautiful relation between the real part of the eikonal and the infrared divergence in its imaginary part. Such relation was proven for N = 8 supergravity by explicit computation and conjectured more generally.

LO radiated momentum
Besides the gravitational impulse, considered in the previous subsections, we are also able to compute the radiated momentum, both in general relativity and maximal supergravity. This observable starts at O(G 3 ) and is related to the energy loss, which has been the subject of our short letter [83] and we present them here just for completeness. In the KMOC setup, the radiated momentum can be obtained either directly by considering the expression in subsection 2.2, or from momentum conservation and the impulse on particles 1 and 2, For us, it was originally easier to obtain the radiated momentum directly, as it only involves the three-particle cut of two-loop diagrams at O(G 3 ) and therefore requires fewer terms in the full soft integrand. We found the following result in D = 4 where we define with the theory dependent coefficient functions f i (σ). This analytic structure is directly inherited from the longitudinal part of the radiative impulse, computed in section 7.4, by momentum conservation. As was pointed out in Refs. [13,18], the homogeneous mass dependence in Eq. (7.38) signals that the result is entirely specified by the probe limit m 1 m 2 . Note that the radiated momentum in Eq. (7.38) is purely longitudinal and yields the energy radiated as gravitational waves. In the center-of-mass (c.m.) frame of the hyperbolic motion, we find From the scattering result of Eq. (7.40), we obtain the energy loss for elliptic (bound) orbits via analytic continuation [11][12][13] of the result where E takes the same general form as Eq. (7.39) [83] and has the expected simplified ν dependence previously observed in Ref. [13]. From our perspective, this is simply inherited from the analytic continuation of the hyperbolic result. As stated previously, the explicit result for radiated momentum in Eq. (7.38) has been obtained in Ref. [83], where the theory specific coefficient functions in general relativity are the same that appear in the impulse computation, c.f. Eq. (7.34). The energy loss for a black hole scattering event can be expanded in small velocity v= √ σ 2 −1 σ and compared to known Post-Newtonian (PN) data, finding agreement with the result known up to 2PN [13,18,137]. We furthermore compared the small velocity expansion of the energy loss in elliptic orbits in Eq. (7.42) to the 3PN expressions for the instantaneous energy flux integrated over one orbit from Refs. [137][138][139][140][141][142][143][144][145] in the large eccentricity limit, i.e. to leading order in large J, again finding perfect agreement with the PN data where our results overlap. After Ref. [83] appeared, Bini, Damour, and Geralico, informed us privately of their computation in the small velocity limit up to O(v 15 ), also agreeing with our result. Additionally, Ref. [24] verified the low-velocity limit up to O(v 7 ) from a world-line EFT perspective, and Ref. [103] reproduced our result with full velocity dependence using methods similar to Ref. [83].
Finally, the radiated energy also appears in the tail term [13,146,147] of the O(G 4 ) radial action, which has been recently obtained by Ref. [51] by an independent computation. Comparing Eq. (7.40) to that expression, we find full agreement.

Comments on universality and relation to eikonal phase
It is interesting to study the high energy limit of the gravitational observables considered in this work. Famously, the leading order observables are universal in this limit [148]. Similarly, the gravitational deflection angle has been observed to have universal properties at O(G 3 ) [75,79,99,100,149], so it is natural to ask whether or not the same is true for the gravitational impulse.
Recombining the radiative impulse in Eq. (7.29) with the conservative impulse in Eq. (7.19), we can study the high-energy (σ 1) limit of the full result. The leading high-energy pieces cancel between radiative and conservative contributions, consistent with previous observations [99].
where it is interesting to note that the leading nonzero terms are independent of the BPS angle φ. On the other hand, taking the limit of our general relativity result we find The logarithmic high-energy divergence in the conservative impulse cancels, as expected [99,100], once radiation reaction effects are included. Interestingly, by comparing (7.44) to the maximal supergravity result (7.43), we find that the universality of the scattering angle (including radiation reaction) described in Ref. [99] does not hold for the full impulse (both in the transverse and longitudinal directions), due to a cancellation between the leading conservative and radiative contributions. However, computing the angle from the impulse requires taking into account products of lower-order terms which restore the previously observed universality of the high-energy limit of the scattering angle [99,100] χ (2) HE,GR = χ Here, we have written the result in terms of the angular momentum J = |b||p| = |b| M ν and note that subleading terms in the large σ expansion differ between both theories.
It is also interesting to consider the relation between the transverse impulse, ∆p 1,⊥ obtained in our calculation and the corresponding quantity derived from the eikonal approach in Refs. [101,103]. Comparing their conjecture for the real part of the two-loop eikonal phase to the transverse impulse (best seen from Eq. (7.31)), all but one velocity dependent factors agree (up to some overall scaling due to the distinction between the two quantities). The only difference is related to the s-dependent term in the first line of (7.31) and is due to the difference between the asymptotic impact parameter b and the eikonal one, b e (see Eq. (B.11)), which yields a correction proportional to χ (0) 3 . Indeed, rewriting our result for the transverse impulse in terms of the eikonal impact parameter we find in N = 8 supergravity where the χ (0) 3 correction cancels the terms proportional to s in the conservative part of the impulse of N = 8 in Eq. (7.19). In general relativity where we have denoted in red the correction due to the change of variables from b to b e . Taking the high energy limit holding b e fixed restores universality in the transverse impulse which now arises from the χ (0) 3 correction introduced by the alternative choice of impact parameter. Note however that the impulse now grows as σ 2 , rather than σ in the transverse part of Eq. (7.44). After this change of variables, we observe that, up to this order, the transverse impulse in N = 8 supergravity agrees with that in Refs. [101,103] which is given in terms of the Real part of the eikonal phase, δ(b e ), as 14 The same relation is true if we compare to the conjectured result for the eikonal phase in general relativity from Refs. [101,103], thus proving their conjecture. It would be interesting to verify it again directly by directly calculating the eikonal phase using the results from this work. This findings suggest a more general relation between the transverse impulse kernel and the eikonal phase which warrants further investigation.
Regarding the longitudinal part of the impulse, or the energy loss, we can study our full velocity dependent expressions in the ultra-relativistic limit σ → ∞ of Eq. (7.39). In N = 8 supegravity, we found the result for a single BPS angle [128] which has the structure of Eqs. (7.39) and (7.34) with the appropriate φ-dependent coefficient functions already defined in Eq. (7.28). The f i in Eq. (7.28) agree with our previous expressions [83] for φ = π/2. As in pure gravity, the ultra-relativistic limit σ → ∞ of the radiated momentum is controlled by the combinations f 1 and −f 2 + f 3 /2 with the leading high-energy term being independent of φ. Similarly, in general relativity Note that the apparent logarithmic divergence cancels in both cases. The high-energy limit of the general relativity energy loss can be compared to the prediction by Kovacs and Thorne [18], based upon the numerical probe calculation by Peters [14]. Our expression agrees structurally with [18], but disagree in the numerical coefficient. After Ref. [83] appeared, we were informed of the numerical computation of the high-energy coefficient by agreeing with our analytic result. Although the high-energy limit does not coincide ins Eqs.(7.51) and (7.52) in its rational prefactors (8 vs. 35/8), we noted in Ref. [83] that the ratio of the logarithmic (log 2) and non-logarithmic contributions is universal. Ultimately, it might not be too surprising that the radiated momentum depends on the theory content, since the number of massless messengers that can be radiated change between the two theories which suggests the bigger overall coefficient in maximal supergravity. Note that, in any case, our results are only valid for σ (GE cm /b) −1 , beyond which perturbation theory breaks down. For large enough σ, according to Eq. (7.52), the radiated energy exceeds the incoming energy, which, of course, is unphysical. Resolving this issue requires to account for destructive interference from multigraviton emissions, which cuts off the spectrum of gravitational waves at high-frequency 15 , as explained in Refs. [76,77,150].

Conclusions
In this work, we have employed the general formalism devised by Kosower, Maybee, and O'Connell (KMOC) to extract classical gravitational observables for the scattering of spinless black holes up to O(G 3 ), or third Post-Minkowskian order. This framework naturally includes radiative effects and goes beyond the much-discussed conservative binary dynamics. The presence of the gravitational interaction between the two massive black holes has two key physical effects, 1) a deflection and momentum shift on the individual black holes which is related to the gravitational impulse, and 2) the emission of gravitational Bremsstrahlung which is related to the radiated momentum. In our previous work, we have already presented the radiated momentum in general relativity and maximal supersymmetric gravity (N = 8 SUGRA). Here, we also present expressions for the impulse (which is related to the scattering angle) in both theories.
In order to render the general KMOC framework a practical computational tool, we have incorporated a number of ideas from collider physics to handle virtual Feynman integrals together with phase-space integration. Starting from generalized unitarity to construct loop integrands from gauge-invariant on-shell quantities, we employ the method of regions to facilitate the classical expansion. The relevant Feynman diagrams can be reduced to a minimal set of master integrals with the help of integration-by-parts identities. Using reverse unitarity we treat virtual and phase-space integrals on the same footing. At the end of the reduction step, we are left with a small set of independent integrals. In order to evaluate the master integrals, we solve a set of (canonical) differential equations, where the main complication is reduced to the computation to the boundary values of the master integrals. Making available all analytic expressions for the soft master integrals, we assemble the classical impulse and radiative momentum observables in both GR and maximal supergravity. Our results include the full radiation effects at O(G 3 ), but we have also reproduced the conservative gravitational impulse in GR, matching known results. From the impulse and the radiated momentum, we can derive the radiative scattering angle and the energy loss. Since our results are valid to all orders in the velocity, we are able to check against different regimes in the literature and compare against the Post-Newtonian computations by expanding our results in small velocity as well as against high-energy expectations. We find agreement with all known results where they overlap.
We have compared our results to the eikonal approach in Refs. [101,103] and found that the transverse impulse, when written in the appropriate variables is directly connected to the eikonal phase. This shows that the conjectured relation between the real and imaginary parts of the eikonal phase, put forward in Ref. [101], is also valid in general relativity. It also suggest suggests a more general relation between the transverse impulse kernel and the eikonal phase which warrants further investigation.
For the classical quantities considered here, we performed the full phase-space integration over all intermediate particles appearing in the KMOC setup, without imposing any further restrictions (or phase space cuts). In principle, the reverse unitarity method can also be adjusted to incorporate additional measurements on the final state particles [84][85][86]. One can envision a similar adaptation to the gravitational setup to measure more exclusive observables, such as the radiated energy spectrum or the angular distribution of the radiated momentum. These quantities depend on more scales and we leave their discussions to future work.
Besides attempting similar computations at higher orders in Newton's constant and the discussion of more exclusive observables, it would also be interesting to generalize the O(G 3 ) computation for spinning observables as well as to include tidal effects. Since all relevant master integrals are known, the only remaining change requires the construction of more complicated loop integrands. Since this step is very mature and can be automatized and streamlined via generalized unitarity, it should be possible to tackle such observables in the near future.

Acknowledgments
We are specially grateful to Paolo Di Vecchia, Carlo Heissenberg, Rodolfo Russo, and Gabriele Veneziano for stimulating discussions and for sharing a draft of their work [103] with us and for comments on our manuscript. We also thank Zvi Bern, Radu Roiban, Chia-Hsien Shen, and Mikhail Solon for helpful comments and collaboration on related projects; and Clifford Cheung

A Fourier transform formulae
The final step of the classical impulse and radiated momentum computation involves the evaluation of Fourier transform integrals of the kind These are simply related by differentiationf µ α (b 2 ) = −i∂f α (b 2 )/∂b µ , so we need only consider f α (b 2 ). It is convenient to use a Sudakov decomposition of the D-dimensional Lorentz vector q µ , where q µ ⊥ points in the (D − 2)-dimensional subspace transverse to u 1 , u 2 . With this parametrization the integral above becomes so that the delta functions localize two of the integration variables and force the momentum transfer into the D − 2 dimensional transverse subspace. Note that the impact parameter is always transverse so b 2 ⊥ = −b 2 ≡ |b| 2 . The remaining Fourier transform is elementary and we obtain

B Gravitational impulse and scattering angle
The gravitational impulse ∆p µ i describes the deflection of gravitationally interacting particles along the entirety of their hyperbolic trajectory. Another basic observable for such a scattering process is the scattering angle χ in the center-of-mass frame. In this appendix we describe the relation between these two observables.
Let us begin by considering a conservative scattering process. In this case, ∆p µ i, cons.
and χ cons. are exactly equivalent and contain the same information. It is well known how to relate the scattering angle in the center-of-mass frame to the impulse. The incoming and outgoing momenta have components and then solving for the three coefficients a i using the stated conditions. The result is which can be expanded perturbatively in G (∆p where χ (n) cons. ∼ O(G n+1 ). Note that magnitude of the c.o.m three-momentum p can be written in terms of the quantities used in the rest of the paper as follows . As a side comment, note that this formula nicely explains the structure of the conservative result in maximal supergravity, where the one-loop scattering angle χ (1) cons. is zero [79], which can be attributed to the "no-triangle" property of this theory [128]. In particular, this implies that the one-loop impulse is purely longitudinal, and the two-loop impulse purely transverse, in agreement with our explicit calculation.
Let us point out that the definition of the impact parameter b is chosen in terms of the initial momenta such that it satisfies b · p 1 = b · p 2 = 0. This choice, however, breaks the symmetry between the initial and final state in the conservative process (i.e. time reversal invariance). Instead, one could choose to modify the impact parameter as follows such that it more symmetric between the initial and final state b eik · p i = −b eik · (p i + ∆p i ). Due to this modification, the magnitude changes |b eik | = |b| cos χcons. 2 , and b eik can be recognized as the so-called "eikonal impact parameter" that naturally arises from semiclassical considerations in the eikonal approach [149]. Almost by definition, the impulse is purely transverse in these variables, that is, proportional to b µ eik , This form of the conservative impulse is the most convenient to compare to results from the eikonal method. More generally, radiative effects imply that the c.o.m of the binary is not an inertial reference frame, so the relation between the scattering angle and impulse is not as straight forward in the non-conservative setup. Momentum conservation ∆p µ 1 + ∆p µ 2 = −∆R µ illustrates that the radiative dynamics are a multi-body process and that the center of mass recoils. One can still write down an analogous formula to Eqs. (B.7)-(B.9) that relates the radiated momentum and radiative scattering angle to the impulse, although generically, the precise meaning of the angle becomes somewhat obscure in the presence of radiation. Said differently, including radiation, the kinematics is that of a five-point process for which there are two t-channel invariants, which translates to two angles. At O(G 3 ), however, the leading recoil effects do not yet affect the transverse part of the impulse and the angle is still given by At two-loops, the unitarity relation includes more terms with both two-and three-particle cuts and is given by (C.6) Similar generalizations also hold to higher-loop order which are, however, irrelevant for the discussion in the present work. Note that these relations were crucial in order to simplify the KMOC kernels and proof the reality properties of the classical observables of interest which led to the results in section 6.

Cutting rules
Alternatively, our computation of phase space integrals can be based on Cutkosky's cutting rules [152], which can be applied to individual diagrams, rather than the full amplitude. For us, the application of the cutting rules is twofold. First, we make use of them to simplify the KMOC kernels in section 6. Second, we extensively utilize these rules to deduce phase-space integrals from virtual integrals. In fact, we actually use the cutting rules for soft-expanded integrals, where massive propagators are linearized and have the form i/(2u i · i ) while their cut versions have the form 2πθ( 0 i )δ(2u i · i ), but the usual proofs of cutting rules, e.g. using Veltman's largest time equation [153], carry through unchanged.
For illustration purposes, we consider a field theory with two massive complex scalar fields Φ 1 and Φ 2 , and a light scalar field φ, whose Lagrangian density is From the point of view of individual Feynman diagrams, the difference between this theory and the gravitational theory considered in the body of the paper is that the latter diagrams contain additional numerator, which leave the discussion below unchanged. The Feynman vertices are always −iκ, for Φ † 1 Φ 1 φ, Φ † 2 Φ 2 φ, and φ 3 couplings. The propagators with momentum k are i/(k 2 − m 2 1 ), i/(k 2 − m 2 2 ), and i/(k 2 ) for the fields Φ 1 , Φ 2 , and φ, respectively. We use thick dashed black lines to denote heavy scalar particles Φ i , and thin dashed lines to denote light scalar particles φ. For discussing loop integration, it is convenient to introduce the notion of "scalar integrals" which have unit numerators with all factors of i removed. For example, a simple two-loop diagram for Φ 1 + Φ 2 → Φ 1 + Φ 2 scattering an the corresponding scalar integral are The scalar diagram on the r.h.s. of (C.8) is essentially the same as the Feynman diagram on the l.h.s., except that dashed lines are replaced by solid lines to indicate that every propagator is understood to be without the i factor in the numerator, and every vertex is simply κ rather than −iκ.

Cutting rules and translation to scalar integrals
Having defined the virtual scalar and Feynman diagrams, we introduce unitarity cuts of Feynman diagrams, where a vertical blue dashed line highlights the cut propagators, e.g. .

(C.9)
Every cut propagator is given by the simple replacement rule i k 2 − m 2 −→ 2π θ(k 0 ) δ(k 2 − m 2 ), (C. 10) which simultaneously imposes the positive energy and the on-shell condition. Instead of performing the full loop integration, in the presence of the on-shell conditions, the remaining integrals are over the Lorentz-invariant phase space of the on-shell states exchanged across the cut. According to the Cutkosky rules, the uncut propagators and vertices on the left hand side of the cut are given by their usual expressions, while those on the right are given by the complex conjugates of their usual expressions. In this notation, the "rightmost" cut is always the same as the uncut diagram up to certain factors, while the "leftmost" cut is equal to the conjugate of the uncut diagram, When using solid lines, the propagators and vertices on either side of the cut are without any factors of the imaginary unit i, while cut propagators are still given by the right hand side of Eq. (C.10), 2π θ(k 0 ) δ(k 2 − m 2 ). Using the cutting rules, the sum of all cuts in a given channel is zero, Translating to scalar integrals without i factors, this reads as (C.17) A clear difference from the above equation for the planar double box, (C. 16), is that the double cut contribution (i.e. the last term) is multiplied by (−1) rather than (−2), because the crossed double box has only one double cut. The u-channel double box has no double cut or triple cut. The u-channel crossed double box has only one triple cut, which evaluates to twice the imaginary part of the virtual integral, analogous to Eq. (C.14).
Note that Eq. (C.15) is also valid for field theories other than scalar φ 3 theory, so Eq. (C.15) holds with numerators multiplying the loop integrand of every diagram in the relation, and therefore also for scattering amplitudes and their unitarity cuts! 16 16 Depending on the properties of the numerator, diagram symmetries might be broken so that one cannot directly obtain Eq. (C.16). After IBP reduction, however, we will mostly be dealing with master integrals without numerators for which the symmetry is restored.

D Master integrals
As explained in Ref. [79] at two-loops there are three irreducible families of master integrals, the III, H and IX families. These families contain a total of 20 unique linearized-propagator master integrals. As explained in Ref. [79] at two-loops there are three irreducible families of master integrals, the III, H and IX families. All results can be found in a computer readable format in the ancillary files of this arXiv submission .
To be more precise, these are integrals that scale like half-integer powers of |q| = −q 2 , before being multiplied by normalization factors like −q 2 .