Conservative and radiative dynamics of spinning bodies at third post-Minkowskian order using worldline quantum field theory

Using the spinning worldline quantum field theory formalism we calculate the quadratic-in-spin momentum impulse $\Delta p_i^\mu$ and spin kick $\Delta a_i^\mu$ from a scattering of two arbitrarily oriented spinning massive bodies (black holes or neutron stars) in a weak gravitational background up to third post-Minkowskian (PM) order ($G^3$). Two-loop Feynman integrals are performed in the potential region, yielding conservative results. For spins aligned to the orbital angular momentum we find a conservative scattering angle that is fully consistent with state-of-the-art post-Newtonian results. Using the 2PM radiated angular momentum previously obtained by Plefka, Steinhoff and the present authors we generalize the angle to include radiation-reaction effects, in which case it avoids divergences in the high-energy limit.

Using the spinning worldline quantum field theory formalism we calculate the quadratic-in-spin momentum impulse ∆p µ i and spin kick ∆a µ i from a scattering of two arbitrarily oriented spinning massive bodies (black holes or neutron stars) in a weak gravitational background up to third post-Minkowskian (PM) order (G 3 ). Two-loop Feynman integrals are performed in the potential region, yielding conservative results. For spins aligned to the orbital angular momentum we find a conservative scattering angle that is fully consistent with state-of-the-art post-Newtonian results. Using the 2PM radiated angular momentum previously obtained by Plefka, Steinhoff and the present authors we generalize the angle to include radiation-reaction effects, in which case it avoids divergences in the high-energy limit.
Recent detections by the LIGO and Virgo collaborations of gravitational waves emitted by binary black hole and neutron star mergers [1-5] have driven demand for high-precision gravitational waveform templates. In the early stage these inspirals typically run over many cycles, making them difficult to model using numerical techniques [6-8]; yet, as the gravitational field is weak this regime is well tackled using perturbation theory. Often this is done in a post-Newtonian (PN) expansion in both G (Newton's constant) and c (the speed of light); however, methods involving the post-Minkowskian (PM) expansion in G are gaining prominence.
The crucial insight driving this shift is that bound orbits are closely related to unbound scattering events, the latter more naturally handled in the PM expansion. A well-studied approach to the bound problem in gravity is reverse engineering a gravitational potential from scattering data [9][10][11][12][13][14][15], which can in turn be used to describe bound orbits. More recent techniques such as the Boundto-Boundary (B2B) correspondence directly relate bound with unbound observables [16][17][18]; scattering observables may also be used as direct input for an effective one-body description of the bound dynamics [19][20][21][22], also of spinning black holes or neutron stars [23][24][25][26][27][28].
The recently developed worldline QFT (WQFT) formalism [114][115][116][117] innovates over these approaches by quantizing worldline degrees of freedom. This leads to a highly streamlined PM setup wherein classical scattering observables are directly computed as sums of tree-level Feynman diagrams. The use of an N = 2 supersymmetric extension to the point-particle action to encapsulate spin degrees of freedom [116,117] circumvents the need for a local co-rotating frame. Recent work on the WQFT has included the double copy [118] and applications to light bending [119]; other closely related approaches involve directly solving the classical equations of motion [120] and Wilson line operators [121].
In this Letter we realize the spinning WQFT's full potential with a state-of-the-art calculation: deriving the quadratic-in-spin conservative momentum impulse ∆p µ i and spin kick ∆a µ i in a scattering encounter between massive bodies at 3PM order, including finite-size effects. Specializing to aligned spins yields the conservative scattering angle θ cons , which we generalize to include dissipa-tive effects using the linear response relation [122][123][124].
Spinning WQFT formalism. -The dynamics of Kerr black holes with masses m i and positions x µ i (τ ) on a curved D-dimensional background metric g µν are described up to quadratic order in spin by the N = 2 supersymmetric worldline action [125,126] The complex Grassmann-valued vectors ψ a i (τ ), defined in a local frame e a µ with g µν = e a µ e b ν η ab and Dψ a i Dτ =ψ a i + x µ ω µ ab ψ i,b , encode spin degrees of freedom (we use the mostly minus metric). The spin tensors S µν i and Pauli-Lubanski spin vectors a µ i are composite fields: i (referred to as π i,µ in Ref. [117]). Reparametrization invariance in τ and U(1) symmetry on the Grassmann vectors respectively imply conservation of p 2 i andψ i · ψ i . Global N = 2 supersymmetry provides two additional fermionic charges: p i · ψ i and p i ·ψ i , which when set to zero together imply the Tulczyjew-Dixon spin-supplementary condition (SSC) p i,µ S µν i = 0 [127,128]. The action (1) extends naturally to include finite-size objects like neutron stars by also including with projectorP ab := η ab − e aµ e bνẋ µẋν /ẋ 2 and Wilson coefficients C E,i , where C E,i = 0 for black holes. The projector ensures supersymmetry for terms up to O(S 2 ): enough to maintain the SSC and preserve lengths of the spin vectors.
The WQFT's distinguishing feature is quantization of both bulk and worldline degrees of freedom. In a weak gravitational field with κ = √ 32πG we expand g µν (x) = η µν + κh µν (x) with the vielbein e a µ = η aν η µν + κ 2 h µν − κ 2 8 h µρ h ρ ν + · · · . Thereafter we no longer distinguish between spacetime µ, ν, . . . and local frame a, b, . . . indices. The worldline fields are similarly expanded around their background values: We also define the Lorentz factor γ = v 1 · v 2 and the relative velocity v = γ 2 − 1/γ. The WQFT is defined by a path integral, with physical observables calculated as operator expectation values: We have included the D-dimensional Einstein-Hilbert action S EH and gauge-fixing term S gf to enforce ∂ ν h µν = 1 2 ∂ µ h ν ν . The stationary phase of the path integral is dominated by solutions to the the physically relevant Einstein and Mathisson-Papapetrou-Dixon (MPD) equations of motion [129][130][131]. This highlights the WQFT's main advantage when studying classical physics: the classical → 0 limit is identified with the sum of tree-level Feynman diagrams.
The WQFT Feynman rules are most naturally expressed in the Fourier domain: (2π) D and ω := dω 2π . Feynman rules for the graviton h µν originating from the bulk Einstein-Hilbert action are conventional, with propagator FIG. 1: The ten types of diagrams contributing to the m1m 3 2 components of ∆p (3)µ 1 and the m 3 2 components of ∆ψ (3)µ 1 , involving I (1;±) -type integrals (12). In the test-body limit m1 ≪ m2 these are the only surviving contributions. All graphs should be considered trees -the dotted lines represent the worldlines on which energy is conserved, instead of momentum. , involving I (2;±) -type integrals (12). We exclude "mushroom graphs" that integrate to zero in the potential region. observables: Diagrammatically this amounts to drawing all tree-level diagrams with a single cut external z µ i or ψ ′µ i line. The diagrams required to calculate both ∆p differ only by the cut outgoing line we display them together. For additional brevity we use only solid lines to represent propagating worldline modes z µ i , ψ ′µ i and ψ ′µ i ; however, it should be assumed that each internal worldline mode could be of all three types (with expressions adjusted accordingly). The third set of diagrams (not drawn) consists simply of mirrored versions of the graphs in Fig. 1 through a horizontal plane, but with the external cut line still on the first (upper) worldline. For the impulse we avoid calculating these contributions directly, instead making use of momentum conservation ∆p (for conservative scattering). We assemble expressions using the WQFT Feynman rules in D = 4 − 2ǫ spacetime dimensions, with the later intention of recovering four-dimensional results in the ǫ → 0 limit. Each retarded graviton (6) and worldline (7) propagator points toward the outgoing line: from cause to effect. As diagrams belonging to each of the three categories carry common overall factors of the masses m α 1 m β 2 the categories themselves are separately gauge invariant. This helpfully breaks the calculation up into gaugeinvariant subcomponents. Diagrams in Fig. 1 carry the maximum allowed power of m 2 , and represent the testbody limit m 1 ≪ m 2 . Integrals are performed over the energies (on the worldlines ω ) or momenta (in the bulk k ) of all internal lines.

The integrals involved in both ∆p
(3)µ 1 and ∆ψ (3)µ 1 are Fourier transforms of two-loop Feynman integrals: n1,n2,...,n7 , i = 1, 2, 3, (11) where δ − (ω) := 2πδ(ω), q µ is the total momentum exchanged from the second to the first worldline and α is an arbitrary power of |q| := √ −q · q. The two-loop integral families are n1,...,n7 v1↔v2 . Each pair (±) is associated with one of the three categories of diagrams. To achieve these representations one must first integrate on the energies carried by any internal deflection z µ i or spin ψ ′µ i ,ψ ′µ i modes on the worldlines. As two-loop integrals of this kind are now well studied -see e.g. Refs. [39,46,132,133] -we relegate full details of how to perform them to Appendix A. The I (1;±) n1,...,n7 integrals -associated with the test-body diagrams in Fig. 1 -are more straightforward, being naturally evaluated in the rest frame v µ 2 = (1, 0). The more involved I (2;±) n1,...,n7 integrals -associated with the diagrams in Fig. 2 -contain the arccoshγ function. To fix boundary conditions we adopt the potential region of integration, which ignores radiation-reaction contributions and may be interpreted as a resummation of the terms arising from a conservative PN expansion v c ≪ 1. We have therefore excluded certain graphs from Fig. 2 the so-called "mushroom graphs" -which integrate to zero within this regime.
Our final results for ∆p The coefficients c We have performed several consistency checks. Firstly, all poles in ǫ = 2 − D 2 arising from the dimensionally regularized two-loop integrals (12) are seen to cancel, thus ensuring finiteness of our results in the limit D → 4. Secondly, conservation of p 2 i ,ψ i · ψ i and the fermionic supercharge p i · ψ i between initial and final states implies a set of consistency requirements: All three of these checks are highly nontrivial: for instance, the third compares parts of ∆ψ confining the motion to a plane. The conservative part of the scattering angle is then given by (see e.g. Ref. [94]): with the full scattering angle (including radiative corrections) given by θ = θ cons + θ rad . The center-of-mass momentum is p ∞ = µγv/Γ where µ = M ν = m 1 m 2 /M is the symmetric mass, M = m 1 + m 2 is the total mass and Γ = E/M = 1 + 2ν(γ − 1), E being the total energy. We decompose the scattering angle as with n and m counting the PM and spin orders respectively. At 3PM order using our results where we have defined δ = (m 2 − m 1 )/M as well as s ± = s 1 ± s 2 and s 2 E,± = C E,1 s 2 1 ± C E,2 s 2 2 . We have checked θ cons both in the test-body limit ν → 0 and up to 4PN order (N 2 LO) for comparable masses against Refs. [26,134]. 1 For aligned spins, this provides a first check on 1 We thank Mohammed Khalil for providing us with an extension of the 4PN scattering angle to include finite-size C E,i coefficients.
the complete quadratic-in-spin conservative dynamics of compact binaries at 4PN order [107,108] together with recent work in the worldline EFT formalism [112]. As explained by Bini and Damour [122][123][124], the conservative scattering angle is generalized to include radiation using the linear response relation: Here J is the total angular momentum in the centerof-mass frame: the derivative is equivalent to one with respect to the orbital angular momentum L = p ∞ |b|.
It has recently been clarified [135] that Eq. (18) applies only using an "intrinsic" gauge choice with respect to BMS symmetry, wherein the radiated angular momen- rad we need only J (2) rad , which was provided by Plefka, Steinhoff and the present authors for arbitrary spin orientations in Ref. [116]. For aligned spins, This yields the radiative part of the scattering angle: The non-spinning part of θ rad has also been confirmed without reference to Eq. (18) -see e.g. Refs. [43,44].
A key criterion of θ rad is that the total scattering angle should remain finite in the high-energy limit. We write θ(E, ν, |b|, γ, s i ) in terms of the energy, symmetric mass ratio, impact parameter, Lorentz factor and spin magnitudes and let γ → ∞, in which case the individual masses are negligible. In this limit: While we know of no spinning extension to Amati, Ciafaloni and Veneziano's result [136] to compare with in the high-energy limit, we do see that a logarithmic divergence appearing in the conservative part of the angle (17) is canceled by the radiative correction (20).
Discussion. -We conclude with a brief discussion of bound observables. Using the B2B dictionary [16][17][18] one may, for instance, recover the aligned-spin periastron advance ∆Φ from our scattering angle: Similarly one may relate the unbound and bound radial actions, from which the scattering angle and periastron advance are respectively given by a derivative with respect to L. At 3PM order θ (3) cancels in Eq. (22); nevertheless, from θ (3) one may reconstruct the leading-PN parts of θ (4) and θ (6) (and similarly for the radial action) [134]. This suffices for a comparison with bound quadratic-in-spin results at N 2 LO: for example, we have reproduced the quadratic-in-spin N 2 LO binding energy for circular orbits [107,108] as was also very recently done in Ref. [112]. For arbitrarily aligned spins there is currently no extension of the B2B map (22). An alternative would therefore be to make an ansatz for a conservative two-body Hamiltonian -for example, building on that used at 2PM order [68,69] -and solve Hamilton's equations for comparison with ∆p µ 1 and ∆a µ 1 , thus extending those results to 3PM. On the other hand, we are hopeful that direct maps between unbound and bound gauge-invariant observables for arbitrary spins will be discovered in the near future. In that spirit, all information is captured by the impulse and spin kick.
There remains much work to be done: for example, extending ∆p to incorporate radiationreaction effects, as we have already done for the scattering angle θ  (20). This requires us to upgrade our two-loop master integrals to account for the retarded pole displacement on the graviton propagator (6) and restore the mushroom graphs to Fig. 2. We are also interested in the eikonal phase, which was computed in Ref. [117] at 2PM order as the free energy of the WQFT, and captures both the impulse and spin kick. Nevertheless, for the time being we believe that we have effectively showcased the spinning WQFT's utility and efficiency.

Appendix A: Two-loop integration
In this Appendix we outline the steps required to perform integrals of the kind appearing in Eq. (11). All scalar two-loop integrals I 1 · · · ℓ µn 1 ℓ ν1 2 · · · ℓ νm 2 ] in arbitrary D dimensions are decomposed onto bases consisting of the vectors v µ i and q µ appearing in their expressions, plus the metric. We find it convenient to introduce where by design dual vectors w µ i satisfy v i · w j = δ ij ; P µν is the metric of the (D − 3)-dimensional space orthogonal to v µ i and q µ . Then, for example The resulting integrals are straightforwardly reduced to scalar-type: n1,n2,n3,n4,n5,n6−1,n7 − |q| 2 I (i;±) n1,n2,n3,n4,n5,n6,n7 .
The highest-rank integrals we encountered in this project had five free indices: I Integration-by-parts relations. -Our next task is to find linear identities satisfied by the scalar integrals, and thus establish a minimal basis -the so-called master integrals. We do this separately for each of the six integral families I (i;±) n1,...,n7 . These linear identities are generated by integration-by-parts relations (IBPs): Following the reverse unitarity approach of Refs. [43,44] we have generalized the definition of the δ − -functions to include derivatives: where b j,k = {−1, 0, 1} on a case-by-case basis and α j are generically functions of γ, |q| and ǫ = 2 − D 2 . The two new indices n 8 , n 9 represent the propagators associated with the delta functions -in general, we always choose masters with n 8 = n 9 = 0. The IBPs separate into two categories depending on whether n 1 + n 2 is even or odd. The integral families also enjoy symmetry relations: n1,n2,n3,n4,n5,n6,n7 = I (i;±) n2,n1,n4,n3,n5,n7,n6 , I (i;±) n1,n2,n3,n4,n5,n6,n7 = I (i;±) n1,n2,n6,n7,n5,n3,n4 , where the first is due to ℓ 1 ↔ ±ℓ 2 symmetries and the second due to shifts in ℓ i by q. Using the Laporta algorithm [137,138], implemented in publicly-available packages such as FIRE [139], LiteRed [140,141] and KIRA [142,143], we reduce to minimal bases: two bases for each family, with even or odd n 1 + n 2 .
3. Insertion of master integrals. -We now provide expressions for the master integrals, beginning with the simpler I (1;±) families to all orders in ǫ = 2 − D 2 :  1,1,1,1,0 where we have set |q| = 1. As promised the dependence on γ and ǫ factorizes. These results are well-established, and may be found in e.g. Ref. [132]. The I (2;±) master integrals do not factorize, and we provide them only up to the order in ǫ to which they are required. For n 1 + n 2 even: 0,0,0,1,1,0,1 = 0 , (A9a) where the first two are known to all orders in ǫ -the integral (A9b) is a product of one-loop integrals. The dilogarithm appearing in the integrals (A9f) and (A9g) is Li 2 (2 − 2γ 2 + 2γ γ 2 − 1): this dilogarithm and arccosh 2 γ cancel from all of our final results between these two integrals. In the non-spinning part of ∆p (3)µ 1 (B1) these eight master integrals are associated with terms proportional to the impact parameter b µ . For n 1 + n 2 odd: I (2;±) 1,0,1,0,1,1,0 = 0 , 1,0,1,1,1,1,1 = 4 −1−3ǫ π −1+2ǫ e −2ǫγE i(−1 + 6ǫ) which are associated with terms proportional to the velocities v µ i . To derive these expressions for the master integrals we set up systems of differential equations (DEs). Using publicly available tools such as Fuchsia [144] and epsilon [145] one may find linear transformations to canonical bases F (x, ǫ) of master integrals that obey where M(x) is a matrix depending only on x = γ − γ 2 − 1. The use of x rather than γ in the DEs was proposed in Ref. [132]: M(x) contains only poles in {x, 1 + x, 1 − x}, i.e. the symbol alphabet. 2 The essential property of Eq. (A11) is factorization of the ǫ-dependence, which enables a straightforwards solution to these DEs for F as Laurent series expansions in ǫ = 2 − D 2 . Integration constants are fixed in the potential region by comparison with the static limit v → 0, i.e. γ → 1. To leading order in v the I (1;±) and I (2;∓) scalar integral families reduce to the same expression: We therefore fix boundary conditions on the I (2;±) integrals by equating them with the I (1;∓) integrals: expressions for the masters in this case having already been given (A8). 4. Fourier transform. -Finally we perform the Fourier transform q-integrals from Eq. (11). These generically take the form The scalar integral is well-known -see e.g. Ref. [114]: where P µν 12 := η µν − w µ 1 v ν 1 − w ν 2 v ν 2 projects to the (D − 2)-dimensional space orthogonal to v µ i . The generalization to higher-rank integrals follows easily by taking derivatives with respect to b µ : One should avoid imposing b · v i = 0 until after these derivatives have been taken -hence our use of the projector P µν 12 .