Classical gravitational scattering from a gauge-invariant double copy

We propose a method to compute the scattering angle for classical black hole scattering directly from two massive particle irreducible diagrams in a heavy-mass effective field theory approach to general relativity, without the need of subtracting iteration terms. The amplitudes in this effective theory are constructed using a recently proposed novel colour-kinematic/double copy for tree-level two-scalar, multi-graviton amplitudes, where the BCJ numerators are gauge invariant and local with respect to the massless gravitons. These tree amplitudes, together with graviton tree amplitudes, enter the construction of the required $D$-dimensional loop integrands and allow for a direct extraction of contributions relevant for classical physics. In particular the soft/heavy-mass expansions of full integrands is circumvented, and all iterating contributions can be dropped from the get go. We use this method to compute the scattering angle up to third post-Minkowskian order in four dimensions, including radiation reaction contributions, also providing the expression of the corresponding integrand in $D$ dimensions.


Introduction
The goal of this paper is to propose an efficient method to compute the scattering angle in a collision of two black holes. To achieve this, we bring together a number of crucial ingredients which we now describe.
We begin by observing that in a scattering process, the heavy black holes can be represented as pointlike particles, exchanging momenta which are much smaller than their masses. In practice, in order to extract classical physics such as a deflection angle, one rescales the momentum of the exchanged gravitons q as [1] q = k, with the wavevector k being kept fixed while taking the classical limit → 0. This motivates the use of effective field theory applied to gravity [2], and in particular of a heavy-mass effective field theory (HEFT) [3][4][5][6][7], which is the appropriate tool if one wishes to describe interactions of particles where the exchanged momenta are much smaller than their masses. Working from the outset with a simplified theory is our first crucial ingredient -the great advantage is that this avoids the soft (or heavy-mass) expansions of the integrands that would be obtained starting from Einstein gravity coupled to matter. The scattering angle is then extracted from the computation of the elastic scattering amplitude of the two black holes, an approach that was initiated fifty years ago in [8] (specifically for the computation of the Newton potential), and then pursued at 2PM order in [9][10][11][12] in conjunction with the unitarity method [13,14]. This has been applied successfully in several works in general relativity at 3PM [15][16][17][18][19], including radiation [20][21][22][23], 4PM [24], and also in theories with higher-derivative interactions [25][26][27][28][29] to compute both corrections to Newton's potential and deflection angles. An important simplification stems from the fact that only terms in the four-point scattering amplitude with a discontinuity in the q 2 -channel need to be computed, as these contribute to long-distance effects [2,30], which makes the unitarity-based method particularly well suited for this task. These amplitude techniques have emerged as potential alternatives to a host of other methods such as the effective one-body formulation [31][32][33][34] and worldline theories [35][36][37][38][39][40], which have led to a wealth of important results  also including spin effects [70][71][72][73][74][75][76][77][78][79][80][81][82][83][84].
In order to be able to compute efficiently, one needs compact expressions for the tree amplitudes entering the unitarity cuts, and here comes the second key ingredient of our procedure. In [85] we have developed a systematic approach to derive compact, Ddimensional expressions of HEFT amplitudes with two massive scalars and an arbitrary number of massless gravitons at leading order in an inverse mass expansion. This is based on a new form of the colour-kinematics/double-copy duality [86,87] for heavymass effective field theories, initially proposed in Yang-Mills and gravity in [88], and further developed in [89]. The key advantage of this method, compared e.g. to the earlier work of [90], is that it provides BCJ numerators that automatically satisfy the Jacobi relations and crossing symmetry. Furthermore, only a subset of all possible cubic graphs contribute to these numerators, which in addition are manifestly gauge invariant (since they can be written in terms of field strengths) and also local with respect to the massless gluons or, in the double-copied theory, gravitons. This is in contradistinction with the double copy for Yang-Mills amplitudes, where the BCJ numerators are in general not gauge invariant. Thus, our novel double copy provides us with particularly compact tree amplitudes which can be fed into the cuts to construct equally compact integrands.
The third key simplification in our approach is of a diagrammatic nature. Specifically, we propose that the scattering angle can be computed from a subset of all possible diagrams in the HEFT, namely only two-massive particle irreducible (2MPI) diagrams, discarding all the others. From this 2MPI amplitude M 2MPI HEFT , one can then construct a corresponding HEFT phase δ HEFT by simply transforming it to impact parameter space, and finally extract from it the scattering angle χ as where J is the total angular momentum of the system. This procedure should be contrasted with the usual eikonal method [91][92][93] -there the S-matrix is believed to exponentiate in impact parameter space, but starting from two loops, the computation of the eikonal phase becomes a delicate task. Conventional unitarity tackles the computation of S-matrix elements, and to extract the eikonal phase at a certain perturbative order one has to remove all terms which reconstruct the exponentiation of the phase at lower loop orders in perturbation theory. Clearly it is highly desirable to be able to tackle the computation of the phase directly without having to evaluate a plethora of iterating terms which pollute the phase, and this is what our proposal provides. Important work in this direction was done in [24], where remarkably the conservative part of the potential at 4PM was computed from the radial action. Related ideas were presented in the recent interesting paper [94], which proposed a different exponential representation of the S-matrix inspired by the WKB formalism applied to quantum field theory, differing from the eikonal one. In particular, the authors of [94] were able to write down expressions of the iterating terms which have to be removed from the complete S-matrix elements (whose calculation is still needed) in order to compute matrix elements of N, where S := e i N . This still requires knowledge of complete amplitudes, but it cleanly explains their structure. Instead, we propose an approach where by computing only a subset of all possible diagrams one can efficiently extract the scattering angle avoiding entirely the computation of iterating diagrams.
We test our proposal by computing the scattering angle at two loops, or 3PM, also including radiation reaction terms. Our seed tree amplitudes are D-dimensional, hence we can produce integrands valid in an arbitrary number of dimensions. We use integration by parts (IBP) relations [95][96][97][98] to expand our result in a basis of seven master integrals (in the four-dimensional case), which we then integrate using the method of differential equations [99] in conjunction with a study of the boundary conditions based on the asymptotic expansion of Feynman integrals [100,101]. Our result for the deflection angle agrees with [15,18,19,102], and specifically the radiation reaction part agrees with [18,21,103,104].
The rest of the paper is organised as follows. Sections 2 and 3 contain a brief review of the kinematics of the scattering process and of the tree-level HEFT amplitudes derived in [85], which will then be needed to perform our two-loop unitarity cuts, respectively. In Section 4 we discuss in detail the HEFT expansion and our diagrammatic method, in particular introducing the 2MPI diagrams needed for the computation of the scattering angle at one and two loops. Sections 5 and 6 describe the calculation of the HEFT amplitudes at one and two loops. In particular, in Section 6 we provide the integrand of the two-loop HEFT amplitude in D dimensions in terms of a basis of eight master integrals, with only seven contributing around four dimensions. We evaluate the master integrals in Section 7, and present our final result for the (integrated) 2MPI HEFT amplitude, the HEFT phase δ HEFT and the scattering angle in Section 8. In Section 9 we present a conjecture valid in D dimensions for the single diagram contributing to the probe limit for an arbitrary number of loops. Finally, we present our conclusions and give a brief outlook in Section 10. Two appendices complete the paper: in Appendix A we present a study of the asymptotic behaviour of our master integrals and the associated boundary conditions, with explicit examples for our one-and twoloop master integrals; and in Appendix B we show explicitly the integrand of the four 2MPI cut diagrams contributing to the two-loop 2MPI HEFT amplitude.

Kinematics of the scattering process
We begin by discussing the kinematics of the 2 → 2 scattering process. We choose the particle momenta so that p 2 1 = p 2 4 = m 2 1 , p 2 2 = p 2 3 = m 2 2 , where particles 1 and 2 are incoming, while particles 4 and 3 are the outgoing deflected particles (with mass m 1 and m 2 , respectively). We parameterise the external momenta in the centre of mass frame, hence the momentum exchange q is spacelike, q µ = (0, q). Specifically, the momenta will be parameterised as follows: (2.1) Since we are considering elastic scattering, we have where p · q = 0 due to momentum conservation. It is also easy to see that which means that only thep i momenta are strictly orthogonal to the momentum transfer q, while e.g. p 1 ·q = q 2 /2. This will be important in the soft/heavy-mass expansion discussed later. Furthermore, our Mandelstam variables are defined as We also definep Furthermore, it is useful to introduce the relativistic velocities of the particles via p i := m i v i , with i = 1, 2 and v 2 i = 1, and their scalar product which is related to the centre of mass energy as Similarly, one can definep i :=m ivi andȳ :=v 1 ·v 2 , withv 2 i = 1. Finally, the deflection angle χ is related to the external kinematics as where P = p 2 + q 2 /4 is the common modulus of the incoming and outgoing threemomenta. This quantity can also be expressed as 9) and using this relation we can rewrite (2.8) as It is also useful to introduce the conserved total angular momentum where b is the asymptotic impact parameter of the scattered objects. Its inverse J −1 is commonly used as an expansion parameter of physical observables such as the deflection angle χ (see for instance our final result (8.9) for this quantity up to 3PM).

Relevant HEFT amplitudes from the new double copy
In this section we summarise the tree-level graviton-matter amplitudes to leading order in an inverse mass expansion, which enter the computation of gravitational scattering of two heavy scalars through unitarity cuts. We call these quantities "HEFT amplitudes". These have been derived and described in detail in the companion paper [85] and correspond to the leading terms in the heavy-mass expansion, which are the only ones relevant for classical physics. In particular, in [85] it was shown that they exhibit a novel double-copy structure with manifestly gauge-invariant numerators. This feature, combined with their very compact forms, makes them ideally suited for the purpose of this paper.
The three-point YM-matter amplitude are directly obtained from the Feynman rule and the GR-matter amplitude is the square of the YM-matter amplitude, At 2PM order we need the two-scalar two-graviton amplitude, which can easily be obtained using our double copy from the corresponding YM amplitude: where F denotes the field strength F µν i = p µ i ε ν i − p ν i ε µ i . Up to four points, there is only one cubic graph for each of these amplitudes. 1 At 3PM we need amplitudes up to five points. In this case there are two cubic graphs for the colour-ordered YM-matter amplitude and three cubic graphs for the GR-matter amplitude. Then the colour-ordered YM-matter amplitude and the GRmatter amplitude in the novel colour-kinematics duality form are and, by the double copy,  1 In the following we omit a factor of (−iκ) n from all graviton-matter and graviton amplitudes. These factors will be reinstated from Section 5 onward. Note that in our conventions κ 2 := 32πG N .
An advantage of this double copy is that the BCJ numerator can be written in a manifestly gauge-invariant form. For the five-point case The operator L(2, 3, 4, . . . , n − 1) was introduced in [88,89] as a group algebra element [105,106] and in [107] as a free Lie algebra element, where P (j 1 j 2 j 3 ...jm) denotes the cyclic permutation j 1 → j 2 → j 3 → · · · → j m → j 1 . For instance In [85] we also presented the six-point HEFT amplitudes, which are key ingredients for computations at 4PM order which were completed only very recently for the conservative dynamics [24,39]. Using the method of [85] it is straightforward to produce the HEFT amplitudes needed at higher PM orders. For instance, the n-point HEFT amplitude enters the probe-limit contribution at (n − 3)-loop order (or (n − 2)PM), as explained in Section 9. This contribution may also be computed from the geodesic motion of a test particle, whereas the contributions beyond the probe limit only require HEFT amplitudes up to n − 1 points. Hence, all ingredients needed at 5PM are already available with the currently known HEFT amplitudes.
The HEFT amplitudes scale universally with the mass as m 2 , where terms proportional to delta functions (arising from the iε prescription) have been dropped in the expansion of the full amplitudes 2 -importantly, these contact terms are captured by products of lower-point HEFT amplitudes and delta functions, and are essential to construct loop integrands. This is at the heart of our HEFT expansion and will be discussed in detail in Section 4.2. In the following we denote the HEFT amplitudes as A n (2 . . . (n−1), v), where from now on we omit the superscript in A as in this paper we only focus on Einstein gravity. The mass dependence of the HEFT amplitude can be read off from the label v orv (e.g. A n (2 . . . (n−1), v) ∝ m 2 , and A n (2 . . . (n−1),v) ∝m 2 ).
We conclude this section with a few comments. Our HEFT amplitudes are valid in a generic number of dimensions D and are manifestly gauge invariant, which brings a number of advantages. First, the unitarity cuts we perform are D-dimensional, avoiding the risk of missing terms by working in four dimensions, and lead to integrands that are valid in D dimensions. Second, manifest gauge invariance allows us to obtain very compact expressions and perform intermediate D-dimensional state sums without introducing spurious poles at any stage. Finally we will present all integrands in D dimensions, while performing the integrations in D = 4 − 2ǫ.

Computational strategy
We now propose an approach that addresses directly the computation of a HEFT scattering phase δ HEFT from which we can derive the bending angle as where J is the total angular momentum of the system introduced in (2.11). This phase is defined as the sum of all two-massive particle irreducible (2MPI) diagrams contributing to the scattering amplitude 3 , where J = 4 (p 1 ·p 2 ) 2 − p 2 1 p 2 2 = 4m 1 m 2 y 2 − 1. Conversely, two massive particle reducible diagrams automatically factorise in impact parameter space, effectively trivialising the exponentiation of the classical HEFT phase. We note that our HEFT expansion makes this manifest at the diagrammatic level, and all iterating, two-massive-particle reducible diagrams can simply be dropped.
Our proposal is very reminiscent of the non-abelian exponentiation theorem of [108,109], which allows the computation of the exponent w of the expectation value of the Wilson loop W = e w in gauge theory from a subset of all diagrams, namely those with a "maximally non-abelian" colour factor. This is in complete analogy with the Wilson loop diagrams contributing to the exponent w, which are also two Wilson line irreducible diagrams, in the sense that cutting two Wilson lines does not break the diagram into two disconnected components, see e.g. [110] for recent applications and [111] for an extension to multiple Wilson lines.
It is useful to compare more closely our approach to that of the recent interesting paper [94]. As mentioned in the Introduction, in that paper it was suggested to write the S-matrix in an exponentiated form as S = e i N , where N is a hermitian operator.
Anticipating our story, we will find that up to two loops, the two-body matrix element of the operator N introduced there is identical to the real part of our 2MPI amplitude -compare our (8.4) to (3.27) of that paper. Our M 2MPI HEFT has in addition an imaginary part which at 3PM arises entirely from the radiation reaction terms, but the real part is identical to the elastic matrix element of the operator N. The advantage of our approach is that it provides a practical way to compute δ HEFT directly, from which one can then extract immediately the scattering angle using (4.1). We will expand on the comparison of our results to those of [94] at the end of Section 4.4.

Diagrammatics of the HEFT expansion for classical dynamics
We now discuss how the HEFT provides us with a systematic expansion of L-loop amplitudes at the diagrammatic level. This will allow us to cleanly isolate contributions to classical physics without the need of expanding complete amplitudes/integrands, which in general can be rather cumbersome. To do this expansion one rescales all the soft graviton momenta ℓ i → ℓ i and the momentum transfer q → q [1], which in turns is equivalent to a (refined) expansion in terms of the inverse masses 1/m 1 and 1/m 2 of the binary system.
Concretely, in order to construct the loop integrands we will employ the unitarity method, where only massless graviton legs are cut. We only need to compute terms in the four-point amplitude which have discontinuities in the q 2 -channel, which contribute to long-distance effects [2,30]. As we will explain, we only need the leading-term in the heavy-mass expansion -that is the tree-level HEFT amplitudes reviewed in Section 3, as well as pure graviton amplitudes. For concreteness we work with massive scalar fields, and the result at leading order in the mass expansion is independent of the spin of the massive particle.
We now describe the expansion of the gravity-matter amplitudes in HEFT. The first observation is related to joining two different HEFT amplitudes into a larger tree. In gravity we have to sum over all possible orderings of gravitons, as implemented in (4.3): The scalar propagator in (4.3) is proportional to We can rewrite this in two ways. The first one is: where Q is the sum of the soft graviton momenta, and we have performed a Taylor expansion, with the dots indicating terms of O(Q). As noted in [17], such a Taylor expansion does not lead to a systematic expansion. Indeed, since (p + q) 2 = p 2 , it follows that p·q = −q 2 /2. This implies that p, when dotted with other momenta, does not have a homogeneous degree in .
A better choice is to use thep =mv variables introduced in (2.1), in which case the propagator takes the form This alternative expansion is useful becausep·q = 0, which follows from the equality (p + q 2 ) 2 = (p − q 2 ) 2 . This is important when considering the combination of terms in (4.3), because it avoids producing terms that carry higher powers of . Crucially the sum of the two diagrams in (4.3) contains the factor 1 2p·Q a + i ε where Q a = ℓ a 1 + · · · + ℓ a i := ℓ a 1 ···a i and This term is of O(m 3 ): two factors ofm 2 from the HEFT amplitudes, and one factor ofm −1 from the delta function. A Taylor expansion of the full tree-level amplitude, without taking into account the iε prescriptions, would miss such terms and produce only the HEFT amplitude with i + j gravitons.
This procedure works recursively for our matter-graviton amplitudes. A generic amplitude is then decomposed as where P(n − 2, h) denotes the partitions of the n − 2 gravitons into h non-empty subsets, and the summation is taken over all the partitions with h = 1, . . . , n − 2. P j denotes the j th subset of graviton indices of a given partition P with length i j and total momentum ℓ P j . The term without any delta function (h = 1) is the HEFT amplitude, and the dots stand for terms which are subleading in the heavy-mass expansion, which we ignore as they correspond to quantum corrections. Each term in the expansion is of order O(m h+1 2 ) and carries a uniform power of , facilitating a systematic separation of classical and hyper-classical terms in our computations.
The expansion (4.8) produces integrands involving inverse powers ofp i · l i and δ(p i · l i ) which is crucial to produce a clean expansion, and in particular leads to loop integrals that have a trivial dependence on q through a universal power that is determined by power counting and otherwise are non-trivial functions of the dimensionless quantityȳ :=v 1 ·v 2 only. This choice is also natural from the point of view of the Fourier transform to impact parameter space, if the (D − 2)-dimensional subspace of q is defined throughp i ·q = 0 and not via p i ·q = O(q 2 ) as often done.

Tree-level examples
We now give a few concrete examples of the HEFT expansion. The first case we consider is the four-point tree-level amplitude, expanded in the heavy-mass limit. The exact four-point amplitude with two massive scalars is Its expansion in terms of 1/m 2 is where we have definedp i =m ivi withv 2 i = 1 and i = 1, 2. The delta function term is of O(m 3 2 ), while the HEFT amplitude is of O(m 2 2 ), and the dots represent terms of O(1) in the mass.
Similarly, the expansion of the five-point amplitude in the heavy-mass limit is where A 5 (234,v 2 ) is given in (3.4). It is easy to see that gauge invariance, locality and crossing symmetry of the gravitons are manifest.

Loops
At one loop, the amplitude integrand can be expanded as Here the black dots represent three-point HEFT amplitudes while higher-point ones are depicted using grey blobs. The wavy lines correspond to the graviton lines, the dashed red lines denote unitarity cuts δ(l 2 i ), while the continuous vertical red lines represent the delta functions δ(2p i · l j ) of the massive lines -these are not unitarity cuts but arise from linear propagators as in (4.7). More concretely the four cut diagrams are obtained from the product of two copies of the expanded four-point amplitude (4.10) -one for particle 1 (lower line), the other for particle 2 (upper line).
The first diagram is two massive particle reducible and gives a hyper-classical contribution (O( −1 ) compared to the classical terms). The last diagram gives rise to quantum corrections, which we will not compute in this paper. Only the remaining two diagrams are relevant for classical scattering, and will be computed in Section 5.
At two loops, the integrand is expanded in terms of a hyper-classical, classical and quantum part. The hyper-classical graphs arē As we will discuss in Section 4.5, these diagrams factorise in impact parameter space, and hence exponentiate. Therefore they do not need to be computed and we can drop them directly at the diagrammatic level.
The classical pieces are obtained from the 2MPI diagrams: and There are additional graphs to consider, usually called "radiation-reaction contributions". These arē In Section 5 we will compute explicitly all these 2MPI diagrams (including radiation reaction). We will also show that the last cut diagram in (4.15), which does not have a three-graviton cut, does not add any new contribution to the four-dimensional 2MPI amplitude, and hence to classical physics.
We also note that the following graph has the correct mass scaling to potentially give classical contributionsm but as we argue in the next subsection it can be dropped since it is not 2MPI, similarly to the diagrams in (4.13).
The heavy-mass expansion can in principle be extended to compute quantum corrections, leading to the following diagrams at two loops: where we do not draw diagrams which are related by crossing of the external legs. In this paper we will not be concerned with such quantum contributions.
At higher loops, more graph topologies need to be calculated, but we stress that only 2MPI diagrams need to be evaluated. As an important example of higher-loop (L > 2) diagram, we will consider the probe limit m 2 ≫ m 1 . The relevant graph is ), for which we will present an all-loop conjecture valid in D dimensions in Section 9.
Finally we can compare our approach to that of [94]. It is instructive to focus on Eq. (2.17) of that paper. In our procedure, the iterated diagrams in the second line of that equation are simply not included (hence there is no need to subtract them). Instead, the cut diagram in the first line involving two cut matter lines and one cut graviton is included in our radiation reaction diagram, and is precisely responsible for the imaginary part in the two-loop phase δ HEFT computed in Section 8. Had we removed it, our procedure would have led to an identical result to the two-loop matrix element of the operator N of [94].

To bar or not to bar: factorisation and exponentiation
In the previous sections we have seen that in our HEFT expansion it is natural to use the variablesm i andp i instead of m i and p i , and that each diagram is weighted by a unique factorm imj which allows for an immediate classification into classical, hyper-classical and quantum contributions at any loop order. This might appear unnatural since the parameterm 2 i = m 2 i − q 2 /4 mixes classical and quantum terms. For the classical contributions this is not an issue since the expansion would only lead to subleading quantum corrections which we can drop, but expandingm i within hyper-classical contributions in general produces classical feed-down terms.
However, we extract the deflection angle from δ HEFT that appears in the exponentiated form of the 2 → 2 scattering amplitude in impact parameter space (IPS) with where the bars indicate that the HEFT scattering phases at the respective loop orders are expressed in terms of barred variables. Now at the level of phases only classical and possibly quantum corrections appear, and all hyper-classical contributions have been repackaged in the exponentiated form of M. Since we are only interested in classical physics, we can drop all quantum corrections, and at this stage replace all quantities with bars by unbarred quantities: Hence in the next sections we will evaluate the appropriate sum of 2MPI diagrams that gives δ HEFT using unbarred variables.
In the remainder of this section we want to show that in the HEFT expansion the exponentiation is manifest at the diagrammatic level, which is extremely useful since this allows us to "subtract" hyper-classical or other unwanted contributions simply by not computing them. This is in contrast to well-known eikonal exponentiation, which usually requires the computation of complete amplitudes, even including quantum corrections [115].
In order to achieve this diagrammatic rearrangement we need to align the Fourier transform to IPS with the expansion in terms ofm i ,p i . Therefore, we define the IPS form of an amplitude as where the Jacobian isJ = 4m 1m2 ȳ 2 − 1 = 4 (p 1 ·p 2 ) 2 −p 2 1p 2 2 . The crucial observation is now that any two massive particle reducible diagram is a convolution integral in momentum space, which after Fourier transform to IPS becomes a simple product of sub-diagrams in IPS: where b = | b |. Therefore, with the choice (4.23) any massive two particle reducible diagrams factorise exactly in IPS 5 .
Important examples of this are the first diagram in (4.12) and all diagrams in (4.13). This makes it clear that the hyper-classical diagram in (4.12) (first diagram that figure) is proportional to the square of the tree-level diagram M (0) , while the first diagram in (4.13) is a triple convolution and hence is proportional to the third power of the tree diagram. For this class of diagrams involving only trivalent vertices one can actually perform an all-loop computation. Taking into account our conventions and appropriate symmetry factors, the diagrams can be resummed as which agrees of course with the well-known result for the leading eikonal found by [93].
Focusing on the second and third diagram of (4.13), we see that in IPS they factorise into a tree diagram, and the sum of two irreducible one-loop diagrams (diagrams two and three in (4.12)) which contribute to classical physics, hence in IPS these reducible two-loop diagrams factorise into 2MPI . On the other hand the diagram (4.17) factorises in IPS into a tree diagram and the fourth diagram of (4.12) (which is a quantum correction), i.e.
Hence, up to two-loop order we can combine the contributions from all two-loop reducible diagrams and all one-loop and tree diagrams as Adding now also contributions from all other two-loop 2MPI diagrams (4.14), (4.15) and (4.16) we find, up to two loops, the exponentiated form of the S-matrix: We can now define the exponent of this equation as δ HEFT , which is a complex quantity. We expect that in general this phase is given by the sum of 2MPI diagrams even beyond two loops, with the deflection angle given as χ = −∂Re(δ HEFT )/∂J, where J is the total angular momentum. Importantly, δ HEFT has only classical and quantum terms but no hyper-classical terms. Since we are only interested in classical physics, we can now replace all barred quantities by unbarred ones, and we can also drop M 2MPI -which we actually never have to compute.

The one-loop HEFT amplitude
In this section we compute the 2 → 2 scalar HEFT amplitude at one loop. In the HEFT approach, the amplitude can be expanded in powers of the massesm 1 andm 2 .
At one loop, this expansion is given by Reinserting powers of , one immediately sees that the first term corresponds to the hyper-classical part of the amplitude while the second and third to the classical part.
Remaining quantum corrections are beyond the scope of this paper.
The hyper-classical part is computed from the following HEFT diagram: As discussed earlier in (4.7), the top part of the figure, involving two three-point vertices, represents the sum of the two possible diagrams Each of the three-point vertices here is O(m 2 2 ), however similarly to [93] in the sum we obtain the combination As a consequence, the massive propagator in (5.2) is effectively put on shell (as depicted by the continuous red line), and hence this term is of O(m 3 2 ). Note that this corresponds to the first term in (4.10) in the expansion of the full four-point tree amplitude. Furthermore, symmetrising the integrand over the two loop legs leads to the remaining delta function in (5.2). In fact the symmetrisation argument of [93] is not needed to argue for the second delta function. In order to obtain the cut integrand for the hyper-classical contribution we simply need to multiply the first term of (4.10) with a similar factor for the lower line in the diagram (5.2), giving directly a product of two delta functions. The factor 1/2! in (5.2) is simply the symmetry factor required in this two-particle cut where identical particles cross the cut.
The sum over internal graviton polarisations in (5.2) was performed using 6 where by f | v , f | η we denote replacing ε µ i ε ν j by v µ v ν and η µν , respectively. We do not evaluate this diagram (or any other two-particle reducible diagrams), since it gives only hyper-classical contributions which we subtract simply by not evaluating them. From now on we will only consider diagrams that make classical contributions, therefore we can drop the distinction between barred and unbarred quantities which becomes irrelevant as the difference will only create quantum terms, as explained in Section 4.5.
The classical part of the amplitude is evaluated from the 2MPI diagrams as is obtained from the diagram below: where the dotted line represents the unitarity cut. The lower half of the graph is the two-massive two-graviton tree amplitude in the HEFT, which is of O(m 2 1 ) and universal for all spins of the massive particle. Hence the diagram at hand is of O(m 2 1 m 3 2 ), and represents a classical contribution.
Using (3.1) and (3.2) for the three-and four-point amplitudes, respectively, we obtain a numerator (5.9) Using (5.6) to perform the sum over the internal polarisations, this becomes (5.10) Hence the relevant integral for the classical contribution to the amplitude is (5.12) The propagators D 1 , . . . , D 4 are given in the following table: Using IBP relations [95][96][97][98] we can reduce (5.11) to the master integral G (1) := G This result is valid in D-dimensions, and is in agreement with Eq. (B.16) of [117], which was obtained from the probe approximation of the deflection angle of a massive particle in the background of a Schwarzschild black hole in D dimensions. The master integral G (1) can be evaluated to hence in four dimensions the amplitude becomes

The two-loop HEFT amplitude
In this section we give explicit results for the 2MPI cut diagrams contributing to the two-loop elastic HEFT amplitude, from which we then extract the deflection angle.
6.1 First three-graviton cut: probe limit The first two-loop diagram we consider is the "fan diagram" below. Its integrand is: where the arguments 1, 2, 3 in the numerators are a shorthand notation for ℓ 1 , ℓ 2 , ℓ 3 .
Using the novel BCJ representation (3.4) of the five-point amplitude, we can write the numerator of (6.1) as Performing the sums over internal polarisations, the numerator becomes  .3) is as in the one-loop case, with the replacements done with respect to ε µ ℓ 1 ε ν ℓ 1 , ε µ ℓ 2 ε ν ℓ 2 , ε µ ℓ 3 ε ν ℓ 3 in the numerator. Using IBP reductions, we can re-express the result in terms of a single master integral All other master integrals cancel in the sum over the three nested commutators. One then arrives at the final result for the contribution from this diagram: where y is defined in (2.6). This result is in agreement with Eq. (B.17) of [117]. The value of G (2) is evaluated in Appendix A.3 to In four dimensions, (6.5) reduces to 3 (y 2 − 1) 2 ǫ . (6.7) 6.2 Second three-graviton cut: beyond probe limit Next we consider the two "zig-zag" diagrams, together with the two radiation reaction diagrams. These are given by ∪ three other terms, (6.8) where the ∪ symbol denotes the operation of merging the corresponding cuts from the different cut integrands, producing an integrand which can then be evaluated with IBP reductions. Traditionally the radiation reaction diagrams (second line of (6.8)) are treated separately from the remaining conservative term, but we find it more natural to combine this contribution together with the zig-zag diagram.
The detailed form of the integrands of these four cut diagrams are given separately in Appendix B. After performing the IBP reduction, we can express the D-dimensional where the coefficients c i are given by The integrals G i 1 ,i 2 ,i 3 ,i 4 ,i 5 ,i 6 ,i 7 ,i 8 ,i 9 are defined as with the propagators We use the following box topologies to denote all the propagators that can appear: (6.13) Two comments on two of our master integrals are in order here.
• The master integral G 1,0,1,1,1,0,0,1,1 has the topology of a bow-tie diagram, which is a product of two one-loop triangles, and hence is finite in four dimensions and proportional to (−q 2 ) −2ǫ . It appears in (6.9) with a coefficient proportional to ǫ and hence it does not contribute to four-dimensional classical physics. If one is interested in D > 4, this integral would have to be included.

Four-graviton cut
For completeness, we also compute the four-graviton cut diagram involving a fourpoint graviton amplitude. However, it will turn out that the contribution of this cut is already accounted for by the first two (zig-zag) diagrams in (6.8) as we now explain. The relevant cut diagram is In the final step we have to identify any overlap with the three-graviton cut considered in Section 6.2 in order to avoid double counting, and isolate any new contribution if present. It is easy to see that the terms in the last three lines are already contained in the integrands corresponding to the zig-zag diagrams (6.8). The terms in the first two lines are new contributions that are not detected by the three-graviton cut and can all be reduced to the bow-tie master integral G 1,0,1,1,1,0,0,1,1 : As discussed earlier this master integral is finite in four dimensions and thus does not contribute to classical physics. Hence the four-graviton does not give any new contribution.

Summary
In conclusion, the complete D-dimensional integrand is obtained by adding (6.5) (and the same term with m 1 and m 2 swapped) and (6. As explained earlier, the four-graviton cut in (6.17) does not give any new contribution to the integrand, hence the final result of our evaluation is given by (6.20).
In the next section we evaluate the relevant integrals appearing in that equation, and will then state our (integrated) result at two loops in Section 8, where we also derive the scattering angle up to 3PM.

Generalities
In this section we discuss the evaluation of the two-loop master integrals using the differential equation approach of [95,99,119]. This leads to a canonical basis for the master integrals encountered in Section 6.2 which we will evaluate in the full soft region. This method was applied to the specific type of integrals which appear in the post-Minkowskian expansion of classical scattering, initially for the potential region in [17] and later for the soft region in [19,21,118]. The following discussion is valid in an arbitrary number of dimensions, while for the evaluation we will work in 4 − 2ǫ dimensions except for the probe limit contributions, which we will evaluate in general dimensions.
We first set up the system of differential equations in terms of a basis of master integrals: To get a canonical basis, we need to perform a linear transformation on the integral basis. In general, if we transform the basis as I ′ j = SI j , the differential equation becomes The canonical basis is defined in such a way that the transformation matrix is proportional to the dimensional regularisation parameter ǫ. This can be realised in two steps: 1. Diagonal transformation S 0 : Since the matrix A has poles at y = 0, −1, +1, we choose the following ansatz for the diagonal transformation, where the parameters b i do not depend on y. The matrix A is then transformed as follows: We can choose the parameters in S 0 such that the matrix A is transformed into for general dimensions D = d 0 − 2ǫ, with A 0 and A 1 being independent of ǫ.
2. Non-diagonal transformation S 1 : Here S 1 is independent of ǫ and satisfies the following equation: Then it is easy to see that the new integral basis is a canonical basis, which satisfies The formal solution of this equation is

Canonical integral basis
We only need to consider the following seven master integrals: where circles on lines denote squared propagators and for I 6 we have also indicated the presence of two numerator factors. We also set the master integral G 1,0,1,1,1,0,0,1,1 to zero as it does not contribute to classical four-dimensional physics.
We set up the system of differential equations as Following the two steps outlined in Section 7.1, we arrive at The canonical basis satisfies the following differential equation, where At two loops, the study of the boundary conditions is well documented in the literature [17,19,118], albeit for different bases and using slightly different methods. As discussed in Section 6.2 the bow-tie topology would have to be included in the system of differential equations if we want to work in general dimensions. In this paper we focus on D = 4 − 2ǫ, and hence we can drop the this topology.
The final step is now to provide appropriate boundary conditions for the solutions to this system of differential equations. A systematic method is provided in Appendix A, based on the study of the asymptotic expansion of Feynman integrals [100,101] in the static limit y → 1. Using the results of Appendix A, we find With these boundary conditions, we finally obtain the solutions of the canonical basis from the differential equations 18) and γ denotes the Euler-Mascheroni constant.
The possible boundary conditions for the master integrals are related to different physical regions (e.g. the potential or the soft region), see for instance [100,120,121] and Appendix A. These regions are related to different asymptotic behaviours of the integrals in a given limit (e.g. the static limit y → 1 or the ultrarelativistic limit y → ∞). In general, a master integral has several regions corresponding to a given limit. For a set of master integrals of a system of differential equations, the asymptotic behaviour is also constrained by the system, and therefore the regions for the master integrals are usually not independent. Here we related the independent regions for the master integrals to the physical regions. After applying the differential equation constraints, we have two independent regions characterised by the boundary values 8 We find that I (R 1 ) 3,LD corresponds to the potential region and I (R 2 ) 2,LD to the radiation region. We will also see in the next section that the real part of the two regions generates the conservative and radiation reaction contributions to the scattering angle, respectively. An interesting observation is that the values of the leading terms in the two regions are all determined from three master integrals I ′ 2 , I ′ 3 and I ′ 4 , where in the static limit the latter can be expressed in terms of the former two.

Final result and scattering angle from the HEFT amplitude
We now combine the results of Sections 6 and 7 and provide an expression for the scattering angle up to two loops.

(8.3)
The terms in the HEFT amplitude of order O(m 4 1 m 2 2 ) and O(m 4 1 m 2 2 ), which contribute only to the real part of the amplitude, have been computed in Section 6.1, see (6.7).
Combining them with (8.2), we arrive at the complete 2MPI two-loop HEFT amplitude: The second line is the radiation reaction part. We observe that (8.4) is identical to the matrix element of the operator N introduced in [94], see Eq. (3.30) of that paper.
Finally, we compute the scattering angle. First we introduce the Fourier transform to impact parameter space: At tree level we find where we used J = P b, with P being defined in (2.9).
At one loop, by performing the Fourier transform we find δ (1) Finally, at two loops our result is Re δ , (8.8) where Re M (2,2MPI) is given in (8.2). The scattering angle is then obtained from (4.1), and our final result up to two loops, or 3PM, is This is in agreement with the result of [15,18,19,102]. The radiation reaction part (third line of (8.9)) was also computed in [18,103,104] and our result agrees with theirs.
We conclude this section with two comments. We note that the imaginary part of our 2MPI amplitude, quoted earlier in (8.3), is very similar (but not identical) to that of [18,19]. There is no reason for these quantities to agree completely since there are subtle differences in how iterating terms are removed in [18,19]. However the relation [19] lim ǫ→0 Re δ (2) HEFT rr = − lim ǫ→0 πǫ Im δ (2) HEFT , (8.10) where rr stands for the radiation reaction part, does hold also in our case, with our Re δ HEFT rr being identical to that of [18,19]. The appearance of the imaginary part (8.3) in the 3PM amplitude is a very interesting fact, related to radiation emission. The radial action (from which the deflection angle is extracted) is believed to be related to the matrix element of the hermitian operator N introduced in [94]. As discussed in Section 4.4, the real part of our δ HEFT is identical to this matrix element -with the extra imaginary part only arising from the first diagram in (2.17) of that reference (the triple cut involving two massive particles and a radiation graviton). Hence at 3PM order we just need to drop the imaginary part of δ HEFT to obtain the matrix element of N and hence the radial action.
Beyond 3PM the situation is more involved. It is reasonable to expect that the correct prescription in our HEFT approach is again to relate δ HEFT to a specific matrix element of N, generalising (2.10) of [94] to the next PM order. Our expectation is that our approach again automatically implements many of the subtractions implied by the definition of N but will require the subtraction of a larger set of radiation cut diagrams. 9 All-loop conjecture for the classical contribution in the probe limit In this section we present a conjecture for the all-loop structure of the leading terms in the 2MPI amplitude in the probe limit. Consider the following "fan" diagram at L loops: The lower part of the diagram is tree-level amplitude in HEFT, which is of O(m 2 1 ) while the upper is of O(m L+2 2 ) at L loops. Therefore the diagram behaves as O(m 2 1 m L+2 2 ), which dominates the classical gravitational scattering in the probe limit m 2 ≫ m 1 .
A similar procedure to the one followed in previous sections at one and two loops allows us to express the L-loop diagram in terms of a basis of master integrals. Leaving a systematic study of general higher-loop diagrams for the future, here we conjecture that the simple fan structure found at one and two loops gives the leading contribution in the probe limit also at higher loops. In particular, all the internal propagators arising from the lower tree-level HEFT amplitude are absent in the master integrals. That implies that there is only one master integral, given by The needed HEFT amplitude has some simple general properties, which allow us to write down a general expression for the fan diagram at an arbitrary number of loops L.
The tree-level amplitude in the HEFT as obtained from the double copy contains a factor 1 (p·v) 2 from each heavy-mass propagator, which will generate a factor 1 y 2 −1 after the IBP reduction, resulting at L loops in a factor of 1 (y 2 −1) L . The terms that do not contain any heavy propagators will have a similar velocity dependence as the tree-level amplitude, which is of O(y 2 ) in the limit y → ∞. Therefore we conclude that the L-loop contribution to the classical potential in the probe limit takes the form We can now compare the tree-level amplitude together with the one-loop and two-loop results in (5.13) and (6.5), respectively. We then observe some interesting pattern for the D-dependent factors: • The overall prefactor is of the form • In the degree-(2L+2) polynomial in y, the coefficient of the y 2 term is proportional to D − 2 + L(D − 3).
• The ratio between the y 2j+2 coefficient and the y 2j coefficient in the polynomial is proportional to D − 2 + L(D − 3) + 2j.
Combining these observations, we expect that the amplitude contributing to the scattering angle in the probe limit at L loops takes the form where c L,i are some undetermined D-independent constants. For example, at three loops, the relevant amplitude in the probe limit is conjectured to be We have verified this pattern up to 9PM by comparing (9.5) to the scattering angle for a massive probe particle in the background of a Schwarzschild black hole (see for example the Appendix of [117], and we conjecture that it holds for any loop order.

Outlook
In this paper we have studied the 2 → 2 scattering amplitude of two heavy particles at one and two loops using a formalism based on heavy-mass effective field theory. We have defined a new HEFT phase as the Fourier transform to impact parameter space of the sum of all 2MPI diagrams in this theory. We have found that this phase is directly related to the scattering angle by taking the derivative with respect to the angular momentum of the binary system.
An important question is to understand more precisely the connection of our work to the results of [24,94], and in particular the relation of our HEFT phase to the radial action (or the real part thereof), which our phase agrees with for the conservative part. This hints at a strong relation between the ab initio removal of iterating terms advocated in [24] and our approach, avoiding the task of explicitly evaluating terms which exponentiate as usually done in the eikonal approach. It would also be interesting to compare to worldline approaches, which tackle directly the computation of an effective action sitting in an exponent [36].
An obvious task consists in generalising our work to higher loops. The use of our novel colour-kinematic duality [85] will be crucial to produce compact expressions for the tree amplitudes entering the unitarity cuts. An additional challenge beyond 3PM is that elliptic functions appear [24,39]. The inclusion of spin effects would also be interesting to discuss within our HEFT approach, as well as the computation of other classical observables in the spirit of [1,122].
Finally, since our approach has produced integrands valid in D dimensions, it would be interesting to work out the 3PM integrated amplitude beyond the probe limit in a general number of dimensions. This requires work on the relevant master integrals, setting up the relevant differential equations around integer dimension D > 4, and computing the corresponding boundary conditions. We leave these questions for future work.

A.1 The general methodology
In order to find the asymptotic behaviour [100] of the master integrals, it is convenient to use Feynman parameters. A general scalar integral with a single dimensionless kinematic parameter x (e.g. a ratio of masses and Mandelstam variables) takes the following form where a = m i=1 a i , each a i is the power of a Feynman propagator, and F , U are the Symanzik polynomials. The numerical evaluations of the integrals are usually performed using sector decomposition [123,124] or Mellin-Barnes representations [125].
Here we use the Cheng-Wu theorem [126] to remove the delta function, as follows: and the hat denotes the omission of the corresponding integration. In general, the integral can be expanded as 9 where the summation is over different regions R, and LD stands for leading term. The parameters w R depend on the number of dimensions and the difference of any two w R is never an integer. In the literature, each non-vanishing I LD is called a region, while the dots in (A.3) denote higher-order terms in each region.
In principle, all the expansion coefficients can be written as an integral over the Feynman parameters. However, the closed integral form of I (R) LD for a general Feynman integral is still not known. It was conjectured that all these integral forms can be obtained by re-parameterising the Feynman parameters α i → α ′ i (α 1 , · · · , α m , x) and then extracting the common factor x w R as Then for the region R we have In practice, the method introduced in [101,127] is more efficient, and we now review it. The key assumption in [101] is that the leading terms in the power expansion of the integral in the parameter x can be obtained by rescaling the Feynman parameters by some powers of x and expanding the integrand directly. 10 In particular, the monomials in the U and F functions have the form x r 0 m i=1 α r i i . We then rescale α i as α i x v i , and denote by U ′ and F ′ the rescaled functions, dropping all subleading terms in x. Furthermore, we introduce vectors containing as components the powers of the monomials in the U ′ and F ′ functions where t and s denote the number of monomials in U ′ and F ′ . More concretely, each r is a vector whose components are the degrees of the monomials with respect to the variables x, α 1 , α 2 , · · · , α i , · · · , α m , where the hat indicate omission of the corresponding variable. We then scan all possible leading terms (in x) by choosing different monomials in the original U and F functions. For each choice of leading term, the new U ′ and F ′ are represented as sums over the following monomials from the original U and F functions: As the terms in the U ′ or F ′ functions are of the same degree in x after the rescaling, we have the constraint equations for the rescaling degrees where v is the vector {1, v 1 , · · · ,v i , · · · , v m } made of the powers of the rescalings by x, i.e. the rescaled variables are If there is no solution to the constraint equations, then such a region does not exist. Indeed, if the constraint equations do not fix the values of v, then there exists some rescaling freedom such that the leading integral I LD has the scaling where w is non-zero. In this case the leading integral always vanishes in dimensional regularisation. If the constraint equations fix the scaling vector v, then using the solution for v from (A.9) it is easy to check that all the other monomials in the original U and F functions are of higher degree in x after rescaling.
We scan over all the U ′ or F ′ functions, and solve the constraint equations (A.9) for the monomials in the U ′ or F ′ function. Finally we are able to find all the regions: where, again, the difference of any two w (R) is never an integer.

A.2 One-loop example
As an example, we consider the one-loop box integral . (A.13) After using Feynman parametrisation, we arrive at (A.14) where y ′ = p 1 ·v 2 m 1 and x = − p 2 ·p 3 Then we scan over all the possible leading terms for the U or F function. For example if we choose the leading terms to be 1, α 3 and α 2 1 , α 2 2 , α 2 α 1 , then by (A.9) the constraint equations for the rescaling degrees are These equations do not determine the rescaling degrees and hence the region with these leading terms does not exist.
If we choose the leading terms to be 1, α 3 and α 2 1 , α 2 2 , α 2 α 1 , xα 3 for the U and F functions, respectively, then the constraint equations are In this case the solution exists, which is given by It is easy to check that the chosen monomials are the leading terms under this rescaling.
After scanning over all the possible leading terms in the U and F functions, we find that the only nontrivial solution for the rescaling degrees is (A.18). Therefore we get, for the asymptotic behaviour, There is only one region for the box integral I

A.3 Two-loop examples
In this section we discuss the asymptotic behaviour of the two-loop master integrals.
In order to evaluate our two-loop HEFT amplitude, only the results quoted in this Section are required. In particular, our two-loop master integrals contain both linear propagators raised to arbitrary powers and delta functions. These are implemented by the following replacements: so that in Feynman parameterisation, the parameter corresponding to a cut propagator is integrated over the domain (−∞, ∞). Such integrals can usually be transformed into contour integrals and evaluated as a sum of residues.
According to the canonical basis in (7.16) and the differential equation in (7.14), the regions I