The Two-Loop Four-Graviton Scattering Amplitudes

We present the analytic form of the two-loop four-graviton scattering amplitudes in Einstein gravity. To remove ultraviolet divergences we include counterterms quadratic and cubic in the Riemann curvature tensor. The two-loop numerical unitarity approach is used to deal with the challenging momentum dependence of the interactions. We exploit the algebraic properties of the integrand of the amplitude in order to map it to a minimal basis of Feynman integrals. Analytic expressions are obtained from numerical evaluations of the amplitude. Finally, we show that four-graviton scattering observables depend on fewer couplings than naively expected.

Scattering amplitudes are ubiquitous in high-energy physics: they connect physical observables and the quantum field theories describing the different forces of Nature. By understanding the structure of amplitudes, we can learn about properties of these theories and their physical implications. Unlike other field theories, such as Yang-Mills', Einstein's theory of general relativity cannot be consistently quantized in its minimal form. Indeed, it was shown over 30 years ago [1][2][3][4] that quantum effects render scattering amplitudes ill defined in the ultraviolet (UV). Since then, our understanding of the UV properties has been refined by the study of scattering amplitudes in this regime, both in Einstein gravity [5][6][7] and in supersymmetric extensions of it such as maximal supergravity [8,9]. New results for amplitudes have also been obtained, but mostly in supersymmetric theories [10][11][12][13][14][15][16][17]. Computations in Einstein gravity are famously involved, and while the one-loop four-graviton amplitudes have been known for decades [11], the two-loop amplitudes remained unknown till now. In this letter, we present them for the first time.
Following the detection of gravitational waves, interest in quantum gravity amplitudes has surged as a means to predict the classical gravitational dynamics of large massive objects in the post-Minkowskian approximation, most notably that of black-hole binaries [18][19][20][21][22][23][24]. Already some time ago, the two-loop scattering amplitudes in string theory were understood to yield the classical scattering angle of massless particles [25], but the validity of this observation was recently questioned [26]. Our results give new insights on the theoretical properties of Einsteins theory of gravity and associated physical phenomena. The amplitudes presented here were already used in ref. [27] to confirm the results of ref. [25].
Our calculation is performed with techniques developed for the computation of amplitudes in the Standard Model of particle physics. They have already been successfully applied to computations of planar scattering amplitudes in QCD, both numerically [28][29][30] and analytically [31,32], and are well suited to address the challenges of a quantum gravity calculation. We use a variant of the unitarity method [33][34][35] suitable for numerical computations, the two-loop numerical unitarity approach [28,36,37], which replaces Feynman-diagram input with numerical evaluations of on-shell tree amplitudes. It bypasses the explicit construction of the integrand of the amplitude, and directly reduces it to a minimal basis of Feynman integrals with unitarity-compatible integration-by-parts relations [38,39]. Analytic expressions, provided in a set of ancillary files, are reconstructed from exact numerical evaluations of the amplitudes.
Four-Graviton Scattering Amplitudes. We consider four-graviton scattering in Einstein gravity. The theory is not renormalizable [1][2][3][4], and we work in the effective field theory proposed in ref. [40]. The Lagrangian L is where we suppress terms not relevant for our two-loop calculation such as higher-order operators and those proportional to the equations of motion [1][2][3]. It is given in terms of the Einstein-Hilbert (EH) Lagrangian L EH , complemented by the Gauss-Bonnet (GB) and the R 3 counterterms [2,3,41,42], denoted L GB and L R 3 respectively, whose role is to cancel the UV divergences inherent to L EH . The different contributions to L are where g = det(g µν ) with g µν the metric tensor, R the Ricci scalar, R µν the Ricci tensor and R µνρσ the Riemann curvature tensor. We work in the 't Hooft-Veltman (HV) scheme of dimensional regularization, with D = 4 − 2 . So that each contribution has the same dimensions, we introduce the dimensionful quantity µ which includes conventional factors in dimensional regularization, µ 2 = (4π) −1 e γE µ 2 0 . The coupling κ is related to Newton's constant G N , κµ − 0 = √ 32πG N . The divergent parts of the bare couplings have been determined previously [2,3,5]. The renormalized couplings c GB (µ) and c R 3 (µ) will be discussed at the end of this letter. We compute graviton scattering on a flat background η µν , parametrized by the linear split g µν = η µν + κh µν , where h µν is the graviton field [43]. Perturbation theory is defined as an expansion in powers of κ. The main results of this letter are the helicity amplitudes for fourgraviton scattering M h (s, t; ) at order κ 6 , as a function of s = (p 1 +p 2 ) 2 and t = (p 2 +p 3 ) 2 , with outgoing momenta p i . We will often suppress dependence on Mandelstam variables. The helicity assignments are specified by h = {h 1 , h 2 , h 3 , h 4 }. It is sufficient to compute amplitudes with h = {±, +, +, +} and {−, −, +, +} since all others are related by symmetry. We define the perturbative expansion of the helicity amplitudes through withκ = κµ − /(4π) and helicity-dependent phases N h given in footnote [44]. That is, we normalize M h so that the coefficients M are Lorentz invariant. The index j in eq. (4) is in one-to-one correspondence with the looporder of the contributing diagrams for L EH . This correspondence breaks down for L GB and L R 3 as can be seen by the power of the coupling in the three-point vertices of each term in eq.
has tree, one-loop, and two-loop contributions, depending on which vertex appears. Schematically, where we include sample diagrams for each contribution. White blobs denote L GB vertices, grey blobs denote L R 3 vertices and L EH vertices have no decoration. The first non-vanishing contributions from L GB and L R 3 appear at O(κ 6 ). The amplitudes M (j) h computed from L in eq. (1) are UV finite, but there remain infrared (IR) singularities [45][46][47][48]. It is known [49] that there are no collinear singularities, and the soft singularities exponentiate. We define and construct finite functions F Comparing eqs. (4) and (6), we can write the F This object captures the new four-dimensional information at two-loops.
Computation. The main obstacles in computing the amplitudes M (2) h are rooted in the involved Feynman rules derived from L in eq. (1). Vertices have many terms with high powers of the momenta, making it hard to construct the integrand of the amplitude. Furthermore, despite the simple kinematics of the process, the reduction to a set of master integrals is challenging because the integrand has high powers of the loop momentum.
The framework of two-loop numerical unitarity [28,36,37] is particularly well suited to address these challenges. The starting point is the following parametrization of the integrand of an amplitude [50], denoted M (k) ( l ), with M Γ a set of master integrands, S Γ a set of surface terms, P Γ the set of propagators ρ j associated with each propagator structure Γ, and l the set of loop momenta. The set ∆ of relevant propagator structures is characterized in fig. 1. The decomposition (8) is an ansatz for the integrand of the amplitude. The undetermined coefficients c Γ,i are constrained from the factorization properties of the integrand in loop-momenta configurations Γ l where the propagators in P Γ vanish: with T Γ the tree amplitudes corresponding to the vertices in Γ. The sum over states runs over D s -dimensional graviton helicity states, and the sum over Γ runs over propagator structures such that P Γ ⊆ P Γ . The system of eqs. (9) is constructed numerically. Assuming we have built the ansatz (8) and can evaluate the product of trees in (9), this reduces the calculation of the amplitudes at a phase-space point to solving the linear system of eqs. (9). Indeed, once all c Γ,i have been determined, we directly obtain the decomposition of the amplitude in terms of master integrals, where the integrals I Γ,i correspond to the master integrands in M Γ . In the following, we discuss the construction of ansatz (8) and the computation of tree-amplitudes for eqs. (9). We first focus on the evaluation of tree amplitudes. We use a fast numerical algorithm provided by Berends-Giele recursion [51]. In the pure Einstein-Hilbert theory, L EH , we use the reformulation in terms of cubic interactions proposed in ref. [52]. For counterterm contributions, vertices are computed using the program xAct [53][54][55]. Our Berends-Giele recursion allows for EH, GB and R 3 tree amplitudes. We use integer values for the state counting parameter D s that are large enough to recover the full momentum dependence, i.e. D s ≥ 6.
Next, we discuss the construction of the ansatz (8). It depends on the power-counting properties of the theory and the kinematics of the process. First, we build the full set of propagator structures ∆, which contains both planar and non-planar contributions, see fig. 1. For each Γ ∈ ∆ we then construct the function space M Γ ∪ S Γ . The elements of the space, m Γ,i ( l ), are polynomials in the components of the loop momenta l . The linear span of the space is controlled by the theory-specific maximal polynomial degree. In Einstein gravity one naively expects that the polynomial degree required is twice that of Yang-Mills. The next step is the construction of the surface terms in S Γ , which integrate to zero. A subset of these can be built from tensor reduction techniques [28]. The rest are constructed from integration-by-parts (IBP) provided that so that no new higher propagator powers are generated in the procedure [38,39]. The f j are polynomials in loopmomenta components, and no summation over the index j is implied. Solutions u ν i to eq. (12) are power-counting independent and referred to as IBP-generating vectors. For each Γ, once a set of vectors is found, surface terms are constructed as follows. Consider a polynomial t r ( l ) in the loop-momenta components and a solution u ν i,s to eq. (12). We then insert t r ( l )u ν i,s in eq. (11) to obtain the surface term where D-dependence may arise from the divergence term. We complete the IBP-generating vectors obtained in ref. [28] for planar topologies with the ones for non-planar topologies. To obtain surface terms with the suitable power-counting, we must use a sufficient set of polynomials t r ( ). Each vector u ν i appears in many surface terms, offering the opportunity for caching in the numerical approach. The set of master integrands M Γ in eq. (8) is the complement of S Γ in the integrand function space.
We are now ready to construct the system of eqs. (9). For each numerical phase-space point, choice of = (4 − D)/2, and value of D s , we can solve for the coefficients c Γ,i , yielding the decomposition (10) in terms of master integrals. To expand the result in , we first reconstruct the dependence of the coefficients on this parameter and D s . They are rational in , and so we compute a sufficient number of samples to apply Thiele's formula [56]. The coefficients of in the numerator of this rational function depend on D s . In pure gravity, they are quartic polynomials in D s . GB counterterm amplitudes have rational D s dependence, with numerators that are cubic in D s and denominators that are simply D s − 2 [57]. The R 3 counterterm amplitudes are D s independent. We determine the D s dependence from enough numerical samples.
Through this procedure, we obtain master integral coefficients as rational functions in , with analytic D s dependence at numerical values of s and t. We set D s = 4 − 2 , as prescribed by the HV scheme, insert the expressions for the master integrals [58][59][60][61], and expand the result in . With modern mathematical tools [62] we can express the amplitudes in a basis B of classical polylogarithms, whose elements are denoted h i ∈ B.
Using one-loop amplitudes we computed within the same framework, we obtain the remainders in (7) at the chosen phase-space point as a linear combination of the h i : Finally, we can reconstruct the full analytic result from a sufficient number of numerical samples. As noted e.g. in refs. [31,32,63,64], it is more efficient to reconstruct the coefficients d i (s, t) of eq. (14). The coefficients d i are rational functions of x = t/s, and the s dependence can be reconstructed from dimensional analysis. Therefore, we can use the univariate Thiele formula to reconstruct the rational functions d i . This process requires around 20 numerical samples for each helicity. Numerical stability issues are sidestepped by employing finite-field arithmetic [65,66]. Combining the results from evaluations over two different finite fields with cardinality of order 2 31 , we lift the results to the field of rational numbers using the Chinese remainder theorem and rational reconstruction techniques [67].
Results. We have computed the four-graviton amplitude for the three independent helicity configurations h = {−, −, +, +}, {−, +, +, +} and {+, +, +, +}, up to order O(κ 6 ) in the effective field theory of eq. (1). The amplitudes are obtained by computing the remainders of eq. (7) and then reinstating the IR singularities. By taking into account the contributions from the GB (up to one-loop) and the tree-level R 3 counterterms, we also obtain the two-loop amplitudes in the EH theory. We note that the evaluation of the remainders requires one-loop amplitudes through O( ). All scattering amplitudes, in the HV scheme, are provided in ancillary files.
We performed several checks on our results. First we verified that all the poles in our amplitudes, which are by construction of IR origin, are accounted for by the universal structure (5). Second, some parts of the different ingredients we require to compute the O(κ 6 ) amplitudes have been obtained previously, giving completely independent checks. One-loop amplitudes in Einstein gravity were computed in ref. [11] through O( 0 ). We confirm the {±, +, +, +} results, the {−, −, +, +} amplitude up to a sign [68] and agree with an independent computation [69]. The counterterm amplitudes were partially known. We reproduce the divergent pieces of the counterterm amplitudes for {+, +, +, +} given in ref. [5]. The GB tree-level and one-loop amplitudes match an independent computation of the {+, +, +, +} and {−, −, +, +} helicities [69]. Regarding R 3 , we reproduce known results for the tree-level amplitudes with a single R 3 insertion [5,70]. Third, we could check some of the O(κ 6 ) amplitudes: the {+, +, +, +} amplitude matches the results of ref. [69] and is consistent with ref. [70], and our results for the {−, −, +, +} amplitude match the behaviour established in [71]. Finally, our amplitudes behave consistently with factorization in the limits where the Mandelstam invariants s, t or u = −s − t vanish.
While the results are too large to print in this letter, we can quote the result for the {−, −, +, +} remainder in the s-channel Regge limit. Defined by s −t > 0, this limit is directly relevant for linking scattering amplitudes to classical dynamics [27]. With our choice of IR subtraction, we find where we introduced L = log(−s/t). This expression is independent of the scale µ introduced in eq. (2), consistent with the following discussion. Finally, we consider the remainders' dependence on the couplings in eq. (3). R (2) {−,−,+,+} is independent of both c GB (µ) and c R 3 (µ), while R (2) {±,+,+,+} depend on the unique combination This observation is tightly connected with the dependence of the remainders on the scale µ introduced in eq. (2). Indeed we find that where M tree,R 3 h is the tree amplitude with a single R 3 insertion, which vanishes for h = {−, −, +, +}. This extends the scale dependence proposed for h = {+, +, +, +} in refs. [5][6][7] to all helicities. The scale dependence in eq. (17) takes a much simpler numerical form than the divergent parts of the couplings in eq. (3). Requiring that the remainders are independent of µ allows to determine the µ-dependence of the coupling c(µ) [72]. This is sufficient for the remainders to be well defined, and it is a weaker condition than requiring that C GB and C R 3 be µ-independent. The fact that remainders display a reduced dependence on the couplings in eq. (3) is interesting for two reasons. First, the same is not true regarding how the couplings contribute to the cancellation of the UV poles. This yields two independent equations, allowing to uniquely fix the divergent part of the couplings C R 3 and C GB . Second, this implies that physical observables related to four-graviton scattering at two-loops depend on fewer parameters than those appearing in the effective field theory. It is likely that this degeneracy is a consequence of the evanescence of the GB counterterm, which would then imply that our observation should extend to twoloop amplitudes of higher multiplicities. This is consistent with the results of ref. [70], which can be shown to imply that the two-loop five-point all-plus amplitude depends on the same combination c(µ) of couplings.
Conclusions. In this letter, we presented the O(κ 6 ) four-graviton amplitudes in Einstein gravity, including contributions from counterterms. The computation of graviton amplitudes is notoriously difficult but our results show that modern field-theory methods, notably the numerical unitarity approach, are able to tackle these challenges. Our results give new insights into the analytic structure of the theory, contributing [27] to the ongoing effort to bridge multi-loop scattering amplitudes and classical gravitational dynamics. We find that the {−, −, +, +} remainder only depends on the coupling κ, while the {±, +, +, +} amplitudes depend on a single additional coupling. This implies that observables constructed from these remainders only depend on two out of the three couplings appearing in the effective field theory.
Multiple future directions are worth pursuing. Given the mild dependence of our approach on the number of scales, a clear next step is to consider amplitudes including massive particles. Another natural extension is towards higher loop corrections. Both will be of direct relevance for exploring the classical gravitational dynamics of large massive objects. The analytic results we present also provide insights into the analytic properties of the amplitudes, stimulating the development of more efficient techniques to tackle calculations at higher loop orders and multiplicities.