Gravitational spin-orbit dynamics at the fifth-and-a-half post-Newtonian order

Accurate waveform models are crucial for gravitational-wave data analysis, and since spin has a significant effect on the binary dynamics, it is important to improve the spin description in these models. In this paper, we derive the spin-orbit (SO) coupling at the fifth-and-a-half post-Newtonian (5.5PN) order. The method we use splits the conservative dynamics into local and nonlocal-in-time parts, then relates the local-in-time part to gravitational self-force results by exploiting the simple mass-ratio dependence of the post-Minkowskian expansion of the scattering angle. We calculate the nonlocal contribution to the 5.5PN SO dynamics to eighth order in the small-eccentricity expansion for bound orbits, and to leading order in the large-eccentricity expansion for unbound orbits. For the local contribution, we obtain all the 5.5PN SO coefficients from first-order self-force results for the redshift and spin-precession invariants, except for one unknown that could be fixed in the future by second-order self-force results. However, by incorporating our 5.5PN results in the effective-one-body formalism and comparing its binding energy to numerical relativity, we find that the remaining unknown has a small effect on the SO dynamics, demonstrating an improvement in accuracy at that order.


I. INTRODUCTION
Gravitational-wave (GW) observations [1][2][3] have improved our understanding of compact binary systems, their properties, and their formation channels [4,5]. A crucial component in searching for GW signals and inferring their parameters is accurate analytical waveform models, in which spin is an important ingredient given its significant effect on the orbital dynamics.
Three main analytical approximation methods exist for describing the dynamics during the inspiral phase: the post-Newtonian (PN), the post-Minkowskian (PM), and the small-mass-ratio (gravitational self-force (GSF)) approximations.
In this paper, we determine the 5.5PN SO coupling for the two-body dynamics, which is the fourth-subleading PN order, except for one coefficient at second order in the mass ratio. Throughout, we perform all calculations for spins aligned, or antialigned, with the direction of the orbital angular momentum. However, the results are valid for precessing spins [181], since at SO level, the spin vector only couples to the angular momentum vector.
The results of this paper and the procedure used can be summarized as follows: 1. In Sec. II, we calculate the nonlocal contribution to the 5.5PN SO Hamiltonian for bound orbits, in a small-eccentricity expansion up to eighth order in eccentricity. We do this for a harmonic-coordinates Hamiltonian, then incorporate those results into the gyro-gravitomagnetic factors in an EOB Hamiltonian.
2. In Sec. III, we determine the local contribution by relating the coefficients of the local Hamiltonian to those of the PM-expanded scattering angle. We then calculate the redshift and spin-precession invariants from the total Hamiltonian, and match their small-mass-ratio expansion to first-order selfforce (1SF) results. This allows us to recover all the coefficients of the local part except for one unknown. However, by computing the EOB binding energy and comparing it to NR, we show that the effect of the remaining unknown on the dynamics is small.
3. In Sec. IV, we complement our results for unbound orbits by calculating the nonlocal part of the gaugeinvariant scattering angle, to leading order in the large-eccentricity expansion. 4. In Sec. V, we provide two gauge-invariant quantities that characterize bound orbits: the radial action as a function of energy and angular momentum, and the circular-orbit binding energy as a function of frequency.
We conclude in Sec. VI with a discussion of the results, and provide in Appendix A a summary of the quasi-Keplerian parametrization at leading SO order. The main results of this paper are provided in the Supplemental Material as a Mathematica file [183].

Notation
We use the metric signature (−, +, +, +), and units in which G = c = 1, but sometimes write them explicitly in PM and PN expansions for clarity.
For a binary with masses m 1 and m 2 , with m 2 ≥ m 1 , and spins S 1 and S 2 , we define the following combinations of the masses: define the mass-rescaled spins the dimensionless spin magnitudes and the spin combinations We use several variables related to the total energy E of the binary system: the binding energyĒ = E − M c 2 , the mass-rescaled energy Γ = E/M , and the effective energy E eff defined by the energy map We also define the asymptotic relative velocity v and Lorentz factor γ via and define the dimensionless energy variable (Note that ε used here is denoted p 2 ∞ in Ref. [178].) The magnitude of the orbital angular momentum is denoted L, and is related to the relative position r, radial momentum p r , and total linear momentum p via We often use dimensionless rescaled quantities, such as The total conservative action at a given PN order can be split into local and nonlocal-in-time parts, such that where the nonlocal part is due to tail effects projected on the conservative dynamics [9,184,185], i.e., radiation emitted at earlier times and backscattered onto the binary. The nonlocal contribution starts at 4PN order, and has been derived for nonspinning binaries up to 6PN order [10,178,179]. In this section, we derive the leadingorder spin contribution to the nonlocal part, which is at 5.5PN order. The nonlocal part of the action can be calculated via the following integral: where the label 'LO' here means that we include the leading-order nonspinning and SO contributions, and where the Hadamard partie finie (Pf) operation is used since the integral is logarithmically divergent at t = t. The time-split (or time-symmetric) GW energy flux F split LO (t, t ) is written in terms of the source multipole moments as [9] F split 3) The mass quadrupole I ij and the current quadrupole J ij (in harmonic coordinates and using the Newton-Wigner spin-supplementary condition [186,187]) are given by [188,189] where the indices in angle brackets denote a symmetric trace-free part. As was shown in Refs. [9,10], the nonlocal part of the action can be written in terms of τ ≡ t − t as Following Ref. [178], we choose the arbitrary length scale s entering the partie finie operation to be the radial distance r between the two bodies in harmonic coordinates. This has the advantage of simplifying the local part by removing its dependence on ln r.
A. Computation of the nonlocal Hamiltonian in a small-eccentricity expansion The integral for the nonlocal Hamiltonian in Eq. (2.5) can be performed in a small-eccentricity expansion using the quasi-Keplerian parametrization [190], which can be expressed, up to 1.5PN order, by the following equations: where a r is the semi-major axis, u the eccentric anomaly, the mean anomaly, n the mean motion (radial angular frequency), K the periastron advance, and (e r , e t , e φ ) the radial, time, and phase eccentricities.
The quasi-Keplerian parametrization was generalized to 3PN order in Ref. [191], and including SO and spinspin contributions in Refs. [192,193]. We summarize in Appendix A 1 the relations between the quantities used in the quasi-Keplerian parametrization and the energy and angular momentum at leading SO order. Using the quasi-Keplerian parametrization, we express the source multipole moments in terms of the variables (a r , e t , t), and expand the moments in eccentricity.
In the center-of-mass frame, the position vectors of the two bodies, x 1 and x 2 , are related to x ≡ x 1 −x 2 via [194] In polar coordinates, with r and φ given by Eqs. (2.6) and (2.8), e r and e φ related to e t via Eqs. (A12) and (A13), whileṙ andφ are given byṙ which are only needed at leading order. We then write the eccentric anomaly u in terms of time t using Kepler's equation (2.7), which has a solution in terms of a Fourier-Bessel series expansion nt + e t sin(nt) + 1 2 e 2 t sin(2nt) + . . . . (2.13) We perform the eccentricity expansion for the nonlocal part up to O(e 8 t ) since it corresponds to an expansion to O(p 8 r ), which is the highest power of p r in the 5.5PN SO local part. However, to simplify the presentation, we write the intermediate steps only expanded to O(e t ).
Plugging the expressions for (r, φ,ṙ,φ) in terms of (a r , e t , t) into the source moments used in the time-split energy flux (2.3) and expanding in eccentricity yields F split LO (t, t + τ ) = ν 2 a 4 r n 6 32 5 cos(2nτ ) + 12 5 e t [9 cos(nt + 3nτ ) + 9 cos(nt − 2nτ ) − cos(nt − nτ ) − cos(nt + 2nτ )] (2.14) the orbit average of which is given by (2.15) In the limit τ = 0, this equation agrees with the eccentricity expansion of the energy flux from Eq. (64) of Ref. [193]. Then, we perform the partie finie operation with time scale 2s/c using Eq. (4.2) of Ref. [9], which reads where the functions A 4PN (e t ) and B 4PN (e t ) in the 4PN part are given in Table I of Ref. [178].

B. Nonlocal part of the EOB Hamiltonian
The (dimensionless) EOB Hamiltonian is given by the energy map where the effective Hamiltonian (2.21) The nonspinning potentials A,D, and Q were obtained at 4PN order in Ref. [10]. The 4.5PN gyrogravitomagnetic factors, g S and g S * , are given by Eq. (5.6) of Ref. [182], and are in a gauge such that they are independent of the angular momentum. Note that the gyro-gravitomagnetic factors are the same for both aligned and precessing spins, since the spin vector only couples to the angular momentum vector at SO level. Hence, even though the calculations are specialized to aligned spins, the final result for the gyrogravitomagnetic factors is valid for precessing spins.
Splitting the potentials A,D, Q into a local and a nonlocal piece, and writing the gyro-gravitomagnetic factors as g S = 2 + · · · + 1 c 8 g 5.5PN,loc yields the following LO nonlocal part of the PN-expanded effective Hamiltonian (2.23) Then, we write the nonlocal piece of the potentials and gyro-gravitomagnetic factors in terms of unknown coefficients, calculate the Delaunay average of H nonloc eff in terms of the EOB coordinates (a r , e t ), and match it to the harmonic-coordinates Hamiltonian from Eq. (2.19). Since harmonic and EOB coordinates agree at leading SO order, no canonical transformation is needed between the two at that order.
This yields the results in In this section, we determine the local part of the Hamiltonian and scattering angle from 1SF results by making use of the simple mass dependence of the PMexpanded scattering angle.

A. Mass dependence of the scattering angle
Based on the structure of the PM expansion, Poincaré symmetry, and dimensional analysis, Ref. [64] (see also Ref. [86]) showed that the magnitude of the impulse (net change in momentum), for nonspinning systems in the center-of-mass frame, has the following dependence on the masses: where each PM order is a homogeneous polynomial in the two masses. For nonspinning bodies, the Q's on the right-hand side are functions only of energy (or velocity v). This mass dependence has been extended in Ref. [182] to include spin, such that where b is the covariant impact parameter defined as the orthogonal distance between the incoming worldlines when using the covariant (Tulczyjew-Dixon) spinsupplementary condition [195,196] [85,86,182,197] for more details.) The scattering angle χ by which the two bodies are deflected in the center-of-mass frame is related to Q via [64] sin where P c.m. is the magnitude of the total linear momentum in the center-of-mass frame and is given by where we recall that Therefore, the scattering angle scaled by E/m 1 m 2 has the same mass dependence as Q. (Equivalently, χ/Γ has the same mass dependence as Q/µ, where Γ ≡ E/M .) For nonspinning binaries, and because of the symmetry under the exchange of the two bodies' labels, the mass dependence of χ/Γ can be written as a polynomial in the symmetric mass ratio ν. This is because any homogeneous polynomial in the masses (m 1 , m 2 ) of degree n can be written as polynomial in ν of degree n/2 . For example, for some mass-independent factors c i . Hence, at each nPM order, χ/Γ is a polynomial in ν of degree (n−1)/2 . When including spin, we also obtain a dependence on the antisymmetric mass ratio δ ≡ (m 2 − m 1 )/M , since where α i are some linear combinations of c i . Thus, we find that the scattering angle, up to 5PM and to linear order in spin, has the following mass dependence: where the X . . .
i are functions only of energy/velocity. Since ν and νδ are of order q when expanded in the mass ratio, their coefficients can be recovered from 1SF results.
This mass-ratio dependence holds for the total (local + nonlocal) scattering angle. However, by choosing the split between the local and nonlocal parts as we did in Sec. II, i.e., by choosing the arbitrary length scale s to be the radial distance r, we get the same mass-ratio dependence for the local part of the 5.5PN SO scattering angle. This is confirmed by the independent calculation of the nonlocal part of the scattering angle in Eq. (4.23) below, which is linear in ν. (In Ref. [178], the authors introduced a 'flexibility' factor in the relation between s and r to ensure that this mass-ratio dependence continues to hold at 5PN order for both the local and nonlocal contributions separately.) Terms independent of ν in the scattering angle can be determined from the scattering angle of a spinning test particle in a Kerr background, which was calculated in Ref. [198]. For a test body with spin s in a Kerr background with spin a, the 5PM test-body scattering angle to all PN orders and to linear order in spins can be obtained by integrating Eq. (65) of Ref. [198], leading to Plugging this into Eq. (3.10) determines all the X i (v) and X δ i (v) functions. Hence, we can write the 5PM SO part of the local scattering angle, expanded to 5.5PN order, as follows: where the X ν ij and X δν ij coefficients are independent of the masses, and can be determined, as explained below, from 1SF results. The coefficient X ν 2 59 could be determined from future second-order self-force results.

B. Relating the local Hamiltonian to the scattering angle
The scattering angle can be calculated from the Hamiltonian by inverting the Hamiltonian and solving for p r (E, L, r), then evaluating the integral where r 0 is the turning point, obtained from the largest root of p r (E, L, r) = 0. E and L represent the physical center-of-mass energy and canonical angular momentum, respectively.
As noted above, we express the scattering angle in terms of the covariant impact parameter b, but use the canonical angular momentum L in the Hamiltonian (corresponding to the Newton-Wigner spin-supplementary condition). The two are related via [85,86] L = L cov + ∆L, which can be used to replace L with b in the scattering angle. We can also replace E with v using Eq. (3.5).
Starting from the 4.5PN SO Hamiltonian, as given by Eq. (5.6) of Ref. [182], determines all the unknown coefficients in the scattering angle in Eq. (3.12) up to that order. Writing an ansatz for the local 5.5PN part in terms of unknown coefficients, such as calculating the scattering angle, and matching to Eq. (3.12) allows us to relate the 10 unknowns in that ansatz to the 6 unknowns in the scattering angle at that order. This leads to where we switched to dimensionless variables. We see that the 5 unknowns (X ν 39 , X ν 49 , X δν 49 , X ν 59 , X δν 59 ) from the scattering angle only appear in the linear-in-ν coefficients of the gyro-gravitomagnetic factors up to order p 4 r , while the unknown X ν 2 59 only appears in the quadratic-in-ν coefficients of the circular-orbit (1/r 4 ) part. All other coefficients have been determined, due to the structure of the PM-expanded scattering angle, and from lower-order and test-body results.

C. Redshift and precession frequency
To determine the linear-in-ν coefficients in the local Hamiltonian from 1SF results, we calculate the redshift and spin-precession invariants from the total (local + nonlocal) Hamiltonian, since GSF calculations do not differentiate between the two, then match their smallmass-ratio expansion to 1SF expressions known in the literature.
An important step in this calculation is the first law of binary mechanics, which was derived for nonspinning particles in circular orbits in Ref. [199], generalized to spinning particles in circular orbits in Ref. [200], to nonspinning particles in eccentric orbits in Refs. [201,202], and to spinning particles in eccentric orbits in Ref. [182]. It reads dE = Ω r dI r + Ω φ dL + i (z i dm i + Ω Si dS i ), (3.17) where Ω r and Ω φ are the radial and azimuthal frequencies, I r is the radial action, z i is the redshift, and Ω Si is the spin-precession frequency.
The orbit-averaged redshift is a gauge-invariant quantity that can be calculated from the Hamiltonian using where T r is the radial period. The spin-precession frequency Ω S1 and spin-precession invariant ψ 1 are given by In evaluating these integrals, we follow Refs. [134,141] in using the Keplerian parametrization for the radial variable , (3.20) where u p is the inverse of the semi-latus rectum, e is the eccentricity, and ξ is the relativistic anomaly. The radial and azimuthal periods are calculated from the Hamiltonian using Performing the above steps yields the redshift and spin-precession invariants in terms of the gaugedependent u p and e, i.e., z 1 (u p , e) and ψ 1 (u p , e). We then express them in terms of the gauge-independent variables where k ≡ T φ /(2π) − 1 is the fractional periastron advance. The expressions we obtain for z 1 (x, ι) and ψ(x, ι) agree up to 3.5PN order with those in Eq. (50) of Ref. [134] and Eq. (83) of Ref. [141], respectively. Note that the denominator of ι in Eq. (3.23) is of order 1PN, which effectively scales down the PN ordering such that, to obtain the spin-precession invariant at fourth-subleading PN order, we need to include the 5PN nonspinning part of the Hamiltonian, which is given in Refs. [177,178].

D. Comparison with self-force results
Next, we expand the redshift z 1 (x, ι) and spinprecession invariant ψ 1 (x, ι) to first order in the mass ratio q, first order in the massive body's spin a 2 ≡ a, and zeroth order in the spin of the smaller companion a 1 . In doing so, we make use of another set of variables (y, λ), defined by where the mass ratio q = m 1 /m 2 . Schematically, those expansions have the following dependence on the scattering-angle unknowns: which can be seen from the structure of the scattering angle in Eq. (3.12). In those expressions, the O(a) part of the redshift depends on the unknown X ν 39 and the difference of the two pairs of unknowns (X ν 49 , X δν 49 ) and (X ν 59 , X δν 59 ), while the spin-precession invariant depends on X ν 39 and the sum of the two pairs of unknowns. This means that solving for X ν 39 requires 1SF result for either z 1 or ψ 1 , while solving for the other unknowns requires both.
Hence, to solve for all five unknowns, we need at least three (or two) orders in eccentricity in the redshift, at first order in the Kerr spin, and two (or three) orders in eccentricity in the spin-precession invariant, at zeroth order in both spins. Equivalently, instead of the spinprecession invariant, one could use the redshift at linear order in the spin of the smaller body a 1 , but that is known from 1SF results for circular orbits only [133]. Incidentally, the available 1SF results are just enough to solve for the five unknowns, since the redshift is known to O(e 4 ) [132,134,203] and the spin-precession invariant to O(e 2 ) [139].
The last unknown X ν 2 59 in the 5.5PN scattering angle appears in both the redshift and spin-precession invariants at second order in the mass ratio, thus requiring second-order self-force results for circular orbits.
To compare z 1 (y, λ) and ψ 1 (y, λ) with GSF results, we write them in terms of the Kerr background values of the variables (y, λ) expressed in terms of (u p , e). The relations between the two sets of variables are explained in detail in Appendix B of Ref. [182], and we just need to append to Eqs. (B16)-(B20) there the following PN terms y(u p , e) = y 0 (u p , e) + a y a (u p , e) + O(a 2 ), (3.26) λ(u p , e) = λ 0 (u p , e) + a λ a (u p , e) + O(a 2 ), We obtain the following 1SF part of the inverse redshift U 1 ≡ 1/z 1 and spin-precession invariant ψ 1

E. Local scattering angle and Hamiltonian
Inserting the solution from Eq. (3.30) into the scattering angle in Eq. (3.12), yields For the gyro-gravitomagnetic factors, which are one of the main results of this paper, substituting the solution (3.30) in Eq. (3.16) yields the following local part:

F. Comparison with numerical relativity
To quantify the effect of the 5.5PN SO part on the dynamics, and that of the remaining unknown coefficient X ν 2 59 , we compare the binding energy calculated from the EOB Hamiltonian to NR. The binding energy provides a good diagnostic for the conservative dynamics of the binary system [174,205,206], and can be calculated from accurate NR simulations by subtracting the radiated energy E rad from the ADM energy E ADM at the beginning of the simulation [207], i.e., To isolate the SO contributionĒ SO to the binding energy, we combine configurations with different spin orientations (parallel or anti-parallel to the orbital angular momentum), as explained in Refs. [208,209]. One possibility is to usē where χ here is the magnitude of the dimensionless spin. This relation subtracts the nonspinning and spin-spin parts, with corrections remaining at order χ 3 , which provides a good approximation since the spin-cubed contribution to the binding energy is about an order of magnitude smaller than the SO contribution, as was shown in Ref. [209]. We calculate the binding energy for circular orbits from the EOB Hamiltonian usingĒ EOB = H EOB − M while neglecting radiation reaction effects, which implies thatĒ EOB is not expected to agree well withĒ NR near the end of the inspiral. We set p r = 0 in the Hamiltonian and numerically solveṗ r = 0 = −∂H/∂r for the angular momentum L at different orbital separations. Then, we plotĒ versus the dimensionless parameter v Ω ≡ (M Ω) 1/3 , (3.36) where the orbital frequency Ω = ∂H/∂L. Finally, we compare the EOB binding energy to NR data for the binding energy that were extracted in Ref. [209] from the Simulating eXtreme Spacetimes (SXS) catalog [210,211]. In particular, we use the simulations with SXS ID 0228 and 0215 for q = 1, 0291 and 0264 for q = 1/3, all with spin magnitudes χ = 0.6 aligned and antialigned. The numerical error in these simulations is significantly smaller than the SO contribution to the binding energy. In Fig. 1, we plot the relative difference in the SO con-tributionĒ SO between EOB and NR for two mass ratios, q = 1 and q = 1/3, as a function of v Ω up to v Ω = 0.38, which corresponds to about an orbit before merger. We see that the inclusion of the 5.5PN SO part (with the remaining unknown X ν 2 59 = 0) provides an improvement over 4.5PN, but the difference is smaller than that between 3.5PN and 4.5PN. In addition, since the remaining unknown X ν 2 59 is expected to be about O(10 2 ), based on the other coefficients in the scattering angle, we plotted the energy for X ν 2 59 = 500 and X ν 2 59 = −500, demonstrating that the effect of that unknown is less than the difference between 4.5PN and 5.5PN, with decreasing effect for small mass ratios.

IV. NONLOCAL 5.5PN SO SCATTERING ANGLE
The local part of the Hamiltonian and scattering angle calculated in the previous section is valid for both bound and unbound orbits. However, the nonlocal part of the Hamiltonian from Sec. II is only valid for bound orbits since it was calculated in a small-eccentricity expansion. In this section, we complement these results by calculating the nonlocal part for unbound orbits in a large-eccentricity (or large-angular-momentum) expansion.
The nonlocal part of the 4PN scattering angle was first computed in Ref. [212], in both the time and frequency domains, at leading order in the large-eccentricity expansion. This was extended in Ref. [178] at 5PN at leading order in eccentricity, and in Ref. [179] at 6PN to next-to-next-to-leading order in eccentricity. In addition, Refs. [213,214] recovered analytical expressions for the nonlocal scattering angle by using high-precision arithmetic methods.

A. Radial action
The radial action function contains the same gaugeinvariant information as the Hamiltonian, and from it several other functions can be derived that describe bound orbits, such as the periastron advance, which can be directly related to the scattering angle via analytic continuation [77,88]. This means that the entire calculation in Sec. III could be performed using the radial action instead of the Hamiltonian, as was done in Ref. [182].
The radial action is defined by the integral and we split it into a local contribution and a nonlocal one, such that I r = I loc r + I nonloc r .

(5.2)
We calculate the local part from the local EOB Hamiltonian, i.e., Eq. (2.20) with the nonlocal parts of the potentials and gyro-gravitomagnetic factors set to zero. We invert the local Hamiltonian iteratively to obtain p r (ε, L, r) in a PN expansion, where we recall that with ε < 0, γ < 1 for bound orbits. Then, we integrate where r ± are the zeros of the Newtonian-order p (0) r = ε + 2/r − L 2 /r 2 , which are given by It is convenient to express the radial action in terms of the covariant angular momentum L cov = L − ∆L, with ∆L given by Eq. (3.14), since it can then be directly related to the coefficients of the scattering angle, as discussed in Ref. [182], and leads to slightly simpler coefficients for the SO part. We obtain for the local part where each term starts at a given PN order, with 0.5PN order corresponding to each power in 1/L. Also, as noted in Ref. [178], when the radial action is written in this form, in terms of Γ, the coefficients I (s) 2n+1 become simple polynomials in ν of degree n .
The coefficients I n for the nonspinning local radial action up to 5PN order are given by Eq. (13.20) of Ref. [178]. The SO coefficients I s n were derived in Ref. [182] to the 4.5PN order, but we list them here for completeness. The coefficients I 0 , I 1 , I s 2 are exact, and are given by , where a b ≡ S/M, a t ≡ S * /M . The other SO coefficients, up to 5.5PN, read The nonlocal part can be calculated similarly by starting from the total Hamiltonian, expanding Eq. (5.1) in eccentricity, then subtracting the local part. Alternatively, it can be calculated directly from the nonlocal Hamiltonian via [179] I nonloc where Ω r = 2π/T r is the radial frequency given by Eq. (A8).
The nonlocal Hamiltonian H nonloc in Eq. (2.19) is expressed in terms of (e t , a r ), but we can use Eqs. (A6) and (A12) to obtain I nonloc r (E, L), i.e., as a function of energy and angular momentum. Then, we replace E with (e t , L) using Eq. (A10), expand in eccentricity to O(e 8 t ), and revert back to (E, L). This way, we obtain an expression for I nonloc r in powers of 1/L that is valid to eighth order in eccentricity, and in which each ε n contributes up to order e 2n .
The result for the 4PN  Here, we calculate the gauge-invariant binding energȳ E analytically in a PN expansion, as opposed to the numerical calculation in Sec. III F for the EOB binding energy.
For circular orbits and aligned spins,Ē can be calculated from the Hamiltonian (2.20) by setting p r = 0 and perturbatively solvingṗ r = 0 = −∂H/∂r for the an-gular momentum L(r). Then, solving Ω = ∂H/∂L for r(Ω), and substituting in the Hamiltonian yieldsĒ as a function of the orbital frequency. It is convenient to expressĒ in terms of the dimensionless frequency parameter v Ω ≡ (M Ω φ ) 1/3 .
The nonspinning 4PN binding energy is given by Eq. (5.5) of Ref. [9], and the 4.5PN SO part is given by Eq. (5.11) of Ref. [182]. We obtain for the 5

VI. CONCLUSIONS
Improving the spin description in waveform models is crucial for GW observations with the continually increasing sensitivities of the Advanced LIGO, Virgo, and KAGRA detectors [216], and for future GW detectors, such as the Laser Interferometer Space Antenna (LISA) [217], the Einstein Telescope (ET) [218], the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [219], and Cosmic Explorer (CE) [220]. More accurate waveform models can lead to better estimates for the spins of binary systems, and for the orthogonal component of spin in precessing systems, which helps in identifying their formation channels [4,5]. For this purpose, we extended in this paper the SO coupling to the 5.5PN level.
We employed an approach [177,182] that combines several analytical approximation methods to obtain arbitrary-mass-ratio PN results from first-order self-force results. We computed the nonlocal-in-time contribution to the dynamics for bound orbits in a small-eccentricity expansion, Eq. (2.24), and for unbound motion in a largeeccentricity expansion, Eq. (4.23). To our knowledge, this is the first time that nonlocal contributions to the conservative dynamics have been computed in the spin sector. For the local-in-time contribution, we exploited the simple mass-ratio dependence of the PM-expanded scattering angle and related the Hamiltonian coefficients to those of the scattering angle. This allowed us to determine all the unknowns at that order from first-order self-force results, except for one unknown at second order in the mass ratio, see Eqs. (3.31)-(3.33). We also provided the radial action, in Sec. V A, and the circular-orbit binding energy, in Eq. (5.11), as two important gaugeinvariant quantities for bound orbits. We stress again that, although all calculations in this paper were performed for aligned spins, the SO coupling is applicable for generic precessing spins.
The local part of the 5.5PN SO coupling still has an unknown coefficient, but as we showed in Fig. 1, its effect on the dynamics is smaller than the difference between the 4.5 and 5.5PN orders. Determining that unknown could be done through targeted PN calculations, as was illustrated in Ref. [98], in which the authors related the two missing coefficients at 5PN order to coefficients that can be calculated from an EFT approach. Alternatively, one could use analytical second -order self-force results, which might become available in the near future, given the recent work on numerically computing the binding energy and energy flux [142,143]. Until then, one could still use the partial 5.5PN SO results in EOB waveform models complemented by NR calibration. Such an implementation would be straightforward, since we obtained the gyro-gravitomagnetic factors that enter directly into the SEOBNR [162][163][164] and TEOBResumS [150,166,169] waveform models, and less directly in the IMRPhenom models [221][222][223][224], which are used in GW analyses.

ACKNOWLEDGMENTS
I am grateful to Alessandra Buonanno, Jan Steinhoff, and Justin Vines for fruitful discussions and for their invaluable feedback on earlier drafts of this paper. I also thank Sergei Ossokine for providing NR data for the binding energy, and thank the anonymous referee for useful suggestions.

Elliptic orbits
For a binary in a bound orbit in the orbital plane, and using polar coordinates (r, φ), the quasi-Keplerian parametrization [190], up to 1.5PN, reads r = a r (1 − e r cos u), (A1) ≡ nt = u − e t sin u, (A2) where a r is the semi-major axis, u is the eccentric anomaly, is the mean anomaly, n is the mean motion (radial angular frequency), K is the periastron advance, and (e r , e t , e φ ) are the radial, time and phase eccentricities. Spin was included in the quasi-Keplerian parametrization in Refs. [192,193]. (See Fig. 2 of Ref. [193] for a geometric picture for some of these quantities.) The (dimensionless) harmonic-coordinates Hamiltonian with LO SO reads By inserting r = a r (1 − e r cos u) into the Hamiltonian at periastron (u = 0) and apastron (u = π), one can solve for the energy and angular momentum (with p r = 0) as a function of a r and e r , i.e., , whereĒ ≡ E − 1/ν < 0 is the dimensionless binding energy, which is negative for bound orbits, and we only include the LO nonspinning and SO terms. These expansions can be inverted to obtain e r (Ē, L) and a r (Ē, L) leading to The radial period T r and periastron advance K can be calculated from the integrals where r p and r a are the periastron and apastron separations calculated from the solution of p r = 0, which yields the PN expansion The three eccentricities (e r , e t , e φ ) agree at LO, and can be related to each other, and to the energy and angular momentum, via
Inverting these expansions, we obtain and e r (Ē, L) the same as in Eq. (A6).
The mean motionn and periastron advance K are given byn The eccentricities e t and e φ are given in terms of energy and angular momentum by Eqs. (A10) and (A11), respectively.