Conservative Black Hole Scattering at Fifth Post-Minkowskian and First Self-Force Order

We compute the 5PM order contributions to the scattering angle and impulse of classical black hole scattering in the conservative sector at first self-force order (1SF) using the worldline quantum field theory formalism. This challenging four-loop computation required the use of advanced integration-by-parts and differential equation technology implemented on high-performance computing systems. Use of partial fraction identities allowed us to render the complete integrand in a fully planar form. The resulting function space is simpler than expected: in the scattering angle we see only multiple polylogarithms up to weight three, and a total absence of the elliptic integrals that appeared at 4PM order. All checks on our result, both internal - cancellation of dimensional regularization poles, preservation of the on-shell condition - and external - matching the slow-velocity limit with the post-Newtonian (PN) literature up to 5PN order and matching the tail terms to the 4PM loss of energy - are passed.

Binary black hole (BH) and neutron star (NS) mergers are today routinely observed by the LIGO-Virgo-KAGRA gravitational wave detectors [1][2][3].With the advent of the third generation of gravitational wave detectors [4][5][6], and LISA's recent approval by the European Space Agency, we anticipate an experimental accuracy increase that will enable unprecedented insights into gravitational, astrophysical, nuclear, and fundamental physics.From these experimental programs emerges the theoretical imperative to reach utmost precision in the gravitational waveforms emitted by these violent cosmic events.To meet this demand, a synergy of perturbative analytical and numerical approaches is needed to solve the classical general relativistic two-body problem.The former encompasses the post-Newtonian (PN) [7][8][9] (weak gravitational fields and non-relativistic velocities) and post-Minkowskian (PM) [10][11][12][13][14] (weak fields) expansions; the latter encompasses modern numerical relativity [15][16][17].Gravitational self-force (SF) (small mass ratio) [18][19][20][21], meanwhile, is a hybrid approach: the perturbative SF equations typically being solved numerically.On the analytical side, the incorporation of perturbative quantum field theory (QFT) techniques has significantly strengthened this program, most recently within the PM expansion.
The SF expansion [18][19][20][21] is a complementary perturbative scheme in which one assumes a hierarchy in the two BH or NS masses, m 1 ≪ m 2 , but works exactly in G.The self-force expansion therefore extends systematically beyond the geodesic motion of a probe mass moving in the background of a heavy BH or NS.One may overlay the PM loop expansion with the SF expansion: the PM problem factorizes into separate gauge-invariant SF sectors that may be targeted individually.Concretely, for the 5PM four-loop problem one finds 0SF (known), 1SF (computed here) and 2SF contributions.The complexity of the Feynman integrals to be performed grows considerably with the SF order.Moreover, overlaying the PM with the SF expansion for the scattering sce-

FIG. 1:
The six top-level sectors of the four-loop planar integral family (5), yielding the m 2 1 m 4 2 5PM-1SF contributions.The δ − (ℓi • ui) can here be interpreted as cut propagators -dotted lines, which in the WQFT context alternatively denote the background worldlines.In this planar four-loop family we have 13 active propagators (7), the active graviton propagators that may become radiative (10) being depicted in red.
In this Letter, we compute the previously unknown 5PM contribution at first order in self-force.Our computation lies at the frontier of present Feynman integration technology.In order to master it, we optimized on all aspects of this high-precision challenge: the integrand was produced using the worldline quantum field theory (WQFT) formalism [14,31,43,46], with partial fraction identities used to perform a "planarization" prior to integration.The integration-by-parts (IBP) reduction employed an improved version of Kira [101,102].
Worldline Quantum Field Theory.-The nonspinning BHs or NSs are modeled as point particles moving on trajectories x µ i (τ ).In proper time gauge ẋ2 i = 1 the action takes the simple form suppressing a gauge-fixing term S gf .We employ a nonlinearly extended de Donder gauge that maximally simplifies the three-and four-graviton vertices, and use dimensional regularization with D = 4 − 2ǫ in the bulk.Both the worldline and gravitational fields are expanded about their Minkowskian (G 0 ) background configurations yielding the propagating worldline deflections z µ i (τ ) and graviton field h µν (x).The incoming data is then spanned by the impact parameter b The quest of solving the equations of motions of Eq. ( 1) in a G-expansion is solved upon quantizing the perturbations z µ i and h µν : the tree-level one-point functions then solve the classical equations of motion [103].The impulse of (say) the first BH, ∆p µ 1 , then emerges from ∆p µ 1 = lim ω→0 ω 2 z µ 1 (ω) , working in momentum (energy) space.The WQFT vertices are given by standard bulk graviton vertices -at 5PM we require the 3, 4, 5 and 6 graviton vertices -and worldline vertices coupling a single graviton to (0,. . ., 5)-worldline deflections [43,53].We access the conservative sector by employing Feynman propagators (in-out) in the bulk and retarded on the worldline (in-in) [46,104], taking the part real and even in velocity v. Non-trivial Feynman loop integrals emerge in WQFT due to the hybrid nature of the theory: the worldlines only conserve the total inflowing energy, as opposed to full four-momentum conservation in the bulk.The (non-spinning) nth PM contribution to the impulse is given by (n − 1)-loop integrals, plus a trivial Fourier transform over the momentum transfer q.
Self-force expansion.-The 5PM contribution to the complete impulse, ∆p µ 1 = ∞ n=1 G n ∆p (n) µ , factorizes into (effectively) three self-force (SF) contributions: , each of which is separately gauge-invariant.In fact, the SF order may be directly read off a WQFT diagram: the power of m i is given by the number of times the i'th worldline is "touched" -see e.g.Fig. 1, which contains integral graphs belonging to the 1SF (m 2 1 m 4 2 ) sector.Simplest to compute are the probe limit results ∆p (5)µ 0SF and ∆p (5)µ 0SF , which describe geodesic motion and are known to all orders in G [105].They encode the m 1 ≪ m 2 and m 1 ≫ m 2 limits respectively, and are related to each other by symmetry.For the conservative dynamics that we focus on here, the leading (1SF) self-force corrections ∆p 1SF will be a central result of this Letter.Integrand generation.-The 5PM integrand is generated with a Berends-Giele type recursion relation employing the automated vertex rules from the action (1), as discussed in Ref. [53].It is not a bottleneck of the computation.In the 1SF sector this yields a total of 363 WQFT diagrams; the probe limit (0SF sector), which we also generate as a test-bed, is comprised of 63 diagrams.After inserting the Feynman rules using FORM [106], ∆p (5)µ may be reduced to a sum of scalar-type integrals by replacing any loop momenta with a free index as [44] FIG. 2: Two examples of nonplanar loop integrals.By applying the partial-fraction identity (8), we may re-express then in terms of the integrals in Figs.1e and 1f respectively, and thus include them in the planar loop integral family (5).

The dual velocities vµ
The momentum impulse is then expressed as linear combinations of scalar integrals depending trivially on the momentum transfer |q| := −q 2 (this being the sole dimensionful quantity in the problem) and non-trivially In anticipation of the subsequent IBP reduction step, we organize the resulting scalar integrals into families.We introduce the following generic 1SF planar integral family, valid at any L-loop order: where {σ} and {n} denote collections of i0 + signs and integer powers of propagators respectively.The worldline propagators D i (σ i ) are and the massless bulk propagators (gravitons) D IJ with I = (0, i, q) are (suppressing a Feynman i0 + prescription) In total we have L linear and L(L + 3)/2 quadratic propagators at L-loop order.We also allow for derivatives of the one-dimensional delta function δ The four-loop family is illustrated in Fig. 1, with the following diagrammatic rules: The optional arrow in Eq. (7b) denotes causality flow.By interpreting the background worldlines as cut propagators, we "close" the loops of the tree-level WQFT diagrams, and may thus import the notion of planarity from regular QFT Feynman diagrams.We note in passing that this matches the velocity cuts of Refs.[64,66].
Remarkably, the entire 5PM-1SF result for the momentum impulse may be expressed in terms of integrals belonging to this planar integral family alone.To achieve such a representation, graphs with a nonplanar structure -such as the two depicted in Fig. 2 -are systematically eliminated in favor of planar ones.This is done using partial fraction identities on the worldline propagators: where ℓ µ 12 = ℓ µ 1 +ℓ µ 2 , and each linear propagator carries an implicit +i0 + prescription.This identity, which may be applied internally within a multi-loop integral containing linearized propagators, has the effect of "untangling" the crossed bulk propagators, and can be applied repeatedly in order to produce a fully planar integrand.
IBP reduction.-The planar integral family (5) splits into two branches: even (b-type) and odd (v-type) under the operation v µ i → −v µ i .These two branches are thus distinguished by the number of worldline propagators: even (b-type) or odd (v-type), including also the number of derivatives on the delta functions.They may be IBP reduced separately, and in the final answer for the impulse they contribute in the directions of b µ and v µ i respectively (hence the name).
Crucially, all γ-dependence in the integrals lies in the linear propagators and delta functions.The analytic complexity therefore depends highly on the the combination of contractions with v µ 1 or v µ 2 in these propagators.At mSF and nPM order we have in the delta functions m loop momenta contracted with v µ 1 and n − m − 1 loop momenta contracted with v µ 2 .This yields at 0SF a trivial dependence on γ of the integral.At 1SF the functions space becomes more complex due to a single loop momentum being contracted with the velocity of the first worldline.At 2SF order we would have two loop momenta contracted with v µ 1 and v µ 2 respectively, and see a significant increase in complexity.
At 5PM-1SF we face four-loop integrals with 13 propagators and 9 irreducible scalar products, cp.Fig. 1, whose reduction to master integrals poses a significant challenge.We use Kira [101,102] to perform this integrationby-parts reduction to master integrals (MIs).We encounter up to nine scalar products in the numerator and up to eight powers (7 dots) of D's in the denominators, i.e. n i/IJ ∈ [−9, 8] in Eq. ( 5).The IBP reductions utilize FireFly [107,108], a library for reconstructing rational functions from finite field samples generated with Kira.
Several new strategies have been implemented to decrease the runtime of numerical evaluations in an IBP reduction.The first key concept builds upon the modification of the Laporta algorithm [109].For every sector with n absent propagators compared to the top-levelsector we generate equations with the total number of al-lowed scalar products reduced by n.This approach yields a remarkable 10 (L−1) runtime improvement compared to the current implementation of the Laporta algorithm in Kira.The incorporation of this feature is planned for a future release of Kira 3.0 [110].
We further observe that the IBP vectors used to formulate equations exhibit a useful feature.To reduce a large number of scalar products on linear propagators, it is sufficient for the IBP system to close by seeding at most two scalar products on propagators associated with a graviton.When reducing a high number of dots on linear propagators it is not necessary to seed dots on the graviton propagators.Implementation of this feature results in an additional factor of 10 in main memory management improvement.The complete IBP reductions took around 300k core hours on HPC clusters.Both the IBP reductions and the impulse were also assembled with the aid of Kira.
-After IBP reduction we find a total of 236 + 234 master integrals (MIs), which are solved using the method of differential equations (DE) [111,112].The needed top sectors of MIs are pictured in Fig. 1.Grouping them in a vector I that obeys d dx I = M (x, ǫ)I, we seek a transformation J = T (x, ǫ)I such that the DE factorizes: where x = γ − γ 2 − 1, which is chosen to rationalize the equation.To simplify this task it is important to first find a basis in which the dependency on γ and ǫ factorizes [113,114].For this it is helpful to allow for derivatives on the delta functions.We employ the algorithms CANONICA [115], INITIAL [116] and FiniteFlow [117].More details on these transformations were given in Ref. [118].The function space of the integrals of the (I) family (5), which are needed for the conservative calculation, is the same as at 4PM order [38], being comprised of iterated integrals of the integration kernels { 1 x , 1 1±x , x 1+x 2 } plus kernels containing K(1 − x 2 ) 2 , K being the complete elliptic integral of the first kind.
Boundary integrals.-From the solution of the ǫfactorized DE the master integrals are determined up to integration constants.We fix these in the static limit (γ → 1, v → 0) using the method of regions [119][120][121] by expanding the integrand in v.The regions are characterized by different relative scalings of the bulk graviton loop-momenta of their spacial and timelike components: referred to as potential (P) or radiative (R) modes.
Only gravitons which may go on-shell can turn radiative, i.e. the three propagators {D 12 , D 13 , D 14 } of Eq. (5c), typeset in red in Fig. 1.We hence encounter four possible regions (PPP), (PPR), (PRR) and (RRR).Our definition of conservative dynamics involves taking the evenin-velocity part, hence we consider only the (PPP) and (PRR) regions which have this scaling.The (PRR) region comes with an overall velocity scaling of (1 − x) −4ǫ , which accounts for the tail effect and all log(1 − x) functions in the final result.The 236 + 234 MIs reduce after IBP reduction of their static limits to only 2 + 1 boundary integrals in the (PRR) and 14+12 integrals in the (PPP) region.We solve the (PPP) integrals up to cuts by applying Eq. ( 6) in reverse; partial fraction identities then constrain their values, making them expressible in terms of Γ functions.Interestingly, two-worldline integrals are not fully constrained by this approach, yet appear in linear combinations such that the unknown factor cancels out in the final result.We are not able to reduce the (PRR) integrals using cuts and their expressions are more complicated, involving hypergeometric functions.Function space.-Surprisingly, the resulting function space is simpler than anticipated.The answer for ∆p (5)µ 1SF in the b µ direction is given by multiple polylogarithms (MPLs) [122][123][124] up to weight 3.These MPLs are defined by with G( 0 n , y) = log n (y)/n! and a i , y ∈ C.Even though we encounter the known elliptic integration kernels in the DEs of these integrals, they only contribute to the answer at O(ǫ), and thus disappear once we take the limit D → 4. In the final result complete elliptic integrals of the first and second kind appear only in the v-direction in the combinations K(1 ).In fact, the v-direction component is entirely determined by lower-order PM results upon momentum conservation.As we shall see below, the function space of the scattering angle therefore consists only of MPLs.
Results.-We begin with the 5PM-1SF momentum impulse ∆p with the basis vectors α (γ) are rational functions (up to integer powers of γ 2 − 1).For the explicit expressions we refer the reader to the ancillary file.The non-trivial γ dependence is spanned by the functions F (ρ) α (γ) that take the surprisingly simple form where the 31 functions f k (γ) are given by MPLs up to weight three, explicitly stated in Table I in the supplementary material, and g k (γ) involve MPLs up to weight two known from the 4PM scattering angle [55].We choose to present our results in terms of y = 1 − x, the five-letter alphabet (shifted with respect to the DEs) then being {0, 1, 2, 1 ± i}.This avoids a proliferation of ζ-values, and renders the small-velocity expansion more natural.Complex arguments always appear in conjugate combinations, such that the imaginary part cancels.We also present details on the 0SF computation that was done as a test bed in the supplementary material.All our results are collected in an ancillary file attached to the arXiv submission of this Letter.The conservative scattering angle θ cons may be extracted from the impulse using and the total mass is M = m 1 + m 2 , with ν = m 1 m 2 /M 2 the symmetric mass ratio.The scattering angle may then be double expanded as where n denotes the PM, m the SF orders and we use the floor function ⌊.⌋.The central result of our work is the 5PM-1SF contribution that takes the form where f k (γ) are the linear combinations of MPLs up to weight three and c k (γ) are polynomials in γ except for integer powers of γ 2 − 1 = γv and γ.Notice here the total absence of elliptic functions.Both the functions f k (γ) and the coefficients c k have a definite parity under v → −v such that the angle has even parity (up to factors of log(v)).They are explicitly stated in Tables I and II of the supplementary material.
Checks.-As validation of our result for the impulse ∆p µ 1 , the following checks were successfully performed: (1) total momentum conservation p 2 1 = (p 1 + ∆p 1 ) 2 , (2) reproduction of the geodesic motion (0SF), (3) agreement in the v → 0 limit with the scattering angle up to 5PN order [125,126]: with the velocity v = γ 2 − 1/γ.Finally, (4) the discontinuity of the scattering angle in the complex plane γ ∈ C is given by the radiated energy at one order lower in the PM expansion [34,54,55,127,128]: with the total angular momentum L = p ∞ |b|.This operation picks out the coefficient of log(γ − ) = log(γ − 1), with the branch cut naturally extending along the negative real axis.Given that it is by definition even-in-v, our conservative angle matches the odd-in-v part of the 4PM radiated energy E rad (L being odd in v).Upon including dissipative effects in the scattering angle, we anticipate a match to the full radiated energy.With our new 5PM-1SF result, we have verified Eq. ( 17) to the corresponding order with the 4PM-accurate loss of energy on the right-hand-side [37,54,84].
Outlook.-In this Letter we have computed the first complete results for scattering observables involving nonspinning black holes and neutron stars at 5PM (G 5 ) order -the 1SF component.This was an exceptionally challenging calculation requiring advances in IBP technology plus high-performance computing.The biggest surprise, given the appearance of elliptic E/K functions at 4PM order, was the total absence of these terms in the 5PM-1SF scattering angle, which consists only of MPLs up to weight three.This happens despite these functions appearing in the corresponding DEs.Having so far focused on the purely conservative sector, the question now arises whether this pattern persists when dissipative effects are also included.It will also be fascinating to see whether the Calabi-Yau three-fold, which appears in the DE of the dissipative effects [118], contributes to the full answer.
Looking further ahead, our main challenge will be to complete 5PM with the missing 2SF component.This represents another leap in complexity.Nevertheless, it is an important task: with a complete knowledge of the 5PM scattering dynamics (including spin, which appears at lower loop orders) our results will full encapsulate the 4PN conservative two-body dynamics.Our scattering angle is in one-to-one correspondence with a hyperbolic two-body Hamiltonian: given recent promising work on mapping unbound to bound orbits in the presence of tails [129], there is a prospect of incorporating our results into future-generation gravitational waveform models.Resummation into the strong-field regime for scattering events using EOB [99,100] will also likely show further improvements with respect to NR.

SUPPLEMENTARY MATERIAL
Probe limit.-In order to reproduce the probe limit using our formalism, we also computed the 63 Feynman diagrams contributing to the 5PM-0SF part of the momentum impulse.Following the same steps on integrand generation and tensor reduction we arrive at the 5PM-0SF probe integral family with the propagators (j = 1, 2, 3, 4 and ℓ ij...k := ℓ i + ℓ j + . . .+ ℓ k ) A difference with respect to the 1SF calculation was that we did not apply partial fraction identities: while we could, in principle, have reduced to a planar integral family similar to the one used at 1SF (5), we found it simpler to use the integral family P R instead.Note that these integrals are independent of γ and evaluate to pure functions of the dimensional regulator ǫ.They were IBP reduced to a set of 8 + 7 MIs using Kira.These MIs also feature as boundary integrals in the 1SF computation and needed to be evaluated.The final result for the probe limit momentum impulse is very compact, and takes the form ∆p The resulting scattering angle agrees with the result of Ref. [105] at 5PM-0SF order.
Basis functions and coefficients.-In Table I we list the 31 basis functions that span the 5PM-1SF scattering angle (11).They are polynomials in MPLs based on the alphabet {0, 1, 2, 1 ± i} and pure weight functions of weight 0, 1, 2 and 3 with multiplicities 1, 2, 6 and 22 respectively.In addition they possess a definite parity under v → −v for which γ 2 − 1 = γv changes sign.In Table II we list the corresponding coefficient polynomials of Eq. ( 14) for the 5PM-1SF scattering angle.For an efficient numerical evaluation of MPLs see PolyLogTools [124] using GiNaC [130].