Linear Response, Hamiltonian and Radiative Spinning Two-Body Dynamics

Using the spinning, supersymmetric Worldline Quantum Field Theory formalism we compute the momentum impulse and spin kick from a scattering of two spinning black holes or neutron stars up to quadratic order in spin at third post-Minkowskian (PM) order, including radiation-reaction effects and with arbitrarily mis-aligned spin directions. Parts of these observables, both conservative and radiative, are also inferred from lower-PM scattering data by extending Bini and Damour's linear response formula to include mis-aligned spins. By solving Hamilton's equations of motion we also use a conservative scattering angle to infer a complete 3PM two-body Hamiltonian including finite-size corrections and mis-aligned spin-spin interactions. Finally, we describe mappings to the bound two-body dynamics for aligned spin vectors: including a numerical plot of the binding energy for circular orbits compared with numerical relativity, analytic confirmation of the NNLO PN binding energy and the energy loss over successive orbits.

The need for accurate waveform templates for comparison with gravitational wave signals coming from the LIGO, Virgo and KAGRA detectors of binary merger events [1][2][3][4][5][6] -and in the future LISA, the Einstein Telescope and Cosmic Explorer [7] -has provoked enormous interest in the gravitational two-body problem. One of the most important physical properties influencing the paths of massive objects following inspiral trajectories, which as they accelerate produce gravitational waves, is their spins. Accurately determining the spins of black holes and neutron stars in binary orbits yields crucial information about their origins: if the spins are approximately aligned with the orbital plane, then this suggests formation of the binary system by slow accretion of matter; if they are mis-aligned (precessing), then this indicates formation of the binary by a random capture event.
However, an excellent alternative approach to the bound two-body problem comes by way of studying twobody scattering: here it is natural to define gaugeinvariant scattering observables in terms of the states at past-/future-infinity, where the gravitational field is weak. It is also natural here to adopt the post-Minkowskian (PM) expansion in Newton's constant G, * gustav.uhre.jakobsen@physik.hu-berlin.de † gustav.mogull@aei.mpg.de which resums terms from infinitely high velocities in the post-Newtonian (PN) series. One may use analytic continuation to directly produce PM observables for bound orbits [34][35][36][37]; alternatively, conservative scattering observables may be used to infer a Hamiltonian for the two-body system [38][39][40][41][42][43][44]. A more sophisticated version of this strategy is to infer an effective-one-body (EOB) Hamiltonian [45][46][47][48][49], which may be extended to include spin [50][51][52][53][54] and resums information from the test-body limit.
In this paper, we upgrade these observables to include radiation-reaction (dissipative) effects, using the Schwinger-Keldysh in-in formalism [76][77][78][79][80] that has recently been incorporated into both the WQFT and PMbased worldline EFT frameworks [63,67]. Our results confirm the radiated four-momentum P µ rad recently predicted with the worldline EFT approach [75]. Given these new observables, we postulate and confirm an extension to Bini and Damour's linear response relation [81][82][83] which allows us to predict terms in the conservative and radiative parts of the full scattering observables, depending on their behavior under the time-reversal operation v µ i → −v µ i . The extension holds for arbitrary spin orientations, and goes beyond linear response.
Most notably, a 3PM quadratic-in-spin Hamiltonian has now been derived using amplitudes-based methods [139], involving spin on one of the two massive bodies only and without finite-size corrections. In this paper, using the conservative scattering observables derived in Ref. [61] we both confirm this result and extend it to include spin-spin effects and finite-size corrections relevant for neutron stars. Quite remarkably, we find that knowledge of a single scattering angle suffices to completely determine the Hamiltonian, also when arbitrarily mis-aligned spin vectors are involved.
Our paper is structured as follows. In Section I we review the dynamics of spinning massive bodies, including their description up to quadratic order in spin in terms of an N = 2 supersymmetric worldline action. We demonstrate how, with a suitable SUSY shift, we can switch between the canonical and covariant spin-supplementary conditions (SSCs). In Section II we review the Schwinger-Keldysh in-in formalism in the context of WQFT, and in Section III put it to use deriving the complete 3PM quadratic-in-spin momentum impulse ∆p µ 1 and spin kick ∆S µ 1 including radiation-reaction effects. We present the results schematically, demonstrate how one may introduce scattering angles for mis-aligned spins, and perform various consistency checks.
Next, in Section IV we upgrade the linear response relation to mis-aligned spin directions, generating both conservative and radiative terms from the full 3PM scattering observables ∆p µ 1 and ∆S µ 1 . In Section V we use the conservative scattering observables, and in particular the scattering angle, to build a complete 3PM quadratic-in-spin Hamiltonian. Finally, in Section VI we discuss unbound-to-bound mappings for the specific case of aligned spins: we generate the binding energy for circular orbits, both numerically and analytically and up to 4PN order, and produce plots of the binding energy as a function of the orbital frequency close to merger -comparing our results with numerical relativity. We also determine the energy radiated per orbit using an ap-propriate analytic continuation [37]. In Section VII we conclude.

I. SPINNING MASSIVE BODIES
A pair of black holes or neutron stars interacting through a gravitational field in D-dimensional Einstein gravity are described by where S EH is the Einstein-Hilbert action (κ = √ 32πG), S gf is a gauge-fixing term and S (i) are the two worldline actions. Up to quadratic order in spin [57,58] where the projector is P i,ab := η ab − e aµ e bνẋ µ iẋ ν i /ẋ 2 i , η ab is the (mostly-minus) Minkowski metric and E i,ab := R aµbνẋ µ iẋ ν i /ẋ 2 i . The finite-size multipole moment coefficients C E,i are defined such that C E,i = 0 for black holes, and The two bodies with masses m i have positions x µ i (τ i ); the complex anticommuting fields ψ a i (τ i ), defined in a local frame e a µ (x) with g µν = e a µ e b ν η ab , encode spin degrees of freedom.
The worldline action (3) enjoys a global N = 2 supersymmetry: with constant SUSY parameters i and¯ i = † i . As shown in Ref. [58], these shifts are generated by the conserved superchargesẋ i · ψ i andẋ i ·ψ i . There is also a U(1) symmetry: generated by the conserved charge ψ i ·ψ i . Lastly, reparametrization invariance of the worldlines in τ i implieṡ is also preserved. Asẋ 2 i = 1 generically along the worldlines this implies that τ i are not the proper times; however, as we are generally only interested in the asymptotic behavior this subtlety will not be important.

A. Background Symmetries
Fields are perturbatively expanded around their background values at past infinity: where p µ i = m i v µ i is the initial momentum; the initial value of the spin tensor is given by The antisymmetrization [ab] includes a factor 1/2 -note that this normalization of the spin tensor differs from that used in Refs. [57,58,61]. The vierbein is similarly expanded as which allows us to drop the distinction between spacetime µ, ν, . . . and local frame a, b, . . . indices. The global N = 2 SUSY in the far past is To fix these symmetries we find it convenient to enforce the covariant spin-supplementary condition (SSC): Using the reparametrization symmetry we also enforce v 2 is the impact parameter pointing from the first to the second massive body. Finally, γ = v 1 · v 2 ; we will also make use of unitnormalized "hatted" variables, e.g.b µ = b µ /|b|.
The total initial angular momentum of the system is where L µν is the orbital component. In this context, we see that the background symmetries (11) correspond simply to invariance of the system's total angular momentum under shifts in the origins of the two bodies b µ i . We can also shift the center of our coordinate system x µ → x µ + a µ , in which case L µν → L µν + 2a [µ P ν] as discussed in e.g. Ref. [140]. The orbital and spin angular momentum vectors, defined specifically in D = 4 dimensions, are invariant under these shifts: where are respectively the total and center-of-mass (CoM) momentum, p µ = (0, p ∞ ). Here E = |P | = M Γ = M 1 + 2ν(γ − 1) is the energy in the CoM frame, M = m 1 + m 2 , ν = µ/M = m 1 m 2 /M 2 are the total mass and symmetric mass ratio; p ∞ = |p ∞ | = µ γ 2 − 1/Γ is the center-of-mass momentum. With the covariant SSC choice (12) the total angular momentum J µ is given by Notice that J µ = L µ + S µ 1 + S µ 2 , which is due to S µ i being defined in their respective inertial frames v µ i rather than the center-of-mass frameP µ .

B. Canonical Spin Variables
We also find it useful to introduce canonical variables [50,51,141] which are designed to ensure that and P · L can = P · S i,can = 0. The canonical spin vectors S µ i,can are given by a boost of the covariant spin vectors S µ i to the center-of-mass frame: where γ i =P · v i is the time component of v µ i in the center-of-mass frame. To ensure preservation of the total angular momentum J µ (16), we have The canonical impact parameter b µ can -in terms of which We then have, with and can confirm that the the canonical spin tensor i,can satisfies the canonical Pryce-Newton-Wigner SSC [142][143][144]: The canonical spin vector is then also given by and has a vanishing time component in the center-ofmass frame: P · S i,can = 0. This will be useful when we construct a Hamiltonian in Section V.

II. WQFT IN-IN FORMALISM
Complete observables including both conservative and radiative contributions are produced from WQFT using the Schwinger-Keldysh in-in formalism [76][77][78][79][80]. Formally this involves doubling the degrees of freedom in our theory: Observables are defined in terms of a path integral including two copies of the action: where A = 1, 2 and we use the shorthand The boundary conditions on h Aµν , z µ Ai and ψ µ Ai are that all fields equate at future infinity, and vanish at past infinity: This entangling of the boundary conditions gives rise to off-diagonal terms in the propagator matrices involving the doubled fields. For full details, see Ref. [63]. Fortunately, when performing calculations there is no need to double degrees of freedom in this way. The key insight of Ref. [63] was that tree-level single-operator expectation values (24) are produced using precisely the same Feynman rules as in the in-out formalism, but with retarded propagators pointing towards the outgoing line. The retarded graviton propagator is where P µν;ρσ := η µ(ρ η σ)ν − 1 D−2 η µν η ρσ and i0 denotes a small positive imaginary part. For the worldline modes z µ i and ψ µ i the retarded propagators are respectively The Feynman vertices are unchanged with respect to the in-in formalism: for example, the single-graviton emission vertex from worldline i is where δ − (ω) := 2πδ(ω). At tree level the WQFT simply provides a diagrammatic mechanism for solving the classical equations of motion in momentum space, and so the use of retarded propagators ensures that boundary conditions are fixed in the far past.

III. RADIATIVE OBSERVABLES
Building on Ref. [61], we compute the momentum impulse and change in ψ µ i : but now also including radiation-reaction effects. Using the definitions of the spin tensor S µν i (9) and spin vector S µ i (14b) we can then also derive the spin kick ∆S µ i : We seek the 3PM components in a PM expansion: where ∆X could be any of these observables: ∆p µ i , ∆S µ i , ∆S µν i or ∆ψ µ i . The relevant Feynman diagrams for both calculations are drawn in Fig. 1. These diagrams make no distinction between conservative and radiative effects. As only the m 2 1 m 2 2 component of ∆p  are specifically affected by the inclusion of radiation-reaction effects we recompute these components; for the rest, we simply bring forwards our previous results from Ref. [61]. Integrands are assembled using the WQFT Feynman rules [58] in D = 4 − 2 dimensions, which involves integration on the momenta or energies of all internal lines; vertices contain either energy-or . Diagrams (a)-(v) were already present in the conservative calculation [61], though their expressions are modified by the inclusion of radiation; the mushrooms (w)-(ff) are purely radiative, and did not appear in the strictly conservative calculation [61]. Diagram (y) includes the same worldline propagator with opposite i0 prescriptions, and so belongs to the K integral family (35). For brevity we use solid lines to represent both propagating deflection z µ i and spin modes ψ µ i .
momentum-conserving δ − -functions, whichever is appropriate. The energy integrals, corresponding to internal propagation of z µ i or ψ i modes, are trivial: conservation of energy at the worldline vertices resolves them immediately.
Each graph has three unresolved four-momenta to integrate over. The first of these integrals is a Fourier transform: where q µ is the total momentum exchanged from the second to the first worldline, and q := d D q/(2π) D . Here we have implicitly defined the momentum-space observables ∆X(q µ , v µ i , S µν i ), which are given as linear combinations of two-loop Feynman integrals: These integrals with retarded propagators were discussed at length in Ref. [63]: propagators D 4 -D 7 are prevented from going on-shell by the requirement that 1 · v 2 = 2 · v 1 = 0, so we can safely ignore their i0 prescriptions.
We also require which accounts for the possibility of a worldline propagator appearing twice, but with different i0 prescriptions -diagram (y) in Fig. 1. The subsequent integration steps were discussed in Ref. [61], and are not substantially different with the inclusion of radiation-reaction effects in the observables. Tensorial two-loop integrals are reduced to scalar-type by expanding on a suitable basis, and then reduced to master integrals using integration-by-parts (IBP) identities. Expressions for these master integrals were provided in Ref. [63], and once the Fourier transform (33) has been performed on the exchanged momentum q µ we are left with the observables in D dimensions. The scalar integrals themselves have simple reality properties: i.e. they are either purely real or imaginary, depending on whether they have an even or odd number of worldline propagators respectively. While this implies that the momentum-space observables ∆X(q µ , v µ i , S µν i ) are complex functions, the Fourier transform (33) introduces additional factors of i, giving rise to purely real observables

The final observables ∆p
(3)µ i and ∆ψ (3)µ i in four dimensions are found by taking the limit D → 4, checking to ensure the cancellation of all poles in the dimensional regularization parameter = 2 − D 2 . We generated our integrand in D = 4 − 2 dimensions in order to account for possible cancellations of poles with the two-loop integrals. Like in Ref. [61] we have verified that the supercharges p 2 i , ψ i ·ψ i , p i · ψ i and p i ·ψ i are conserved. This means that the following identities: are satisfied.

A. Results
We find it convenient to decompose the observables into four gauge-invariant parts each: The split into conservative 'cons' and radiative 'rad' pieces is done with respect to the integrals (34), (35): the potential and radiative regions [63]. Meanwhile the split into (±) sectors is defined with respect to behavior under a time-reversal operation: which flips the signs on the timelike vectors v µ i (and the momenta p µ i = m i v µ i ) but not the spacelike vectors b µ and a µ i , leaving γ = v 1 · v 2 invariant. Under this operation, S µν i = µν ρσ p ρ i a σ i changes sign. For the impulse, we have ∆p where is a universal prefactor, and v = γ 2 − 1/γ. The vectors are given by with basis elements even/odd under time reversal: where i, j, k = 1, 2 and1 = 2,2 = 1,L µ = L µ /|L| and b µ = b µ /|b| being unit-normalized vectors. Except for denominators of (γ 2 − 1) n/2 the remaining scalar components f (±) n are polynomials in γ, C E,i , so we refrain from providing explicit expressions in the text; instead, we refer the reader to the ancillary file attached to the arXiv submission of this paper for full expressions.
Let us remark on certain properties of this result. ∆p thus arise from the overall factors of iπ in the purely imaginary master integrals -see Ref. [63]. We also note the behavior under v → −v: in each case, the radiative components pick up the opposite sign from the conservative components. Finally, the function I(v) is familiar: it appears in the 2PM radiated angular momentum (57). As we shall see in Section IV, ∆p i,cons can be inferred directly from lower-PM observables, using a generalization to the linear response relation [81][82][83].
From the impulse ∆p µ i we straightforwardly recover the four-momentum radiated from the scattering event: which vanishes if we consider only conservative scattering. Here we agree with a recent 3PM worldline EFT result for P µ rad obtained by Riva, Vernizzi and Wong [75]. We also agree with our own previous result for the leading-order radiated energy in the CoM frame E rad =P · P rad , produced in collaboration with Plefka and Steinhoff [57], in which performing the required integrals necessitated a PN expansion -now we no longer need to do so. While knowledge of both θ and P µ rad allows one to reconstruct ∆p µ i for aligned spins, this is not true in general -thus, in this work we fill in the missing pieces from Ref. [75]. However, we note that a corresponding expression for J µ rad at 3PM order is still lacking (the leading-order 2PM is known, see Eq. (57) below).
The 3PM spin kick takes a similar form as the impulse: but in this case ∆S is associated with real integrals and ∆S (3;+)µ 1 with imaginary integrals (36). The vectors d (±)µ n are given by with basis elements ρ µ ± even/odd under time reversal: 1 Only the basis elements of this list with at least one factor of a µ 1 are relevant for the kick ∆S µ 1 . Again, we refer the interested reader to the ancillary file for fully explicit results.

B. Scattering Angles
Conservative dynamics with spin vectors aligned to the scattering plane are described by a single angle: where p µ is the CoM momentum (15b). However, generic spins and radiative effects both require a generalization of this simple parameterization. Generic spins result in non-planar motion and a non-zero spin kick; radiative effects in loss of total four-momentum P µ rad -we will present a more generic parametrization below (58).
For generic mis-aligned spin directions the impulse and spin kick are parameterized in spherical coordinates in terms of several angles. We focus on the following two, including radiation-reaction effects: The conservative counterparts of these angles, θ cons and φ cons , are defined by instead inserting ∆p µ i,cons on the right-hand side. In this case the particle label on the angles is superficial as ∆p µ 1,cons = −∆p µ 2,cons . For aligned spins at 3PM order these two definitions are equivalent to each other: θ i = φ i , which using Eq. (48) holds to all PM orders for strictly conservative scattering. Up to linear order in spin the angles equate even for misaligned spin vectors: , as dependence on the spin vectors S µ i only enters through the spin-orbit terms L · S i . Surprisingly though, and in contrast to the use of spherical coordinates, we will find that only one of these angles suffices to fully describe the conservative impulse and spin kick.
In the conservative case the interpretation of θ cons and φ cons is simple. First, φ cons measures the total scattering angle in the CoM frame in the plane spanned byb µ and p µ . Second, θ cons measures the total angle between the initial and final momentum (which may point out of the initial plane). Including radiation-reaction we may still compute θ i and φ i using Eqs. (49) although their physical interpretation as scattering angles is less clear. We also note that, while θ i is independent of the choice of SSC, φ i is not due to the manifest dependence on b µ , which transforms under SUSY shifts. To put it another way: the notion of initial plane of scattering depends on the SSC. For this reason we will mostly focus on θ cons and later use it to parameterize the Hamiltonian in Section V.
We expand the angle θ i in G and spins: Here i and j take the values 1,2. The coefficients θ (n;A) are functions only of γ, ν and C E,i . Note that only the final coefficients θ (n;2,4,i,j) k depend on the particle label. The coefficients θ (3;A) are provided in Appendix B and full expressions for the angles in the ancillary file. Expansion coefficients for θ cons , i.e. θ (n;A) cons , are defined in an equivalent manner.
For aligned spins we verify our results for θ 1 = θ 2 against several results from the literature. First, we re-produce the result of our earlier work [61] where the radiative part of θ i was computed using linear response (see Section IV) and has subsequently been extended to all spin orders by Alessio and Di Vecchia [145]. Second, we match our results for the probe limit and comparable-mass PN results [52,146,147]. Finally, for mis-aligned spin vectors in the high-energy limit where we let γ → ∞ while keeping E and a µ i constant we recover a finite result: Here we use the notation a µ ± = a µ 1 ± a µ 2 . Cancellations between the conservative and radiative pieces are essential to ensure the finiteness of this result. Note the dependence on the particle label in the second term of the second line, which disappears for Kerr black holes.
Finally, let us discuss parameterizations of the full radiative observables. One may introduce Lorentz transformations Λ i µ ν that transform the initial momenta and spin vectors to the final ones [51,53,148]: The same transformation acts on both p µ i and S µ i : one sees this naturally given the requirement that p 2 i , S 2 i and p i · S i must all be explicitly conserved. Conservation of p 2 i (S 2 i ) implies that the final momentum (spin vector) is given by a boost (rotation) of the initial one. Conservation of p i · S i implies that the boost and rotation may be combined into a single Lorentz transformation.

IV. BEYOND LINEAR RESPONSE
The Bini-Damour linear response relation is used to infer linearly radiative contributions to the scattering angle θ from (angular) momentum loss at lower-PM orders. For aligned spins [81][82][83] where J rad and E rad are respectively the angular momentum and energy losses. In Ref. [61] this was used to deduce the radiative part of the quadratic-in-spin 3PM scattering angle θ rad for aligned spins, which we have now re-confirmed with our full calculation of ∆p . At 3PM order, given that E rad ∼ O(G 3 ) the second term plays no role: the linear response is entirely accounted for by the radiated angular momentum J rad ∼ O(G 2 ).
Using our newly derived observables we have checked and can confirm that at 3PM order the linear response relation (53) generalizes for mis-aligned spin vectors to Again, the part of this formula carrying P µ rad vanishes at 3PM order, but we include it to maintain the analogy with Eq. (53); in Eq. (38a) we defined ∆p (+)µ i as the part of the impulse even under a time-reversal operation, where v µ i → −v µ i . As a special case of Eq. (54) we derive a linear response relation for θ i (49a) at 3PM order: where θ i,rad = θ i − θ cons and the (+) superscript was defined in Eq. (39).
For the spin kick we learn about the odd part ∆S The J ν -derivative is equivalent to an L ν -derivative (16): when taking these vectorial derivatives, we ignore all constraints (e.g. L · p i = 0) and treat the vectors involved (L µ , p µ i , a µ i ) as independent. All dependence on b µ should be re-expressed in terms of L µ by inverting Both of these linear response relations involve the full vectorial radiated angular momentum J µ rad . We require it only up to 2PM order: where v = γ 2 − 1/γ is the relative velocity, a µ 3 = a µ 1 + a µ 2 and the complex vector is ζ µ =L µ +ib µ ; the universal prefactor I(v) was given in Eq. (41).
For zero or aligned spins, one may straightforwardly show that the new linear response relation (54) reduces to the Bini-Damour formula (53) by inserting [36]: which holds up to the desired 3PM order. This schematic form of the impulse shows that ∆p µ 1 is fully characterized by θ and P µ rad : as ∆p (+)µ = p ∞ sin θb µ knows only about the scattering angle, the linear response relationship yields no information concerning P µ rad . We use the fact that, for aligned spins, J µ rad = −J radL µ . However, in this section we generalize beyond simple linear response. The linear response relation (54) forms part of a more general pair of relationships that allow us to reconstruct conservative and radiative parts of the scattering observables at higher-PM orders: for the impulse, and for the spin kick; we will define the split ∆X = ∆X cons + ∆X rad of the full observables into conservative and radiative parts below. We interpret all observables as real functions of the initial kinematic vectors: the orbital angular momentum vector L µ (14a), the spin vectors S µ i = m i a µ i , and the momenta p µ i = m i v µ i .

A. Derivation
We define the conservative part of a single-operator expectation value as the average of its value evaluated in the in-in and out-out prescriptions: The expectation O out−out is computed using precisely the same Feynman rules as O in−in , but with advanced propagators pointing towards the outgoing line instead of retarded, both on the worldlines and in the bulk. At 3PM order, one may verify by explicit calculation that this definition of the conservative dynamics coincides precisely with evaluating integrals only in the potential regionthe approach previously taken for the conservative 3PM spinning dynamics in Ref. [61].
Specializing to the impulse ∆p µ i , we therefore have ∆p µ i,cons While ∆p µ i,in−in is given in terms of background parameters defined at past infinity (− subscript), ∆p µ i,out−out is evaluated in terms of parameters defined at future infinity (+ subscript). As we prefer to express ∆p µ i,cons in terms of initial variables we insert Thus p µ i+ and S µ i+ are given using our pre-existing knowledge of the impulse ∆p µ i and spin kick ∆S µ i . We can infer ∆L µ -and therefore L µ + -up to 2PM order from the known 2PM angular momentum loss J µ rad (57). Using Eq. (16) , which we can rearrange to find ∆L µ -ignoring the (linear) momentum loss P µ rad ∼ O(G 3 ). In the non-spinning case ∆L µ = −J µ rad , i.e. the change in the orbital angular momentum vector is given precisely by the total loss of angular momentum.
Finally, to obtain Eq. (59) we use the fact that which simply tells us that, having computed ∆p µ i,in−in , we may easily derive ∆p µ i,out−out by continuing p µ i → −p µ i . This works because the time-reversal operation induces a change of sign on the i0 prescription of the propagators (27) and (28). For the graviton propagator (27) sgn(k 0 )i0 = (k · v i )i0: the sign on the energy component of k µ is defined by the direction of either velocity vector v µ i . For the worldline propagators (28), ω → −ω: the z i propagator (28a) remains the same, but with i0 → −i0, while the ψ i propagator (28b) also picks up an overall sign. This overall sign is compensated for in the WQFT Feynman rules for by the vertices, which are themselves invariant under time reversal except for S µν i = µν ρσ v ρ i S σ i , which flips as S µν i → −S µν i . Each time we propagate an internal spin mode we pick up a factor of S µν i , which compensates for the additional sign. The derivation of Eq. (60) for the spin kick proceeds similarly, although in this case as ∆S µ i is defined indirectly via ∆p µ i and ∆S µν i (31) we obtain the different relative signs. One can see why this is necessary by examining the spin kick at 1PM order: which is of course purely conservative. It is odd under the time-reversal v µ i → −v µ i , which agrees with Eq. (60b): we have ∆S (1)

B. Interpretation
We can now derive ∆p (m;−)µ i,cons at any PM order m. To do so, we insert the PM decomposition (32) into Eq. (59a) and perform a Taylor-series expansion of the right-hand side, picking out the desired PM order m: Taking the difference between this formula and its counterpart with p µ i → −p µ i gives ∆p i,rad . However, at 3PM order an additional simplification is possible: using the fact that the conservative and radiative observables have opposite behaviors under v → −v. As J µ rad is the only non-zero radiative observable at 2PM order, it follows that all other contributions to the linear response relation cancel out at 3PM order, leaving us with Eqs. (54) and (56) as proposed earlier. This we have checked carefully by direct calculation.
We anticipate that Eq. (60) will be useful for future 4PM computations, as we will not need to calculate the complete radiative observables: for the impulse we need only calculate ∆p (4;+)µ 1,cons and ∆p (4;−)µ 1,rad directly. This cuts down on the regions within which we need to evaluate the master integrals: only the conservative sector for the real integrals, and the radiative sector for the pseudoreal integrals. However, in the radiative sector this is predicated on our knowing the 3PM angular momentum loss J (3)µ rad , which currently we have only up to the leading 2PM order (57) in the spinning case. For non-spinning bodies, the 3PM angular momentum loss has been determined and used to infer contributions to 4PM scattering observables [140]; a similar concept will certainly apply in the presence of spin.

V. HAMILTONIAN
Let us now focus on the strictly conservative part of the dynamics, encoded by ∆p µ i,cons and ∆S µ i,cons . Computing a 3PM quadratic-in-spin Hamiltonian maps these unbound observables into bound dynamics, which in the spinning context is especially useful given the current lack of direct analytic continuations between bound and unbound observables with generic mis-aligned spin directions. In line with recent literature on Post-Minkowskian dynamics [127,130,139] we work in the CoM-frame with canonical variables p(t) and x(t) describing the relative momentum and position of the two bodies respectively. These dynamical variables satisfy canonical Poisson brackets: The spin of each body is described by the spin vectors S i (t), with the canonical SSC described in Sec. I B.
The Hamiltonian H is fixed in isotropic gauge, meaning that it does not depend on x(t) · p(t). It takes the general form The potential is expanded in spin structures: where a i = S i /m i . In each case the first index counts the spin order; subsequent indices count the specific structures involved. Note that the symmetric spin structures O (2,a,1,2) = O (2,a,2,1) are counted twice and their coefficients are equal. Finally, we PM-expand each component: These coefficients c (n;A) fully encode the Hamiltonian.
We fix the coefficients c (n;A) (p 2 ) by matching observables computed from the Hamiltonian H with scattering observables from the WQFT. Hamilton's equations for the dynamical variables arė and we solve them perturbatively up to third order in G: The zeroth-order scattering trajectories are where p µ = (0, p ∞ ), S µ i,can = (0, S i,∞ ) and b µ can = (0, b can ). The dimensionless parameter ξ is defined as Inserting the expansions of the dynamical variables (74) into Hamilton's equations (73) we get perturbative equations of motion at each PM order. The spatial components of the impulse and spin kick in the CoM frame are then given by where these are the conservative 3-vector components of ∆p µ 1 and ∆S µ i,can . The change in the canonical spin vector is given in terms of the covariant spin kick by having used the definition of S µ i,can (18), and assuming conservative scattering.
Quite remarkably though, we find that knowledge of θ cons (49) suffices in order to fully fix all coefficients in the Hamiltonian. The impulse and spin kick may then be expressed in terms of the coefficients c (n;A) (p 2 ) and derivatives thereof. The derivatives come about when we evaluateẋ by differentiating H with respect to p: where c (n;A;m) := ∂ m c (n;A) /(∂p 2 ) m -the observables are written as functions of c (n;A;m) . We then match those expressions to the explicit WQFT-derived results and solve for the coefficients. We print the results until linear in spins here, and the results quadratic in spins in Appendix C: Here φ A particular subtlety here is that θ (n;A) can is given in terms of the previously-defined background variables p ∞ , γ, E i = p 2 ∞ + m 2 i evaluated at past infinity, e.g.
instead of the dynamical momentum p(t): we interpolate to the full dynamical coefficients c (n;A;m) simply by replacing one with the other. We have also introduced the differential operator and ]. The angle and its coefficients are given in Appendix B, together with full expressions for the Hamiltonian coefficients given in the ancillary file attached to the arXiv submission of this paper. We have checked this Hamiltonian numerically against the recent results obtained in Ref. [139], which also included 3PM quadratic-in-spin terms. Our results complement those by adding S 1 S 2 contributions to the Hamiltonian together with finite-size effects (C E,i terms) in the S 2 i sector. We have also verified that the PN-expansion of this Hamiltonian correctly reproduces 4PN results in the isotropic gauge [29]. We did so by PN-expanding the c (n;A) (p 2 ) coefficients in powers of p 2 .
Finally, let us observe that: having expressed the coefficients of the Hamiltonian in terms only of the scattering angle φ can , this implies that the full conservative scattering observables ∆p i,cons may themselves be expressed in terms of this angle. We would obtain the precise relationship by solving Hamilton's equations again for the impulse and spin kick, but this time plugging in expressions in terms of φ cons . In contrast to the Hamiltonian, the relations thus obtained are gaugeinvariant and will be an intriguing topic of future studies.

VI. UNBOUND-TO-BOUND MAPPINGS
Let us now discuss how our results may be applied to describe bound orbits, which the now-complete 3PM quadratic-in-spin Hamiltonian (69) gives us partial access to. This will allow us to determine the binding energy, which together with the radiative fluxes may be used to inform complete gravitational waveform models. In this section we specialize to spin vectors aligned with the orbital angular momentum vector: where χ i are the directed spin lengths and S µ i,can = S µ i . For Kerr black holes m i |χ i | are the radii of the ring singularities, and −1 < χ i < 1. Using Eq. (19) we see that for aligned spinsL µ =L µ can ; however, their magnitudes differ and so using Eq. (19) we introduce where E = (E − M )/µ is the reduced binding energy, δ = (m 2 − m 1 )/M and Using axial symmetry the Hamiltonian E = H(r, p r , λ, χ i ) depends -besides the masses m i and finite-size coefficients C E,i -on the radial coordinate r and momentum p r , where in axial coordinates The axial momentum p φ = |L can | is a constant of motion. In all of these expressions we leave the dependence on masses m i and finite-size coefficients C E,i implicit.

A. Numerical PM binding energy
In Fig. 2 we plot the reduced binding energy E = (H − M )/µ for circular orbits as a function of the orbital frequency GM Ω leading up to merger. It is compared with a Numerical Relativity (NR) simulation provided by the Simulating eXtreme Spacetimes (SXS) collaboration [149], extracted in Ref. [150]. Our plots are also determined numerically: within the Hamiltonian we set p r = 0 and, usingṗ r = 0 = −∂H/∂r, we solve for λ(r) for different orbital separations r and with specific values of χ 1 , χ 2 . The reduced binding energy E is plotted against the orbital frequency: with the number of orbits leading up to merger provided by NR -see Refs. [48,71] for more details. The conclusion of these plots is somewhat disappointing: the quadratic-in-spin part of the Hamiltonian yields little improvement over the spin-orbit contribution. However, there is a far more noticeable improvement when going from 2PM to 3PM order, which suggests that producing a 4PM spinning Hamiltonian will be a worthwhile endeavor. A similar improvement of the 4PM (hyperbolic) Hamiltonian over the 3PM seen in the nonspinning case [71] also encourages us in this direction; however, it will also be important to resum in the testbody limit by feeding these results into a suitable EOB model [45]. It is also worth noting that, as these plots are generated for circular orbits, they do not showcase PM results in the best possible light: it is anticipated that PM-based results will perform better for highly elliptical orbits, where the velocity at closest approach between the massive bodies is large [71].

B. Analytic PN binding energy
Working in the PN expansion we may derive precise analytic formulae for the binding energy E and periastron advance ∆φ. Following closely the discussion in Ref. [147] (see also Refs. [35,148]) our starting point is the radial action for unbound orbits: where Pf denotes the partie finie of the radial action; the energy constraint E = H(r, p r , λ, χ i ) can be solved for the radial momentum p r -but we refrain from doing so explicitly. The innermost point r min is given by the root of p r (r min , E, λ, χ i ) = 0. This unbound radial action is related to the bound radial action i r by analytic continuation: which sends λ → −λ and χ i → −χ i . For unbound dynamics the reduced binding energy E > 0 (γ > 1), but now we consider values E < 0 (0 < γ < 1) as we see in Fig. 2. The reduced binding energy E for circular orbits is derived by setting the bound radial action to zero: Using Eqs. (91) and (88) we may write E as a function of x, the spins χ ± , χ E,± , and ν. Rather than the Hamiltonian, we prefer to derive the radial action -and thus the binding energy E -from the gauge-invariant conservative scattering angle θ cons . It is given by a λ-derivative of the unbound radial action w r : and by analytic continuation it is related to the periastron advance for bound orbits [34][35][36][37]: We PM-expand the angle in λ as whereθ (n) depends on χ i through the ratios χ i /λ. The analytic continuation in λ (90) is trivial for all terms in w r except the one coming from the non-spinning part ofθ (1) , wherein the dependence on λ is log(λ). With this exception the odd-in-G terms in the PM expansion disappear in the difference (90), and the analytic continuation of log(−λ) − log(λ) leaves behind a finite piece. We therefore have Here γ should be re-expressed in terms of E. The integral on λ is elementary but left unresolved because of the remaining λ-dependence insideθ (2n) . Naively this result indicates that our 3PM results have no relevance for the mapping to bound results. However, this obstacle is conveniently circumvented in the PN expansion by use of the so-called impetus formula with PM-coefficients f k : This is fed into the definition of the scattering angle by way of the radial action (89): This integral has been performed up to high orders in G in Ref. [43] and the scattering angle is then expressed in terms of f k . Knowledge of the 1PM, 2PM and 3PM scattering angles suffices in order to fully determine f k≤3 , which may in turn be used to reconstruct the leading-PN terms of the angle at higher PM orders. In order to reconstruct the quadratic-in-spin binding energy up to 4PN order, we find it necessary to reconstruct the leading-PN and sub-leading-PN parts ofθ (6) andθ (4) respectively, at both linear and quadratic order in spin. Following this procedure, we discover that the binding energy is At each order in spin, we give the terms up to and including next-to-next-to-leading order in the PN expansion. With the non-spinning terms appearing up to 2PN order, the spin-orbit and spin-spin terms appear up to 3.5PN and 4PN order respectively.

C. Radiated energy
Finally, we may also determine the energy radiated per orbit in the CoM frame [37]: where E rad =P · P rad derives from the full radiative momentum impulse ∆p µ 1 (44). The same analysis was also done successfully by Riva, Vernizzi and Wong [75], who confirmed that this result for the radiated energy for bound orbits agrees with the corresponding 4PN terms from Ref. [32]. To leading-PN at each spin order, we find that Like in the binding energy (98), we see that the linearin-spin terms appear at 1.5PN above the non-spinning; the quadratic-in-spin terms at 2PN above. However, we emphasize that, unlike in the case of the binding energy, this result for the radiated energy does not encode the full leading-PN result. To see why, we recall that in the non-spinning case the leading-PN radiated energy may be derived from Einstein's quadrupole formula: While we do reproduce the λ −3 contribution, gaining access to the λ −5 and λ −7 contributions via PM-based scattering calculations would necessitate 5PM and 7PM calculations respectively: a seemingly impossible task! However, a more encouraging conclusion was reached in Ref. [36]: that by instead using the PM-scattering data to fix the form of the (gauge-dependent) instantaneous fluxes of energy and angular momentum, the leading-PN information could instead be extracted from a future 4PM derivation. Such an approach might be especially beneficial in the spinning case, given the unbound-tobound mapping's current limitation to aligned spin vectors. This would be similar to our present use of the conservative scattering observables to reconstruct a Hamiltonian, a prospect that we leave for future work.

VII. CONCLUSIONS
In this paper we have for the first time provided complete expressions for the impulse ∆p µ 1 and spin kick ∆S µ 1 at third Post-Minkowskian (3PM) order to quadratic order in the spins of two scattering bodies, including finitesize corrections. The computation relied on our use of the spinning, supersymmetric WQFT formalism [57,58] and its extension to utilize the Schwinger-Keldysh in-in formalism [63,[76][77][78][79][80]. These results upgrade the previously obtained conservative observables provided by the present authors [61], and include knowledge of the total radiated four-momentum P µ rad which has also recently been computed using worldline EFT methods [75]. We also wrote down a scattering angle that encapsulates the motion for arbitrarily mis-aligned spin directions.
Next, we demonstrated how both conservative and radiative parts of these observables may be reconstructed using an extension of Bini and Damour's linear response relation [81][82][83], incorporating the full 2PM radiated angular momentum J µ rad . These relations build on a split into conservative and radiative parts as an average of the full in-in and out-out observables. It will be exciting to see how these relations may help produce results at 4PM order. For spin effects this will require a 3PM computation of J µ rad with spin, but one may already explore non-spinning applications of the formula. Similar studies have already been initiated in the non-spinning case in Ref. [140].
Using the conservative parts of our results -already known from Ref. [61] -we constructed a two-body Hamiltonian mapping our unbound results to bound motion. This Hamiltonian describes the conservative twobody dynamics up to 3PM order and to quadratic order in their spins -it complements the Hamiltonian of Ref. [139] by adding the terms O(S 1 S 2 ) and finite-size effects. We note that the coefficients of the Hamiltonian are more complicated than the scattering observables. While the observables can easily be reduced to a number of polynomials in γ, the Hamiltonian coefficients depend on the center-of-mass variables in a complicated manner. This is partly due to its gauge dependence, but also the necessity of using a canonical spin-supplementary condition (SSC). This canonical Pryce-Newton-Wigner SSC [142][143][144] we showed is related to the covariant SSC by a supersymmetry transformation.
It is an interesting study to explore whether gauge choices other than the isotropic gauge could lead to simpler coefficients in the Hamiltonian. Quite intriguingly, we found it possible to fix all coefficients of the Hamiltonian by matching to a single scattering angle defined for generic spins. This in turn leads to the exciting result that the conservative dynamics for generic spins can be described by a single scalar. The expressions for the impulse and spin kick in terms of the scattering angle thus obtained are gauge invariant. They deserve further study and preferably a direct relation highlighting the gauge invariance. Such relations are similar to the eikonal relations explored in Ref. [127].
We also studied mappings to bound orbits. This began with using the complete 3PM Hamiltonian to produce numerical plots of the binding energy for circular orbits leading up to merger, in comparison with Numerical Relativity (NR) simulations [149]. We successfully reproduced the known 4PN quadratic-in-spin binding energy, and the leading-PM radiated energy for bound orbits. Unfortunately, these plots did not show a significant improvement of the quadratic-in-spin Hamiltonian over its spin-orbit counterpart; however, the effect of going from 2PM to 3PM order was more significant. This calls for the future determination of the 4PM spinning Hamiltonian, which -similar to the non-spinning case, and due to the presence of tails -encounters non-localities that distinguish between bound and unbound dynamics [68-71, 103, 104]. To inform realistic waveform models, it will also be important to incorporate knowledge of the test-body limit by way of the effective-one-body (EOB) formalism [45,46].
Finally, we determined the energy radiated per orbit in the CoM frame from an appropriate unbound-to-bound mapping [37]. Current limitations in this mapping restrict us to considering only aligned spins; furthermore, this approach does not reproduce the full leading-PN result -which in the non-spinning case may be derived from Einstein's quadrupole formula. To overcome both of these limitations, and following the suggestion of Ref. [36], we believe that in the future it will be more profitable to focus on reconstructing the (gaugedependent) instantaneous momentum and angular momentum fluxes. In this case, a complete 4PM result would suffice to reconstruct the leading-PN form of the radiated energy. Alternatively, we hope that improved unbound-to-bound mappings for spinning bodies will further alleviate these issues.
A natural continuation of this work will therefore be to progress upwards in the perturbative series -both to higher PM orders and higher spin orders. While higher PM orders will present a challenge for the integration step, there has recently been a promising development in this area: the first complete analytic result for the 4PM momentum impulse including radiation-reaction effects [70], wherein loop integrals with retarded propagators were also used. With a similar basis of master integrals, it will be possible to tackle spin effects at 4PM order. While higher PM orders present a challenge regarding the integration steps, higher spin orders rather challenge the construction of the integrand. In this case, it will be necessary to upgrade the spinning N = 2 supersymmetric worldline action to include more supersymmetry -a tantalizing prospect that we leave for future work.
Acknowledgments. -We thank Alessandra Buonanno, Gregor Kälin, Jan Plefka, Muddu Saketh, Benjamin Sauer, Chia-Hsien Shen and Justin Vines for enlightening discussions. We are especially grateful to Jan Steinhoff for his assistance in comparing with the 4PN quadratic-in-spin Hamiltonian, and to Mohammed Khalil for helping us produce the numerical binding energy plot. We also thank the organizers of the "High-Precision Gravitational Waves" program at the Kavli Institute for Theoretical Physics (KITP) for their hospitality. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Projekt-nummer 417533893/GRK2575 "Rethinking Quantum Field Theory". It was also supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. Here we discuss the relationship between the canonical and covariant expansions of the conservative scattering angle θ cons . As the angle is independent of the choice of SSC, only its expansion coefficients change when we expand in canonical rather than covariant variables. The canonical expansion is defined analogously to the covariant expansion in Eq. (50): θ (n;2,2,i,j) can = Γ θ (n;2,2,i,j) The covariant coefficients θ (n;A) cons are given in Appendix B.