Gravitational spin-orbit and aligned spin$_1$-spin$_2$ couplings through third-subleading post-Newtonian orders

The study of scattering encounters continues to provide new insights into the general relativistic two-body problem. The local-in-time conservative dynamics of an aligned-spin binary, for both unbound and bound orbits, is fully encoded in the gauge-invariant scattering-angle function, which is most naturally expressed in a post-Minkowskian (PM) expansion, and which exhibits a remarkably simple dependence on the masses of the two bodies (in terms of appropriate geometric variables). This dependence links the PM and small-mass-ratio approximations, allowing gravitational self-force results to determine new post-Newtonian (PN) information to all orders in the mass ratio. In this paper, we exploit this interplay between relativistic scattering and self-force theory to obtain the third-subleading (4.5PN) spin-orbit dynamics for generic spins, and the third-subleading (5PN) spin$_1$-spin$_2$ dynamics for aligned spins. We further implement these novel PN results in an effective-one-body framework, and demonstrate the improvement in accuracy by comparing against numerical-relativity simulations.


I. INTRODUCTION
The burgeoning field of gravitational-wave (GW) astronomy has already shown its potential to revolutionize our understanding of our universe [1], gravity [2], and the nature of compact objects [3,4], such as black holes (BHs) and neutron stars. The detection of compactbinary GW sources and the accurate inference of their parameters is contingent on having accurate theoretical predictions for their coalescence. As a result of this, a variety of techniques, both analytical and numerical, have been developed to understand the coalescence of binary compact objects, with the final goal of providing faithful waveform models that can be used in GW data analysis.
In parallel to PN formalisms, the small-mass-ratio approximation, based on gravitational self-force (GSF) theory, has also seen rapid development (see Ref. [55] and references therein for a review). As suggested by the name, the expansion parameter in this limit is the mass ratio of the two bodies q = m 1 /m 2 1. The lead-ing order in this approximation is given by the geodesic motion of a test-body in a Schwarzschild or Kerr background. Successive corrections, which can be interpreted as a force moving the body away from geodesic motion, are due to the perturbation of the background sourced by the small body's nonzero stress-energy tensor. This selfforce effect on the motion of a nonspinning body has currently been numerically calculated to first order in q for generic orbits in Kerr spacetime [56]. In a recent breakthrough [57], the second-order-in-q binding energy in a Schwarzschild background has been calculated and compared to predictions from the first law of binary blackhole mechanics [58]. Meanwhile, much activity has led to the analytic calculation at very high PN orders (but at first order in q) of gauge-invariant quantities, such as the Detweiler redshift [59][60][61][62][63][64][65][66][67][68] and the precession frequency [67,[69][70][71][72][73][74], including effects of the smaller body's spin. This has quite naturally led to related activity in confronting and validating the PN and GSF approximations [58,75,76] in the domain which both are valid, i.e., for large orbital separations and small mass ratios, as well as in constructing EOB models based on both approximations [77][78][79][80].
Recently, there has also been rapid advance in understanding and employing post-Minkowskian (PM) techniques, using a weak-field approximation GM/rc 2 1 in a background Minkowski spacetime, with no restriction on the relative velocity of the two bodies [81][82][83][84][85][86]. This approximation most naturally applies to the weakfield scattering of compact objects, in which possibly relativistic velocities can be reached. Recent advances in PM gravity and in our understanding of the scattering of compact objects have been spearheaded by modern on-shell scattering-amplitude techniques, developed originally in the context of quantum particle physics (see, e.g., Ref. [86] and references therein).
Scattering amplitudes were used in Ref. [83] to calculate the nonspinning 2PM (O(G 2 ), one-loop) scatter-ing angle, reproducing with astonishing efficiency the decades-old results of Westpfhal [87,88] obtained by classical methods; an equivalent canonical Hamiltonian at 2PM order was derived from amplitudes in Ref. [84]. The scattering angle plays a key role in PM gravity: it encodes the complete local-in-time conservative dynamics of the system (at least in a perturbative sense) and it can be used to specify a Hamiltonian in a given unique gauge [82], which can in turn be used for unbound as well as bound systems (with potential relevance for improving waveform models [89]); see in particular Refs. [90,91]. In Refs. [85,86], the scattering angle and a corresponding Hamiltonian have been obtained at 3PM (two-loop) order for nonspinning systems, and the results have been confirmed and expounded upon in Refs. [24,[92][93][94].
The PM approximation for two-spinning-body systems was first tackled only very recently, with the SO dynamics at the 1PM and 2PM levels first derived by classical means in Refs. [95,96]. These results have since been confirmed by amplitudes methods in Ref. [97], which also gave the 1PM and 2PM dynamics for the S 1 S 2 sector, rounding out the current state of the art for generic-spin PM results beyond tree level. Several other works have also considered amplitudes methods in relation to spinning two-body systems, also beyond the SO and S 1 S 2 sectors (beyond the dipole level in the bodies' multipole expansions), in particular for special cases such as bodies with black-hole-like spin-induced multipole structure and/or for the aligned-spin configuration (in which the bodies' spins are [anti-]parallel to the orbital angular momentum); see, e.g., [45,98,99] and references reviewed therein.
These works demonstrate that the study of gravitational scattering continues to provide novel results and useful insights on the relativistic two-body problem, with implications for precision gravitational-wave astronomy yet to be explored. A particularly powerful example of such an insight concerns the nontrivially simple dependence of the scattering-angle function on the masses [100] (see also [86,90,101]). This was exploited in Refs. [23,24] to obtain almost all the 5PN dynamics (with the exception of 2 out of 36 coefficients in the EOB Hamiltonian; see also Refs. [18,21]) from first-order self-force calculations (while appropriately dealing with nonlocalin-time tail terms). This approach has also been used in Ref. [26,27] to obtain most of the 6PN dynamics. An extension of this approach to spinning systems was used by the current authors in Ref. [102] to obtain the next-tonext-to-next-to-leading order (N 3 LO) SO PN dynamics.
In this paper, we provide details for the calculation of the N 3 LO-PN SO dynamics presented in Ref. [102], which completes the PN knowledge at 4.5PN order together with the NLO S 3 dynamics from Ref. [41] (see also [45]). Furthermore, we extend our analysis to include a derivation of the N 3 LO S 1 S 2 effects, contributing at 5PN order, for the case of spins aligned with the orbital angular momentum. We note that partial results of the N 3 LO-PN SO and N 3 LO S 1 S 2 dynamics have previously been presented in Refs. [33,37], where all terms at G 4 were calculated within the powerful effective field theory framework using Feynman integral calculus. The latter of these references gives further results for all quadraticin-spin terms at N 3 LO.
Our derivations are organized in the following procedures. 1. We argue that the scattering angle for an alignedspin binary has a simple dependence on the masses (when expressed in terms of appropriate geometrical variables), which extends the result of Ref. [100] for nonspinning binaries. This mass dependence implies that the 4PM part of the scattering angle, which encodes the N 3 LO PN dynamics, is determined by terms up to linear order in the mass ratio. We use analytic results for the test-spin scattering angle to fix all terms at zeroth-order in the mass ratio, leaving the linear terms to be fixed by firstorder GSF results.
2. Assuming the existence of a PN Hamiltonian at the desired 4.5PN SO and 5PN S 1 S 2 orders, and making use of its associated mass-shell constraint with undetermined coefficients, we calculate the scattering angle and match it to the constrained form from step 1. This procedure fixes its lower orders in velocity at 3PM and 4PM orders, leaving but half of the linear-in-mass-ratio coefficients to be determined by GSF calculations. We construct the bound-orbit radial action from the scattering angle (via the Hamiltonian dynamics), noting its simple dependence on the bodies' masses.
3. From the radial action, we calculate the redshift and spin-precession invariants and compare them with GSF results available in the literature to determine the remaining coefficients of the scattering angle. Vital to this step is the first law of spinning binary mechanics [58,103,104], which is used to relate the radial action to the redshift and precession frequency, and for which we herein discuss an extension to arbitrary-mass-ratio aligned-spin eccentric orbits.
(Although we work with aligned spins throughout, we note that the aligned SO result actually fixes the SO Hamiltonian also for precessing spins [102].) The paper is organized as follows. Sections II, III and IV discuss points 1, 2 and 3, respectively. In Sec. V, we implement the new PN results in the scattering angle in an EOB model, and use it to compare our results against NR simulations. We conclude in Sec. VI with a discussion of results and potential future work. Finally, Appendix A contains expressions for tail terms in the radial action, while Appendix B contains explicit expressions for a certain mapping between variables used to connect redshift and precession-invariant results from the radial action to GSF results in the literature, which have been previously erroneously (yet innocuously) reported in the literature.

Notation
We use the metric signature (−, +, +, +), and use units in which the speed of light is c = 1. For a binary of compact objects with masses m 1 and m 2 , we use the following combinations of the masses with m 1 < m 2 . We often make use of the rescaled versions of the canonical spins S 1 and S 2 , i.e., and define the following combinations of spins The relative position and momentum 3-vectors are denoted by r and p, respectively. Using an implicit Euclidean background, it holds that where n = r/r with r = |r|, and L is the orbital angular momentum with magnitude L.

II. THE MASS DEPENDENCE OF THE SCATTERING ANGLE
Here we argue that the structure of the PM expansion, applied to the conservative orbital dynamics of a two-massive-body system, leads to simple constraints on the dependence of the scattering-angle function on the bodies' masses, at fixed geometric quantities characterizing the incoming state. We closely follow the arguments given in Sec. II of Ref. [100] for the nonspinning case, considering only the local-in-time, conservative part of the dynamics, while generalizing to the case of spinning bodies, finally, in the aligned-spin configuration.
The motion of a two-point-mass system (the nonspinning case) is effectively governed by the coupled system of (i) geodesic equations for the worldlines of the two point masses, using the full two-body spacetime metric (with a suitable regularization or renormalization procedure), and (ii) Einstein's equations for the metric, sourced by effective point-mass energy-momentum tensors. In the case of spinning bodies, to dipolar order in the bodies' multipole expansions, the geodesic equations are replaced by the pole-dipole Mathisson-Papapetrou-Dixon (MPD) equations [105][106][107], where, for the ith body (i = 1, 2), p µ i (τ i ) is the linear momentum vector, S µν i (τ i ) is the antisymmetric spin (intrinsic angular momentum) tensor, andẋ µ i (τ i ) is the tangent to the body's worldline x i (τ i ). The constraint (2.1c), the "covariant" or Tulczyjew-Dixon spin supplementary condition [108][109][110][111], combined with (2.1a) and (2.1b), uniquely determines a first-order equation of motion for the worldline,ẋ µ . The corresponding effective energy-momentum tensor, sources Einstein's equations, In the PM scheme, an iterative solution to these equations is obtained as an expansion in G of the worldlines, momenta and spins, and of the metric, where η µν is the Minkowski metric, which we henceforth use instead of the full metric g µν for all 4-vector manipulations (index raising and lowering, dot products and squares of vectors, etc.). At the leading orders in (2.4), given by the solutions to (2.1) with g = η, each body moves inertially in flat spacetime, Here, y µ i are constant displacements from the origin at τ i = 0, and u µ i are constant 4-velocities, with u 2 i = −1, so that τ i are the (Minkowski) proper times, and p 2 i = −m 2 i where m i are the constant rest masses. The zeroth-order spin tensors S µν i0 are also constant, and, being orthogonal to u iµ , have been parametrized in terms of a constant mass-rescaled (Pauli-Lubanski, "covariant") spin vector, with dimensions of length, the magnitude of which would measure the radius of the ring singularity of a corresponding (linearized) Kerr black hole. We identify the zerothorder geometric (mass-independent) quantities, y µ i , u µ i and a µ i , with those characterizing the asymptotic incoming state, along with the masses m 1 and m 2 .
Inserting (2.6) into (2.2) (with g = η) yields the zerothorder stress-energy tensor, which serves as a source for the first-order metric perturbation h 1µν in the linearization of (2.3). The solution for the trace-reversedh µν distance of the field point x from the (zeroth-order, flat geodesic) worldline x i0 = y i + u i τ i in its rest frame, and ∂ µ is the flat covariant derivative. (Note that the result for the first-order field (2.8) would be the same whether we used the physical retarded Green's function or the time-symmetric Green's function, given the nature of the zeroth-order source, constant momentum and spin along a flat-spacetime geodesic.) A key property to be noted here is that h 1 is linear in the masses m i , while having a more intricate dependence on the geometric quantities y µ i , u µ i and a µ i . (It is linear in the spins a µ i here only because we are working to linear order in the spins, to dipolar order in the multipole expansions.) In the next step of the iterative scheme, one uses g = η + h 1 in the bodies' equations of motion (2.1) to solve for the first-order perturbations in (2.4) [for which it is sufficient to integrate the RHSs of (2.1a) and (2.1b) along the zeroth-order motion (2.6), and to regularize by simply dropping the divergent self-field contribution]. Importantly, one finds that x µ i1 , p µ i1 /m i and S µν i1 /m i are each linear functionals of h 1µν (x), and are thus linear in the masses. From Poincaré symmetry, it follows that these results can depend on the positions y i only through the vectorial impact parameter b µ = y µ 1 − y µ 2 , where the y µ i here are chosen along the two zeroth-order worldlines by the conditions u 1 · b = u 2 · b = 0 (at mutual closest approach). For example, the impulse (net change in momentum) for body 1, ∆p µ 1 = Gp µ 11 (τ 1 → ∞) + O(G 2 ), is given by 1 1 Results equivalent to the first two lines of Eq. (2.9) were first is the asymptotic relative Lorentz factor, and Π µ ν = µραβ νργδ u 1α u 2β u γ 1 u δ 2 /(γ 2 − 1) is the projector into the plane orthogonal to both u 1 and u 2 . Here, as below, we work to linear order in each spin, a 1 and a 2 , keeping the cross term. We note again in (2.9) the simple dependence on the masses, with an overall factor of m 1 m 2 , at fixed geometric quantities b µ , u µ i and a µ i . In continuing the iterative PM solution, the O(G n ) terms in the bodies' degrees of freedom (2.4) correct the source (2.2) for the field equation (2.3), determining the O(G n+1 ) metric perturbation in (2.5); the latter, via the bodies' equations of motion (2.1), determines the O(G n+1 ) corrections in (2.4). As in Ref. [100] we are assuming here a systematic use of the time-symmetric Green's function, to pick out the conservative sector of the dynamics. It becomes evident from the structure of these expansions that the O(G n ) metric perturbation h µν n in (2.5) can be expressed as a homogeneous polynomial of degree n in the masses, where the h µν ··· on the RHSs are functions only of the (asymptotic incoming) geometric quantities (y µ i , u µ i , a µ i ) and the field point x. The first line of (2.11) matches (2.8). Similarly, the O(G n ) corrections x µ in , p µ in /m i , S µν in /m i for the body degrees of freedom (2.4) will be homogeneous polynomials of degree n in the masses; this is the crucial point for the following analysis (and for its conceivable extensions beyond the aligned-spin case). The zeroth-order quantities (2.6) are (taken to be) independent of the masses, as is the zeroth-order metric h 0 = η; they, along with the masses, both (i) fully parametrize the asymptotic incoming state and (ii) can be used to parametrize all the higher-order corrections.
Let us now specialize to the case of aligned spins, in which both spin vectors a µ i are (anti-)parallel to the orbital angular momentum, all of which remain constant throughout the scattering, while the orbital motion is derived in Ref. [95], and the last line results from an expansion in spins of the all-orders-in-spin results for black holes from Ref. [112], both references having worked from purely classical considerations; see also [113,114] for derivations from quantum scattering amplitudes.
confined to the fixed plane orthogonal to the angular momenta (just as for the nonspinning case). This entails u 1 · a i = u 2 · a i = 0 and b · a i = 0. Choosingẑ µ (witĥ z 2 = 1) to be the direction of the orbital angular momentum (∝ − µνρσ u ν 1 u ρ 2 b σ ), let us write a µ i = a iẑ µ for the constant rescaled spin vectors (equal to their incoming values), where the scalars a i are positive for spins aligned withẑ µ and negative for anti-aligned. Crucially, in this case, the only nontrivial independent Lorentz-invariant scalars that can be constructed from the vectors u µ i , a µ i and b µ are the magnitude b = (b 2 ) 1/2 of the impact parameter and the two spin lengths a 1 and a 2 , all three with dimensions of length, and the dimensionless Lorentz factor γ = −u 1 · u 2 . Now consider the extension to higher orders in G of the impulse ∆p µ 1 (2.9), which equals −∆p µ 2 (under the conservative dynamics) as the total momentum p µ 1 + p µ 2 is conserved. Its magnitude Q := (∆p 1µ ∆p µ 1 ) 1/2 must be a Lorentz-invariant scalar. In the aligned-spin case, given the previous discussion, and due to Poincaré symmetry and dimensional analysis, it must be a function only of the dimensionless scalar γ and the dimensionlength scalars b, a 1 , a 2 , Gm 1 and Gm 2 . Given also the conclusion from above that, in (2.4) with i = 1, p µ 1n /m 1 is a homogeneous polynomial of degree n in the masses, with the leading n = 1 result seen in (2.9), it follows that the magnitude Q of the impulse must take the following form through fourth order in G (through 4PM order), Furthermore, Q must be invariant under an exchange of the two bodies' identities, (m 1 , a 1 ) ↔ (m 2 , a 2 ). At 1PM order, this tells us that Q 1PM (γ, a 1 /b, a 2 /b) is symmetric under a 1 ↔ a 2 , and thus Q 1PM a2 , so that the third line of (2.12b) in this case is proportional to a 1 + a 2 . Indeed, the explicit expression for Q 1PM is given by the magnitude of the aligned-spin specialization of (2.9) (divided by 2Gm At 2PM order, the 1 ↔ 2 symmetry tells us that each of the two functions in the second line of (2.12a) determines the other, This function, like Q 1PM , is in fact fully determined by the (extended) test-body limit of Q/(m 1 m 2 ) -the limit where one of the masses, say, m 1 , goes to zero, while keeping fixed m 2 , a 2 and a 1 (and γ and b). The result for Q/m 1 in this limit can be consistently determined by solving the pole-dipole MPD equations (2.1) for a spinning test body in a stationary Kerr background; we will present explicit results from this procedure below in terms of the scattering-angle function. This test-body limit, with m 1 → 0, determines all of the functions Q nPM m n−1 2 with no powers of m 1 , for all n, and the 1 ↔ 2 symmetry also tells us that . They are however still constrained by the exchange symmetry as follows. Firstly, which implies that the third line of (2.12b) for Q 3PM m1m2 (like for Q 1PM above) is proportional to a 1 +a 2 . Secondly, so that one of these two functions determines the other.
2 Note that this is the expansion to linear order in the spins of the result (80) from [112] for a two-black-hole system, 13) to all orders in the spin-multipole expansion at 1PM order.
Taking all of these constraints from exchange symmetry, we can eliminate all of the Q's with more m 1 's in the subscript for those with more m 2 's, while those with the same number of m 1 's and m 2 's must be symmetric under a 1 ↔ a 2 . First focusing on the nonspinning (a 0 ) part of (2.12a), this becomes  Remarkably, through 4PM order, this is just linear in the mass ratio ν at fixed M . Precisely the same manipulations go through for the a 1 a 2 terms, replacing a 0 with a 1 a 2 in all the subscripts and with an overall factor of a 1 a 2 /b 2 on the right-hand side.
Next consider just the 1PM and 2PM terms of the SO (a 1 ) part of (2.12a), after accounting for the exchange symmetry in the same way as in the previous We recognize in the second line the following spin combinations often used in the PN and EOB literature, We will find it convenient to rescale each of these by the total rest mass M , defining where b stands for background (or big) and t stands for test (or tiny). The (first) reason for these labels is that, in the extended test-body limit [m 1 → 0 at fixed m 2 (or M ) and fixed a 1 and a 2 ], we see that a b → a 2 becomes the spin-per-mass of the big background object with mass M = m 2 , and a t → a 1 becomes the spinper-mass of the tiny spinning test body with negligible mass (with a further reason explained below). Note that a b + a t = a 1 + a 2 . Now extending (2.22) to 4PM order, from (2.12a) accounting for exchange symmetry, using our new notation, we find , all still functions only of γ. We see that (2.25), like (2.21), is linear in the symmetric mass ratio ν (at fixed M , a b and a t ). Now, just as in Eq. (2.14) of [100] -following from conservation of the total momentum p µ 1 + p µ 2 and simple geometry and kinematics (which is identical for the nonspinning and aligned-spin cases) -the scattering angle χ, by which both bodies are deflected in the system's center-of-mass (cm) frame, is related to the magnitude Q of the impulse by sin where p ∞ (called "P c.m. " by Damour) is the magnitude of the bodies' equal and opposite spatial momenta in the cm frame, at infinity, Here, E is the total energy in the cm frame, (2.28) determined by the asymptotic Lorentz factor γ and the rest masses. Note also the definition of the asymptotic relative velocity v as used e.g. in [45,101,102], We will find it convenient to define yet another variable equivalent to γ or v, namely which, like v 2 , can serve as a PN expansion parameter, and unlike v, is real for both unbound and bound orbits, [the squared momentum per mass of the effective test body], while our p ∞ is Damour's "P c.m. ".) We will also find it convenient to define a notation for the dimensionless ratio Γ (Damour's "h") between the total energy and the total rest mass, with Γ > 1 (γ > 1) for unbound orbits, and Γ < 1 (γ < 1) for bound orbits. Then p ∞ = µγv/Γ = µ √ ε/Γ. With this notation in order, we can take our simplified result for the impulse magnitude Q (2.12a) [namely the sum of (2.21), its analogous a 1 a 2 version, and the SO part (2.25)], insert it into (2.26), and solve for the aligned-spin scattering angle χ. After this process, χ/Γ turns out to be linear in ν in the same way that Q is, thanks to the facts that the sine function is odd in its argument and that Γ 2 is linear in ν. The result can be expressed as follows, where each X ν m G k takes the form with × standing for the "cross term" a 1 a 2 , and with the special constraints All the X's on the right-hand side of (2.33b) are dimensionless and are functions only of the dimensionless ε = γ 2 − 1; they can be expressed in terms of the above Q(γ)'s alone. We see that the 1PM and 2PM terms in (2.33) are independent of the symmetric mass ratio ν and are thus fully preserved in the (extended ) test-body limit ν → 0 (at fixed M , or equivalently m 1 → 0 at fixed M , and at fixed a 1 , a 2 , b and γ), while the 3PM and 4PM terms are linear in ν. This allows us to deduce the complete 1PM and 2PM results for χ/Γ from its test-body limit, and the complete 3PM and 4PM results from first-order self-force (linear-in-mass-ratio) calculations.
The special constraints (2.34) are consequences of the 1 ↔ 2 symmetry, as seen in the G 1 ν 0 and G 3 ν 1 SO terms in (2.25). This is a prediction of the above arguments which our considerations below will be able to test, rather than to rely on. For the case of the G 3 ν 1 SO terms, which we will determine (in a PN expansion) below from matching to first-order self-force calculations, we will allow X 1 3b and X 1 3t to be independent -in fact, X 1 3b will be determined by the redshift invariant in a Kerr background and X 1 3t by the spin-precession invariant in a Schwarzschild background -and we will find from the matching procedure that they are indeed equal through the considered PN orders. The fact that the complete content of Eqs. (2.33) holds through N 2 LO in the PN expansion can be seen in Eqs. (4.32) of Ref. [101].
The ν 0 terms in (2.33) can be determined by solving the MPD equations of motion (2.1) for a spinning (poledipole) test body in a stationary background Kerr spacetime. An integrand for the test-spin-in-Kerr aligned-spin scattering angle function, to all PM orders, was derived in Ref. [115]; see, e.g., their Eq. (66) (which also includes pole-dipole-quadrupole terms for a test black hole). The results of the integration are as follows, to all orders in ε (to all PN orders at each PM order), extending Eq. (5.5) of Ref. [101] to 4PM order in the spin-orbit and bilinear- 3 In Ref. [102], the expression of the result (2.33) for the mass dependence of the scattering angle differed in that (i) we did not pull a factor of 1/ √ ε out of the X's for every factor of 1/b, (ii) we used v instead of ε, and (iii) we used a + and δ a − in place of a b and at, with a ± := a 2 ± a 1 and δ := (m 2 − m 1 )/M ; the equivalence of the two expressions is apparent since (2. 35) in-spin terms. The nonspinning parts are the SO parts are and the bilinear-in-spin parts are The ν 1 terms in (2.33), at 3PM and 4PM orders, can be determined in a PN expansion (here, an expansion in ε) from first-order self-force results (as well as from consistency with lower orders), as we will explicitly demonstrate below for the spin parts. We will use the known nonspinning coefficients through 4PM-3PN order [101], Note that, through 2PM order and up through the SO terms, the first two lines of the right-hand side of (2.33a), with (2.36a) and (2.36b) plugged into the first two lines of (2.33b), correctly give either (i) the aligned-spin scattering angle for a spinning test body with rescaled spin at in a Kerr background with mass M and rescaled spin a b , or (ii) the rescaled aligned-spin scattering angle χ/Γ for the arbitrary-mass two-spinning-body system, using the "spin maps" (2.24); this is a further reason for the labels at and a b . This gives a different "EOB scattering-angle mapping," an alternative to Eq. (3.16) of [101], which produces the 1PM and 2PM SO terms in the two-body scattering angle from its extended test-body limit. [Note however that this different mapping fails at quadratic order in the spins, while Eq. (3.16) of [101] still holds, according to all known results.] noting the transcendental ζ(2) contribution in the last term (the 4PM-3PN term). We will parametrize the SO coefficients as with i = b, t, and the bilinear-in-spin coefficients as We have included all the same powers of ε present in the ν 0 coefficients (2.36), up to the orders in ε which will contribute at the N 3 LO PN level. (We have also factored out γ = √ 1 + ε in the SO terms and π in the 4PM terms, following the patterns at ν 0 .) For these X 1··· kn , which are all pure numbers, k gives the PM order, and n gives the maximum PN order (N n LO) which determines that coefficient. This labeling and the consistency and sufficiency of this ansatz for the scattering angle will become evident in the matching between the scattering angle and a canonical Hamiltonian described in the following section.
Finally, it is important to note that the impact parameter b appearing everywhere in this section is the distance orthogonally separating the two spinning bodies' asymptotic incoming worldlines as defined by the "covariant" or Tulczyjew-Dixon condition [108][109][110][111], Eq. (2.1c) above, for each body-the so-called "proper" or "covariant" impact parameter b ≡ b cov [45,101,116]. This is crucial to the above argument because only with the covariant condition (2.1c) (or something equivalent to it at 0PM order) does it hold that the first-order field (2.8) is linear in the masses. Below, we will also work with the canonical orbital angular momentum L ≡ L can = p ∞ b can , where b can is the impact parameter orthogonally separating the asymptotic incoming worldlines defined by cm-frame Newton-Wigner conditions [117,118] for each body. This coincides with the conserved canonical orbital angular momentum L appearing in a canonical Hamiltonian formulation of aligned-spin two-body dynamics [119,120].
[Note that, for the aligned-spin case, the covariant/Pauli-Lubanski spin vectors m i a µ i used above coincide with the canonical spin vectors S µ i (spatial vectors in the cm frame) which would be associated with the cm-frame Newton-Wigner conditions, and thus so do the alignedspin (signed) magnitudes, S i = m i a i .] As shown in [101,112], the canonical L =: L can is related to the covariant b by Solving this for b, inserting the result into (2.33) [or (2.39)], and re-expanding to bilinear order in the (mass-rescaled) spins a 1 and a 2 , one obtains the final parametrized form for the aligned-spin scattering angle function χ(E, L; m i , a i ) used in the following matching calculations. Let us finally rewrite the scattering angle to include both the ν 0 and ν 1 terms in single coefficients (or which could allow mass-dependence differing from that deduced above), and which would accommodate general quadratic-in-spin terms, with sums over i and j implied, Our prediction for the mass-ratio dependence of the kPM The ν 0 coefficients X 0A k (ε) from the extended test-body limit are given explicitly in (2.36), and the ν 1 coefficients X 1A k (ε) which we will determine from self-force results are parametrized in a PN expansion in (2.37). Note that we will also be able to use the self-force results to test the fact that there are no ν 1 terms at 1PM and 2PM orders in this parametrization of the scattering angle. The fact that there are no ν 2 or higher terms through 4PM order cannot be probed with first-order self-force results, but has already been confirmed by arbitrary-mass PN results through N 2 LO. Our prediction for the mass-dependence will yield new arbitrary-mass results at the N 3 LO PN level once we have fixed the PN expansions of the coefficients X 1A k from first-order self-force calculations.

III. FROM THE UNBOUND SCATTERING ANGLE TO THE BOUND RADIAL ACTION VIA CANONICAL HAMILTONIAN DYNAMICS
Besides the mass dependence of the scattering angle function established in the previous section, and the inputs of test-body results (discussed above) and first-order self-force results (discussed below), the other central ingredient in our derivation is the assumption of the existence of a (local-in-time) canonical Hamiltonian governing the aligned-spin conservative dynamics in the cm frame, for generic (both bound and unbound) orbits, with the Hamiltonian having well-defined (regular, polynomial) PN and PM expansions. Through the desired 4.5PN order in the SO sector and 5PN S 1 S 2 one, we can safely ignore nonlocal-in-time (tail) contributions in the final dynamics/scattering angle. While these do appear at the 4PN level in the nonspinning sector [14] (see e.g., Ref. [121] for a translation into a nonlocal-in-time scattering angle), they only start appearing at 5.5PN order in the spinning one. This can most easily be seen in the first line of Eq.(68a) in Ref. [122], where the linear-in-spin tails are a relative 1.5PN order from the leading quadrupolar contributions to the tail. [As mentioned at the very end of this section, we find it necessary to include tail terms at 4PN order in the nonspinning sector to make contact with available results in the GSF literature.] Our ultimate goal in this section is to take the gaugeinvariant scattering-angle function χ for unbound orbits, parametrized in the previous section, and derive from it a parametrized expression for the gauge-invariant radialaction function I r which characterizes bound orbits, from which we can derive all the bound-orbit gauge invariants to be compared with self-force results in Sec. IV B below.
We do this by passing through the gauge-dependent canonical Hamiltonian dynamics. It is to some extent true that this process (as we implement it here) can be bypassed by using relationships between gauge invariants for unbound and bound orbits found in [91], but not entirely. Those relationships yield We begin in Sec. III A by discussing canonical Hamiltonians for aligned-spin binaries, the resultant equations of motion, and their gauge freedom under canonical transformations, in a PM-PN expansion. We fix a unique gauge by imposing simplifying conditions not on the Hamiltonian function H itself, but on its corresponding "mass-shell constraint" (or "impetus formula" [90]), which is simply a rearrangement of the expression of the Hamiltonian, in which the squared momentum is given as a function of the Hamiltonian H (of the energy E = H). In Sec. III B, we describe how the scatteringangle function can be derived from the canonical massshell constraint, or vice versa (with our gauge-fixing for the mass shell), and derive the explicit relationships between the scattering-angle coefficients and the mass-shell coefficients. Finally, in Sec. III C, we compute the radial action I r , and point out a hidden simplicity in its dependence on the mass ratio, when expressed in terms of appropriate (covariant rather than canonical) variables, which is a simple consequence of the mass dependence of the scattering angle χ and the relationship between χ and I r discovered in [91].
where we note that the canonical orbital angular momentum L is a constant of motion. Such a Hamiltonian is not unique, but is subject to a type of gauge freedom, namely under canonical transformations: diffeomorphisms of the phase space which preserve the canonical form (3.2) of the equations of motion. In a quite general gauge (one which encompasses all gauges encountered in previous PN or PM alignedspin Hamiltonians), the Hamiltonian takes the following form through quadratic order in the spins, through 4PM order, where is the total squared canonical linear momentum. Here, H 0 is the 0PM (free) Hamiltonian, and the functions c k , c i k and c ij k encode respectively the nonspinning, spinorbit, and quadratic-in-spin gravitational couplings at the kPM orders. The c's are assumed to have regular Taylor series around L 2 = 0 and p 2 = 0. We will work here with the standard (gauge) choice for the free Hamiltonian in the cm frame, such that, as r → ∞, the magnitude p 2 of the canonical linear momentum corresponds to the two bodies' physical equal and opposite spatial momenta in the cm frame. The expression (3.3) of the Hamiltonian can be solved, working perturbatively in G, for p 2 (r, E, L; m i , a i ), where E ≡ H(r, p r , L; m i , a i ) is the total energy; one finds where the 0PM part p 2 ∞ is found by (exactly) inverting which we recognize as the same p ∞ from (2.27). The functions f k , f i k and f ij k are determined by (and carry all of the information of) the c ··· k coefficients in the Hamiltonian (3.3). Importantly, the f ··· k functions will have regular limits as γ 2 − 1 = ε → 0 (as p ∞ → 0) and as L 2 → 0, given our assumption that the c ··· k functions were regular as p 2 → 0 and L 2 → 0. The quantities γ, ε and Γ are all defined in terms of the energy E and the rest masses just as in the previous section.
As discussed in Ref. [101] (through N 2 LO in the PN expansion, and as we have explicitly verified through N 3 LO), it is possible to find a perturbative canonical transformation which brings the Hamiltonian (3.3) into a "quasi-isotropic" form, i.e., a form in which the c's depend only p 2 and not on L 2 /r 2 . Furthermore, the freedom in canonical transformations [among Hamiltonians of the form (3.3)] is completely fixed once one imposes this quasi-isotropic-Hamiltonian condition and uniquely specifies a 0PM Hamiltonian H 0 , as we have done in (3.5). For such a quasi-isotropic Hamiltonian, one finds that the corresponding "mass shell constraint," the expression for p 2 (3.6), has nonspinning and SO coefficients f k f i k which are independent of L 2 /r 2 , but its quadratic-in-spin coefficients f ij k have terms at zeroth and first orders in L 2 /r 2 . However, there also exists a different (non-quasi-isotropic) gauge for the Hamiltonian (3.3) (one with L 2 /r 2 terms in c ij k ) such that its mass shell constraint (3.6) is quasi-isotropic, with the f k , f i k and f ij k all depending only on E (and the masses) and not on L 2 /r 2 . Because both the scattering angle and the radial action are more directly related to the f coefficients in the mass shell, we will find it convenient to adopt this quasiisotropic-mass-shell gauge (which is also unique with a given choice for H 0 ), specializing (3.6) to the form Regrouping in terms of powers of r instead of powers of G, we have where we definẽ with f ··· −1 = f ··· 0 = 0, and we need to extend the sum to k = 6 (while dropping the nonspinning f 5 and f 6 ). Our starting point for the following calculations will be this ansatz for the mass shell constraint, which is fully equivalent to an ansatz for a Hamiltonian of the form (3.3) modulo gauge freedom. Our fundamental assumption is the existence of such a canonical Hamiltonian. We will find that the coefficients f ··· k (E; m i ) are uniquely determined by the expansion of the scattering-angle function to kPM order.

B. The scattering angle
As shown in [82], the scattering angle χ(E, L; m i , a i ) for an unbound orbit can be found directly from the canonical mass-shell constraint as follows. The constraint (3.9) can be solved for the radial momentum p r (r, E, L; m i , a i ), and then the scattering angle is given by the integral where r min is the largest real root of p r = 0. In the direct evaluation of this integral, it would matter that thef k in (3.10) depend on L (in the SO terms). But let us define an antiderivative of π + χ with respect to L to be "the unbound radial action," which is essentially a partie finie of the radial action integral for unbound orbits, The eikonal phase [83,97,123,124] is W/ (up to a constant). For the expression of W in terms of thef k , it does not matter that thef k depend on L. That expression will be identical to the L-antiderivative of the nonspinning scattering angle expressed in terms of the nonspinning f k , with f k →f k , so this reduces the evaluation of the integral for the spinning case to the nonspinning problem, using the coefficient mapping (3.10). The results of the nonspinning integral (for χ, from which constructing W is trivial) have been tabulated at high orders, e.g., in [125]. One finds whereχ k are the entries of Table 1 in [125] with f k →f k ; the first few read 14) 12 , The scattering angle χ is then given by with the L-derivative acting also inside thef k in (3.10).
To obtain W or χ through quadratic order in spins and through 4PM order, O(G 4 ), counting both the G k in (3.13) and the 1/G 2 in (3.10), we must include parts of the contributions up tof 6 and up toχ 8 . The resultant explicit expression of the scattering angle χ in terms of the f ··· k coefficients up to 4PM order and quadratic order in spins is We see that the kPM coefficients f ··· k first enter in the G k terms; however, they do not enter those terms at the leading orders in p ∞ (in the PN expansion of each PM coefficient). Recalling that all of the f 's are finite as p ∞ → 0 (ε → 0), we see that, within each set of square brackets multiplying G k , the lowest orders in p ∞ do not depend on f ··· k , rather only on the lower-PM-order f 's (with some exceptions at G 1 and G 2 ). Similarly, for the scattering-angle coefficients at even higher orders in G (some of which will be relevant below), the lower orders in their PN expansions will be determined by coefficients from lower orders in G already appearing here.
This gives the scattering angle χ in terms of the massshell coefficients f k , f i k , f ij k , as an expansion in the canonical orbital angular momentum L. Equating that expression to a parametrization of χ of the form (2.39) in terms of the covariant impact parameter b, using the translation (2.38) while re-expanding in spins, one can solve for the f coefficients in the mass shell in terms of the X coefficients in the scattering angle (or vice versa), order by order in the PM expansion. Recall p ∞ = µ √ ε/Γ. Rewriting ∆L = L − p ∞ b from (2.38) as a sum over (effective) spins, the results for the f 's through 2PM order are as follows: nonspinning, spin-orbit, and quadratic in spin, with symmetrization over i and j understood. These 1PM and 2PM results are exact (to all orders in ε). With our predicted mass-ratio dependence from the previous section, we have, for k = 1, 2, X k (ε, ν) = X 0 k (ε), a i X k i (ε, ν) = a b X 0b k (ε) + a t X 0t k (ε), and a i a j X k ij (ε, ν) = a 1 a 2 X 0× k (ε) + O(a 2 1 , a 2 2 ), all independent of ν, and the X 0··· k (ε) from the extended test-body limit are given explicitly by (2.36). Though it is not immediately obvious here, each of these f 's has a finite limit as ε → 0, as is required by our Hamiltonian ansatz. We will need the expansions of the f ··· 1 up to O(ε 3 ), and of the f ··· 2 up to O(ε 2 ). Along with f ··· 3 up to O(ε 1 ) and f ··· 4 at O(ε 0 ), we will then have a complete mass-shell constraint (3.8) up to N 3 LO in the PN expansion, which could be solved for the corresponding canonical Hamiltonian (3.3).
At 3PM and 4PM orders, one can also solve for the f 's in terms of the X's, obtaining exact expressions analogous to the above. But we will now work in a PN expansion, an expansion in ε, while enforcing our predicted mass-ratio dependence [which (3.18) did not]. For the nonspinning coefficients, using the known results (2.36a) and (2.37a) for the X's, we find (3.19) through the orders that contribute to the N 3 LO PN level.
Here again we note the finite limits as ε → 0. For the spinning contributions, we must enforce that all the f 's have finite limits as ε → 0, which will fix some of the unknown coefficients in our parametrization (2.37) of the ν 1 parts of the scattering angle, or relationships between them, from consistency with the lower-order f 's and X's [recall the discussion following (3.16)]. At the SO level, this determines or constrains the lower-PN-order scattering-angle coefficients, and expressions for f i 3 and f i 4 which are explicitly regular as ε → 0 and depend on the remaining unknowns X 1i 32 , X 1i 33 , and X 1i 43 , with i = b, t, kn , k is the PM order, and n is the relative PN order.

C. The radial action
For a bound orbit (γ 2 − 1 = γ 2 v 2 = ε < 0), the same canonical mass-shell constraint (3.8) governs the motion. The (gauge-dependent) radial momentum function p r (r, E, L; m i , a i ) is still given by (3.24) but now p 2 ∞ = (µ/Γ) 2 ε is negative. As a result, p 2 r (r) has two positive real roots r = r ± between which p 2 r is positive, with r + being the largest real root, and the trajectory oscillates between these radial turning points r ± . The canonical radial action function I r (E, L, m i , a i ) is defined as the integral of p r dr over one period of the radial motion, and it is a gauge-invariant function, from which one can derive several other gauge-invariant functions physically characterizing bound orbits [90,91]. Like the "unbound radial action" W (the L-antiderivative of the scatteringangle χ) (3.12), the bound radial action I r (E, L, m i , a i ) encodes the complete gauge-invariant information content of the canonical Hamiltonian (governing both unbound and bound orbits) (at least up to the N 3 LO PN level) -though in a subtly different way, concerning orders in the PM-PN expansion of I r versus that of W . It was shown in [91] that the periastron-advance angle, Φ = 2π + ∆Φ = −2π∂I r /∂L, the angle swept out by a bound orbit during one period of the radial motion, is related to the scattering angle, π + χ = −2π∂W/∂L, by where the right-hand side requires an analytic continuation from E > M (unbound, for which χ is real) to E < M (bound, for which χ is complex), as detailed below. It follows from a straightforward extension of their argument that a particular L-antiderivative of this relation holds, giving the bound radial action I r in terms of the unbound radial action W , as can also be verified by explicit calculation.
Consider the unbound radial action in the form (3.13), after replacingχ 1 using (3.14) and (3.18a), (3.28) In continuing this from the unbound case, ε > 0, p 2 ∞ > 0, to the bound case, ε < 0, p 2 ∞ < 0, the second term with 1/ √ ε becomes imaginary, as do all of the terms in the sum with k odd, having odd powers of p ∞ = (µ/Γ) √ ε. Note, from (3.14) and (3.10), and from the fact that all of the f 's have regular Taylor series in ε about ε = 0, that all of theχ k are still real for the bound case, and that thẽ χ k are unchanged by (L, a i ) → (−L, −a i ). Thus, plugging the continuation of (3.28), with √ ε = i √ −ε, into (3.27), we see that all of the odd-k terms are canceled; after using ln L − ln(−L) = ln(−1) = −iπ (choosing the branch which yields the physically sensible result), we are left with which is real for bound orbits. Only theχ k with k even (k = 2l) remain, and those with k odd are gone (except forχ 1 ). This may make it seem as though we have lost information in passing from W to I r , but in fact we have not, as long as we are sure to keep all terms in the consistent PN expansion of I r (at least up to the N 3 LO PN level); this is due to relationships between theχ k as discussed below (3.16). As we will make clearer below, the complete PN expansion of I r up to N 3 LO is contained in its PM expansion up to O(G 6 ) for the nonspinning terms and up to O(G 8 ) for the spin-orbit and quadratic-in-spin terms. This can be computed directly from (3.29), recalling that theχ k are the entries of Table 1 of [125] with f k →f k , as in (3.16) above, with thef k given by (3.10). We need again the contributions from To reach the all the G 8 quadratic-in-spin terms, we must take the sum in (3.29) up to l = 6, involving parts ofχ 12 .
This process yields the radial action I r through the N 3 LO PN level as an expansion in the inverse canonical orbital angular momentum L ≡ L can . To express the results of that process, it will be advantageous to use the covariant orbital angular momentum L cov , which we define for the bound-orbit case by L cov := L − ∆L, (3.30) with ∆L(E, a i ) still given by the last two lines of (2.38) or by (3.17), in which we note that everything is still real for bound orbits [unlike in the second line of (2.38), where we would need to continue to imaginary b to keep In fact, the expression of the radial action (mostly) in terms of L cov is simply related to the expression of the scattering angle in terms of L cov , as follows. Taking the form (2.39) for the scattering angle and eliminating b in favor of L cov = (µ/Γ) √ εb, and then using (3.12a), being sure to match up the constant of integration with (3.28), we find where these X k ··· (ε, ν) are precisely the same coefficients from the scattering angle in (3.31). These coefficients up through k = 2l = 4 are those we gave or parametrized above in (2.36) and (2.37), with (2.41). Recollecting them here, while using the constraints (3.20) and (3.22) obtained in matching between the scattering angle and the canonical mass shell, we have the G 2 coefficients which are independent of ν and are known exactly, and the G 4 coefficients which are linear in ν, As mentioned above, for the complete expression of the radial action at the N 3 LO PN level, we need the low orders in the PN expansions of X 6 ··· and (for the spin terms) X 8 ··· . We have obtained these from the procedure to compute the radial action described in the paragraph containing (3.29) and the following paragraph, in which the inputs are the f ··· k up to k = 4 found in the previous subsection, finally changing variables using (3.30) to bring the result into the form (3.31). At G 6 , we find the nonspinning spin-orbit, and bilinear-in-spin, At G 8 , spin-orbit, and bilinear-in-spin Note that the X 6 ··· coefficients in (3.34d) and (3.34e) are exactly quadratic in ν, in spite of the fact that the f 's from which they are constructed, in (3.21) and (3.23), are cubic in ν. Less surprisingly, the X 6 ··· are cubic in ν, and more surprisingly the X 4 ··· are linear in ν and the X 2 ··· are independent of ν. This is all in fact a simple consequence of (i) the link (3.33) between the scatteringangle coefficients X k ··· and the radial-action coefficients, and (ii) the (straight-forward) extension of the predicted mass-ratio dependence (2.41) to kPM order: X k ··· is a polynomial of degree k−1 2 in ν. This is the spinning analog of the "hidden simplicity" of the mass dependence of (the local-in-time part of) the radial action (which is the complete radial action through the N 3 LO PN level) emphasized in Ref. [24]; here in the spin terms, this is crucially dependent on expressing I r in (3.33) in terms the covariant L cov rather than the canonical L.
Finally, we can make the PN order counting explicit by restoring factors of 1/c. Through N 3 LO, (3.33) reads with all the coefficients, to the orders in ε = γ 2 − 1 = O(c −2 ) contributing here at N 3 LO, relative O(c −6 ), given explicitly by (3.34). These depend on the remaining unknowns X 1A kn from the parametrization of the scattering angle, at kPM order and relative nPN order.
In all the above manipulations, it was consistent to keep the nonspinning, spin-orbit, and bilinear-in-spin terms all through the same relative PN orders, here relative 3PN order, N 3 LO. However, in matching to self-force results, due to certain changes of variables discussed below, the treatment of the N 3 LO spin-orbit and bilinearin-spin terms will require the inclusion of the 4PN nonspinning terms. We thus need to add to (3.35) the 4PN nonspinning part of the radial action for bound orbits, which includes contributions from the nonlocal-in-time tail integrals. We present in Appendix A the additional terms at 4PN order, which have been computed from (3.25) applied to the 4PN EOB Hamiltonian derived in [15], valid in an expansion in eccentricity (about the circular orbit limit) to sixth order. Replacing the first two lines of (3.35) with (A1) yields the final form of the radialaction function which we will use to compute the gaugeinvariant quantities to be compared with self-force calculations.

IV. THIRD-SUBLEADING POST-NEWTONIAN SPIN-ORBIT AND SPIN1-SPIN2 COUPLINGS
The remaining unknowns in the parametrization of the scattering-angle function (2.37) can be fixed with available self-force results. The key feature here is the existence of a Hamiltonian/radial action allowing us to connect the scattering-angle to the redshift and spinprecession invariants that, in the small-mass-ratio limit, can be matched to expressions independently calculated in GSF literature. A vital step in this calculation is the first law of BBH mechanics, which we extend to alignedspins and eccentric orbits.

A. The first law of BBH mechanics
The first law of BBH mechanics [58] was first derived for nonspinning point particles in circular orbits in Ref. [58], then generalized to spinning particles on circular orbits in Ref. [103], to nonspinning particles in eccentric orbits in Refs. [104,126], and to precessing eccentric orbits of a point-mass in the small mass-ratio approximation [127]. In the following, we briefly review the arguments leading to these incarnations of the first law for binaries, making explicit how they apply to generic mass-ratio aligned-spin systems on eccentric orbits.
Let us follow Ref. [103] and start out with an action S for the binary, where the compact objects are approximated by effective point-particles moving along worldlines x µ i (τ i ), 2) and the gravitational action S grav is given by the Einstein-Hilbert one with appropriate gauge-fixing and boundary terms. Here Λ cµ i are frame transformations between the coordinate frame and a body-fixed frame (labeled by c = 0, 1, 2, 3) that is Lorentz-orthonormal (Λ ic µ Λ idµ = η cd ). We take τ i to be the (full-metric) proper times from now on. The equations of motion are obtained by varying the action with respect to the dynamical variables X A = {x i , S µν , Λ cµ i , λ µ i , g µν }, leading to Eqs. (2.1)-(2.3), see, e.g., Refs. [109,120]. The dots in Eq. (4.2) represent nonminimal (curvature) couplings to the worldline that may carry undetermined coefficients. These terms also include couplings of quadratic and higher orders in spin related to spin-induced multipole moments of the body [109].
Let us write the action as an integral of a Lagrangian L over coordinate time t as We can vary the Lagrangian L not only with respect to the dynamical variables X A , but also vary certain constants appearing in the action, e.g., the masses C B = {m 1 , m 2 }. Furthermore, taking the dynamical variables X A on-shell (fulfilling their equations of motion) after variation, we arrive at (using summation convention for with a total time derivative (td). Now, if one performs a transformation of the dynamical variables X A → X A , which may depend on the C B , then on-shell it holds Also allowing for changes of the Lagrangian of the form L = L + (td), we arrive at where the subscripts indicate quantities that are kept fixed during differentiation and with . . . an appropriate on-shell averaging that removes the total time derivatives. For generic bound orbits, one can average the conservative motion in Eq. (4.6) over an infinite time in order to remove total time derivatives, which can be traded for a phase-space average in regions where the motion is ergodic; see, e.g., Refs. [127,128]. For the aligned-spin case where the motion is confined to a plane, all oscillatory behavior can be removed by an average over a single orbit [104] (defined as an oscillation cycle of the radial distance r); this is the averaging used in the present paper. Further specializing to circular orbits, the radial distance is constant and hence the average becomes trivial [103]. Finally, note that another benefit of the averaging in Eq. (4.6) is that it helps to make expressions manifestly gauge-invariant [127], which is important when matching PN Hamiltonians to (eccentric-orbit) self-force results.
It is straightforward to generalize the discussion from Lagrangians L to Hamiltonians H . Hamilton's dynamical equations for some pairs of canonical variables (q c , p c ) are equivalently encoded by Hamilton's action principle, Noting that the dynamical variables are now X A = {q c , p c }, and that the kinematic pq-terms in L are independent of the C B , we see that either Lagrangian in Eq. (4.6) can be replaced by minus a Hamiltonian (i.e., it can be applied also to canonical transformations between two Hamiltonians). The rather general on-shell relation (4.6) is interesting on its own, aside from facilitating the derivation of the first law of binary dynamics as demonstrated below.
We are now in a position to elaborate on the redshift variables z i [58,103,104], where the first equality is the definition of z i adopted by us and the second equality is a consequence of the definition of L (4.3) together with the original point-particle action (4.2), dt L ∼ −m i dt dτ i /dt. We note that this relation holds to all orders in spin if the coefficients in the nonminimal couplings (the dots) in Eq. (4.2) are normalized such that no further explicit dependence on the masses m i arises [68]. Now, several nontrivial transformations of the original action (4.1) are performed to arrive at a PN Hamiltonian (see, e.g., Refs. [39,103,126]): a transformation to SO(3)-canonical (Newton-Wigner) variables for the spin degrees of freedom, integrating out the orbital/near-zone metric or tetrad field (calculating the "Fokker action"), reduction of higher-order time derivatives via further variable transformations, a Legendre transform to the Hamiltonian H, specialization to the COM system, and eventually reducing nonlocal-intime tail contributions to local ones. However, all of these transformations fall into the class of transformations (X A , L) → (X A , L ) discussed above, so we may apply Eq. (4.6) (with L → −H) to Eq. (4.8) and conclude that the redshift variables z i can be obtained from a PN Hamiltonian H via Beside the redshift, let us introduce the (averaged) spin precession frequency Ω i as another important observable [103], (4.10) The (instantaneous, directed) precession frequency Ω inst Si can be read off from the equations of motion for the SO(3)-canonical spin vectors S i i generated by the Hamiltonian H, Indeed, this describes a precession of the spin vector; it is straightforward to see that the spin length From now on, as in previous sections, we simplify the discussion to nonprecessing (aligned or anti-aligned) spins, so that Ω inst i S i and d S i /dt = 0. That is, the spin degrees of freedom become nondynamical and can be dropped from the set of dynamical variables. 5 We can now include the spin lengths into our set of constants, We have now arrived at the important Eqs. (4.9) and (4.13) for the (gauge-invariant) observables z i and Ω i , that could be used to relate a PN Hamiltonian H to self-force results [67,129]. But here, for the purpose of matching to self force, we perform a canonical transformation to different phase-space variables that simplify explicit calculations and connects to the radial action introduced above.
As a first step in that direction, we choose the (nonprecessing) motion to be in the equatorial plane θ = π/2, removing the polar angle θ and its canonical conjugate momentum p θ from the phase space; the Hamiltonian is now of the form discussed in Sec. III A. Furthermore, since we consider a system where the Hamilton-Jacobi equation is separable, one can construct a special canonical transformation (for bound orbits) where the constant action variables are the new momenta [130], with the COM orbital angular momentum of the binary p φ ≡ L = const conjugate to the azimuthal angle φ. The advantage of these variables for our purpose is that the averaging . . . over one radial period becomes trivial due to the integral over one radial period in their definition. The canonical conjugates to I r , I φ are the so-called angle variables q r , q φ and evolve linear in time, i.e., their angular frequencies Ω r =q r , Ω φ =q φ are constant [130]; overall Hamilton's equations of motion for the new, canonically transformed, Hamiltonian H (I r , I φ = L; C B ) read Recalling that C B = {m i , S i }, we can apply Eq. (4.6) (with both Lagrangians replaced by Hamiltonians) for the canonical transformation to action-angle variables as well. Equations (4.9) and (4.13) then turn into where the averaging over one radial period is inconsequential and can be dropped. Collecting Eqs. (4.15) and (4.17), we see that the differential of the COM energy E ≡ H can be written as In analogy to the first law of thermodynamics for the differential of the internal energy, this can be called the first law of conservative spinning binary dynamics for nonprecessing bound orbits (covering eccentric orbits and generic mass ratios). It also resembles the first law of BH thermodynamics, which provides a relation for the differential of the Arnowitt-Deser-Misner (ADM) energy dm i of an isolated BH and can be generalized to other compact objects as well [131]. Recall that Eq. (4.18) is valid to all orders in spin, if the coefficients of possible nonminimal coupling terms denoted by dots in Eq. (4.2) are normalized such that no additional dependence on m i arises. It would be interesting to consider these coefficients as part of the constants C B in future work.
Since the fundamental function introduced in the last section that generates observables for bound orbits is the radial action I r (E, L; m i , S i ), we consider the first law (4.18) in the form where we have introduced 20) As a consequence of the first law, we hence obtain Now the redshift variables can be calculated, from a given radial action I r , as the ratio of proper and coordinate times, which manifestly agrees with the (inverse of the) Detweiler-Barack-Sago redshift invariant calculated in GSF literature [59,132]. The spin-precession frequency Ω i is given by Ω i = Φ i /T r from which we obtain the spin-precession invariant [69] (4.24)

B. Comparison with self-force results
Starting from the radial action (3.33), we calculate the redshift z 1 and spin-precession invariants ψ 1 of the small body using Eqs. (4.23) and (4.24). To compare with results available in the literature, we express them in terms of the gauge-invariant variables 6 which are linked to (ε, L) via Eqs. (4.15) and (4.22). The expressions we obtain for z 1 (x, ι) and ψ(x, ι) agree up to N 2 LO with those in Eq. (50) of Ref. [67] and Eq. (83) of Ref. [129]. The full expressions up to N 3 LO are lengthy, which is why we provide them as a Mathematica file in the Supplemental Material [133]. Next, we expand U 1 ≡ z −1 1 and ψ 1 to first order in the mass ratio q, first order in the massive body's spin a 2 , and zeroth order in the spin of the smaller companion a 1 , withâ = a 2 /m 2 . In performing that expansion, we make use of the gauge-independent variables y and λ, which are related to x and ι via To These results can be directly compared with the GSF results in Eq. (4.1) of Ref. [65], Eq. (23) of Ref. [134] and Eq. (20) of Ref. [67] for the redshift, and Eq. (3.33) of Ref. [72] for the precession frequency. At N 2 LO, as expected, our expressions depend on the scatteringangle coefficients. Upon matching these with the abovementioned equations in the literature, we get the following four constraints (at each order in eccentricity): Again, the special constraint (2.34) is satisfied by X 1b 33 and X 1t 33 . Considering the S 1 S 2 dynamics, the relevant constraints can be obtained from the linear-in-spin correction to the spin-precession invariant, which in terms of the remaining unknown coefficients X 1× ij reads Importantly, we have checked that all the above results can be reproduced by starting from a Hamiltonian ansatz (rather than a radial action), constraining it via the mass-ratio dependence of the scattering angle (calculated via (3.11)), and obtaining the redshift and spinprecession invariants through Eqs. (4.9) and (4.13).

V. EFFECTIVE-ONE-BODY HAMILTONIAN AND COMPARISON WITH NUMERICAL RELATIVITY
In this section, we quantify the improvement in accuracy from the new N 3 LO SO and S 1 S 2 corrections using numerical relativity (NR) simulations as means of comparison. We do this using an EOB Hamiltonian, calculated using the scattering angle obtained above, since the resummation of PN results it grants is expected to improve the agreement with NR in the high-frequency regime.
The EOB Hamiltonian is calculated from an effective Hamiltonian H eff via the energy map where Λ = (r 2 + a 2 ) 2 − ∆a 2 with ∆ = r 2 − 2GM r + a 2 . The Kerr spin a is mapped to the binary's spins via a = a 1 + a 2 , and the potentials are taken to be i.e., we factorize the PN corrections to the Kerr potentials. The zero-spin corrections A 0 (r), D 0 (r) and Q 0 (r) are given by Eq. (28) of Ref. [54] and are based on the 4PN nonspinning Hamiltonian derived in Ref. [15]. The SO corrections are encoded in the gyro-gravitomagnetic factors g S , and g S * , while the S 1 S 2 corrections are added through A SS , B SS pr , and Q SS . For those PN corrections, we choose a gauge such that g S , and g S * are independent of L [46,47,135]; we write an ansatz such that, up to N 3 LO, for some unknown coefficients α ij and α * ij . For the S 1 S 2 corrections, A SS and B SS pr start at NLO and are independent of p r , while Q SS starts at NNLO and depends on p 4 r or higher powers of p r , i.e., we use an ansatz of the form To determine those unknowns, we calculate the scattering angle from such an ansatz using Eq. (3.11) (which entails inverting the EOB Hamiltonian for p r in a PN expansion, differentiating with respect to L, and integrating with respect to r). We then match the result of that calculation to the scattering angle calculated in the previous section and solve for the unknown coefficients in Binding energy versus the velocity parameter vω for the S1S2 contribution to the EOB (left panels) and PN-expanded (right panels) binding energies for mass ratios q = 1 (top panels) and q = 1/3 (bottom panels). The NR error is indicated by the shaded regions.
Expressing the PN-expanded Hamiltonian in terms of v ω yields, for the SO part, (5.12) The same steps can be performed numerically to obtain the EOB binding energy without a PN expansion.
To examine the effect of the new N 3 LO terms on the binding energy, we isolate the SO and the S 1 S 2 contributions to the binding energy by combining configurations with different spin orientations (parallel or antiparallel to the orbital angular momentum), as explained in Refs. [138,140]. For the SO contribution, we use (5.13) while for the S 1 S 2 contribution, we use (5.14) In Fig. 1, we plot the SO contribution to the EOB and PN-expanded binding energies versus the velocity parameter v ω for spin magnitudesâ = 0.6. We also plot the NR results by combining the binding energies of configurations with different spins using results from Refs. [138,139]. From the figure, we see that, adding each PN order improves agreement of the EOB binding energy with NR, especially in the high-frequency regime, with better improvement for equal masses than for unequal masses. In contrast, the PN binding energy, plotted using Eq. (5.11), seems not to converge towards NR in the high-frequency regime, with little difference between the N 2 LO and N 3 LO SO orders. Figure 2 shows the S 1 S 2 contribution to the EOB and PN binding energies. As in the SO case, adding the new N 3 LO significantly improves agreement of the EOB binding energy to NR, especially for equal masses, but there is little difference between PN orders for the PN binding energy.
Note that Figs. 1 and 2 should not be interpreted as a direct comparison between PN and EOB dynamics since our results were obtained for simplicity using exact circular-orbits, which leads to a very different behavior than for an inspiraling binary; Refs. [80,137,138], for example, show that EOB results are significantly better than PN when taking into account the binary evolution. Let us also stress that while the EOB and PN curves are based on the same PN information, the EOB Hamiltonian represents a particular resummation of the PN results. We leave the exploration of other resummations and a calibration to NR for future work.

VI. CONCLUSIONS
GW astronomy allows a multitude of applications in fundamental and astrophysics [1][2][3][4] that rely on accurate waveform models for inferring the source parameters. In this paper, we improved the PN description of spinning compact binaries using information from relativistic scattering and self-force theory, which is an extension of the approach introduced and used in Refs. [23,24,100] for the nonspinning case. We started by extending the arguments from Ref. [100] to show that the scattering angle for an aligned-spin binary has a simple dependence on the masses. This allowed us to determine the SO and aligned S 1 S 2 couplings through N 3 LO in a PN expansion using GSF results for the redshift and precession frequency of a small body on an eccentric orbit in a Kerr background. This result is neatly encapsulated in the gauge-invariant aligned-spin scattering-angle function, given explicitly in Eq. (4.38). The derivation presented here provides the full details for the recently reported result at SO level in Ref. [102], while extending the analysis to aligned S 1 S 2 couplings.
Using these new PN results, we calculated the circularorbit binding energy, the EOB gyro-gravitomagnetic factors, and implemented these results in an EOB Hamiltonian. To illustrate the effect of the new N 3 LO terms, we compared the binding energy with NR simulations (see Figs. 1 and 2,) showing an improvement over the N 2 LO. These results could be implemented in state-ofthe-art SEOBNR [48][49][50][51] and TEOBResumS [52,53] waveform models used in LIGO-Virgo searches and inference analyses [4].
While it is arguable whether PM results already provide a useful resummation of the PN ones [89], the present work shows that, with the crucial contribution of GSF theory, advances in PM theory already allow one to advance the PN knowledge in the spin sector. We thus beseech further research to explore synergies between GSF, PM, and PN theory, along the lines of Refs. [23,24,26,27,45,102] and the present paper. One could, for instance, extend the results in this paper to N 3 LO S 2 couplings, i.e. at quadratic order in each spin. This is an important step to complete the aligned-spin 5PN dynamics for BBHs. However, we leave such a calculation for future work, since it would require currently unavailable GSF results.
One can envision further important work at the inter-face between the PM and GSF approximations. With knowledge of first-order GSF theory, one can in principle determine the full 3PM and 4PM scattering angle in a completely independent way from techniques employed, e.g., in Ref. [85]. To this end, one could calculate the PM expansion of GSF gauge-invariant quantities for bound orbits directly (e.g., expansions in u p valid at all orders in the eccentricity e). This enterprise would have to take great care in the inclusion of tail terms in the dynamics, as well as in the analytical continuation of such results to scattering systems. Should these quantities be calculated, one could exploit the method herein presented to fix the 3PM and 4PM scattering angles without further PN re-expansions. Even better would be a direct GSF treatment of scattering orbits and the scattering angle. This is likely to first come in the form of numerical calculations at first-order in the mass ratio. It will however be worth exploring whether "experimental mathematics" techniques can be used to obtain analytic expressions for the 4PM scattering angle by pushing such numerical calculations to extreme precision (see, e.g., Ref. [62] for an example along these lines in the GSF literature).
Finally, we stress that it is paramount to check our results with more established PN calculations (e.g., with the EFT approach, as was done partially at N 3 LO in Refs. [33,37]), as they have been obtained with a sofar completely unexplored method in the spinning sector that is begging to be further scrutinized.
where e is the eccentricity and p the (dimensionless) semilatus rectum. This defines the turning points of the orbit to be at χ = 0, π. Note that here we use p instead of u p ≡ 1/p from the text since it makes the equations below simpler. To determine the constants of motion E, L as functions of (p, e) we setṙ = 0 at the turning points. While these simultaneous equations can be solved fully, we give their expansion in a, which will be sufficient for this work,