{'fraction_non_alphanumeric': 0.04545366627337963, 'fraction_numerical': 0.01870370012185451, 'mean_word_length': 4.131202858277094, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Next-to-next-to-leading order spin-orbit effects in the equations of motion of compact binary systems

We compute next-to-next-to-leading order spin contributions to the post-Newtonian equations of motion for binaries of compact objects, such as black holes or neutron stars. For maximally spinning black holes, those contributions are of third-and-a-half post-Newtonian (3.5PN) order, improving our knowledge of the equations of motion, already known for non-spinning objects up to this order. Building on previous work, we represent the rotation of the two bodies using a pole-dipole matter stress-energy tensor, and iterate Einstein's field equations for a set of potentials parametrizing the metric in harmonic coordinates. Checks of the result include the existence of a conserved energy, the approximate global Lorentz invariance of the equations of motion in harmonic coordinates, and the recovery of the motion of a spinning object on a Kerr background in the test-mass limit. We verified the existence of a contact transformation, together with a redefinition of the spin variables that makes our result equivalent to a previously published reduced Hamiltonian, obtained from the Arnowitt-Deser-Misner (ADM) formalism.


I. INTRODUCTION
Gravitational wave detectors will soon enter a new era, with the advanced versions of Virgo and LIGO expected to start operating around 2015, and the construction of KAGRA [1][2][3]. Primary targets for these laser ground-based interferometers, or for possible future space-based interferometers, are inspiralling binaries of compact objects -neutron stars and/or black holes -, which emit gravitational radiation during their late inspiral, merger, and ringdown phases. Hunting for the faint signal and best separating it from the noise using matched-filtering techniques (e.g. [4,5]) requires a very precise modelling of the expected waveform, which the post-Newtonian (hereafter PN) approximation aims at providing for the inspiral phase [6].
Up to now, such high-precision PN templates are available for the non-spinning case up to the 3.5PN order, i.e. the order 1/c 7 in the formal expansion in powers of 1/c 2 (with c being the speed of light). However, observational evidence points toward the existence of fast-rotating black holes, endowed with significant angular momentum (or spin), for stellar size [7][8][9][10][11] as well as for supermassive black holes [12][13][14].
It is thus important to complete our knowledge of the PN predictions for the gravitational waves emitted by such systems by including the spins of each of the two bodies. This implies first taking those spins into account in the dynamics of the binary. For maximally spinning objects, the leading-order linear-in-spin (which we shall call spin-orbit) effect arises at 1.5PN order when regarding the spin as a 0.5PN quantity [see Eq. (1.1) below]. The next-to-leading and next-to-next-to-leading contributions show up at respectively 2.5PN and 3.5PN orders. The effect of spacetime curvature on the motion of a spinning test particle was obtained in the seminal work by Papapetrou [15][16][17], after an earlier derivation by Mathisson (republished in [18]). Barker and O'Connell [19,20] determined the leading-order spin-orbit and spinspin effects in the two-body dynamics. More recently, Kidder, Will and Wiseman [21,22] computed the corresponding contributions to the radiation field and, most importantly, to the orbital phase of the binary, to which the templates are crucially sensitive. Effective field theory methods [23] were also used to rederive the leading-order spin-orbit and spin-spin contributions to the dynamics [24]. The problem in the limit of a spinning test particle on a Kerr background has been addressed in [25,26] (see also [27] for a Hamiltonian model neglecting the gravitation damping force).
By means of a PN iteration of the Einstein field equations in harmonic coordinates, the equation of motion at the next-to-leading order was first investigated by Tagoshi, Ohashi and Owen [28,29]. It was then confirmed and completed (as well as extended to the radiation field) by Blanchet, Buonanno and Faye [30,31]. The results for the evolution equations were retrieved by two independent calculations, using a Hamiltonian approach in ADM coordinates [32] on the one hand, and using effective field theory methods [33,34] on the other hand. The ADM computation was later generalized to the many-body problem [35] and extended to the next-to-leading order spin1-spin2 and spin square (e.g. spin1-spin1) effects in Refs. [36,37]. Next-to-next-to-leading order spin-orbit effects in the Hamiltonian of the binary were first computed in [38], resorting to the ADM scheme adapted to matter sources composed of spinning point particles [39]. In a subsequent work, spin1-spin2 interactions terms were also added at the 4PN order [40]. Effective field theory methods progressed concurrently by computing the 3PN spin1-spin2 and spin1-spin1 contributions [41][42][43], and the 4PN spin1-spin2 interactions [44].
Our aim in this work is to extend the approach of Ref. [30], 1 based on the PN expansion of the metric in harmonic coordinates, in order to compute the spin-orbit next-to-next-toleading order 3.5PN contributions to the equations of motion. We shall thus provide an essential validation of the result obtained by the ADM method [38]. Our final purpose is to derive the spin-orbit terms in the emitted GW energy flux and in the phase of the binary at the same approximation, extending therefore Ref. [31] to next-to-next-to-leading order, so as to provide more accurate PN predictions to the template-based data analysis.
We adopt the same convention for the PN order counting of the spins as in Paper I; namely, we redefine our spin variable with respect to the "true" spin angular momentum following : S ≡ cS true = Gm 2 body χ , (1.1) where m body is the mass and χ the dimensionless spin parameter (G denotes Newton's constant). For a maximally rotating compact object, we have χ ∼ 1 and S ∼ Gm 2 body , so that our spin variable can be counted as of Newtonian order in this case. For slowly rotating objects however, we have χ ∼ v surf /c where v surf is the rotation velocity of the body surface, which implies that S acquires an additional factor 1/c. In that case, spin-orbit contributions effectively appear at 2PN order, and spin-spin ones at 3PN only. Our counting will always assume rapid rotation. The leading-order spin-orbit terms will thus carry an explicit factor 1/c 3 ; in the present article we shall obtain the next-to-leading and next-to-next-to-leading spin-orbit corrections ∼ 1/c 5 and ∼ 1/c 7 .
The paper is organized as follows. In Sec. II, we review the pole-dipole formalism for modelling compact objects as point particles with spins. In Sec. III, we describe the parametrization of the PN metric in terms of a set of elementary potentials and present the equations of motion as deduced generically from metric components. The core of our work consists of the computation of the required potentials in Secs. IV and V, first using Hadamard's regularization (as well as the "pure Hadamard-Schwartz" prescription [45]), and next completing some of our calculations by means of the more powerful dimensional regularization. Sec. VI contains our results for the precession equations and the acceleration. In Sec. VII, we explain the various checks that have been performed to validate them: construction of a conserved energy, check of the global Lorentz invariance, test-mass limit, and recovery of the results of Ref. [38] in the ADM Hamiltonian formulation. We conclude shortly in Sec. VIII.

II. EFFECTIVE POLE-DIPOLE FORMALISM
The starting point of our calculations is the model of pole-dipole particles developed by Mathisson [18], Papapetrou [16,17], Tulczyjew [46,47] and generalized by Dixon [48][49][50] and Bailey & Israel [51]. This model allows an effective description of the dynamics of bodies that accounts for their spin angular momentum by means of singular Dirac deltafunctions, making analytical computations tractable. We use the same version of the model as in Paper I, staying linear in the spins, but we summarize here the main formulae for completeness and refer to Paper I for more details.
The central assumption is that each particle is described by a stress-energy tensor made of two parts, a monopolar one and a dipolar one: T µν = T µν M + T µν D . These are built with respectively a Dirac delta function and a gradient of a delta function, integrated over the world line of the particle, according to: . (2.1b) Here τ is the proper time measured along the world line, described itself by the particle position y µ (τ ); δ (4) denotes the four-dimensional Dirac delta function, g stands for the determinant of the spacetime metric g µν ; u µ = dy µ /(cdτ ) is the four-velocity of the particle (satisfying u µ u µ = −1), p µ its linear momentum, and S µν an antisymmetric tensor that represents the spin of the particle. The linear momentum p µ and spin tensor S µν obey the Mathisson-Papapetrou equations of evolution: with D/dτ ≡ cu ν ∇ ν , and R µ νρσ denoting the Riemann tensor. As is well known, a choice of supplementary spin condition (thereafter SSC) is necessary, in order to obtain the correct number of degrees of freedom (see [52] for a summary of the various choices in use in the litterature). We adopt the following covariant SSC: With this relation in hand, one can combine the equations (2.2) to obtain the link between the four-velocity u µ and the linear momentum p µ ; hence we deduce the conservation laws: for the mass and for the magnitude of the spin defined by m 2 c 2 = −p µ p µ and S 2 = S µν S µν /2 respectively, and thus conserved along the particle's world line. We now restrict the evolution equations themselves to linear order in spins, neglecting any term O(S 2 ), quadratic or of higher order. One can readily check the proportionality relation between u µ and p µ in that approximation: (2.5) In particular, the SSC now reads S µν u ν = O(S 3 ). Notice that the monopolar part of the stress-energy tensor has the same form as the one of an ordinary point-like particle, but it does depend on the spin implicitly through the metric. Equations (2.2) become: We see that in the linear approximation, the spin tensor is parallely transported along the motion, which is actually non-geodesic due to the coupling of the spin tensor with the Riemannian curvature. In the following, we shall use the 3-dimensional form of the energy-momentum tensor (2.1). Using a 3 + 1 splitting of spacetime, the particle's position and coordinate velocity are denoted y µ = (c t, y(t)) and v µ (t) = (c, v(t)) (where v µ = cu µ /u 0 , with u 0 = 1/ −g ρσ v ρ v σ /c 2 ), and the spin tensor S µν (t) is considered a function of time. From now on, we use boldface letters to denote three-dimensional vectors. We have is the three-dimensional Dirac delta function. Expliciting the covariant derivative in the dipolar stress-energy tensor, we obtain an alternative form in terms of ordinary Christoffel symbols: Note that in (2.7)-(2.8) the factor √ −g has the generic field point (c t, x) for argument.
Because of the SSC, we may work only with the spatial components S ij of the spin tensor, and eliminate the S 0i components according to: (2.9) In Paper I, a spin vector was used instead of a spin tensor. First, a spin covector was defined, consistently with the SSC S µν u ν = O(S 3 ), by: where ε µνρσ is the Levi-Civita tensor such that ε 0123 = −1/ √ −g. Then, a spatial spin vector S i FBB was constructed as S i FBB = γ ij S j , with γ ij being the inverse of the spatial part of the metric, i.e. γ ik g kj = δ i j . In the present paper we shall mostly work with the spatial components S ij of the spin tensor. This presents two advantages: first, it somewhat simplifies the algebra and the index structure by getting rid of the Levi-Civita tensors, and second, it makes the check of the Lorentz invariance more straightforward, since S ij is directly the spatial components of the tensor S µν . We provide in Appendix A the explicit link between S ij and the spin variable S i FBB .

III. POST-NEWTONIAN METRIC AND EQUATIONS OF MOTION
As in Paper I, we use harmonic (or De Donder) coordinates, defined by the gauge condition ∂ ν ( √ −gg µν ) = 0. The PN iteration of Einstein's field equations developed up to 3.5PN order in Ref. [53] is valid for a generic matter source, the only hypothesis being that the stress-energy tensor T µν must have compact support. We will need here the full 3.5PN metric, which is parametrized in terms of elementary potentials as: 2 From now on, spatial indices i, j, . . . will be raised and lowered using the Kronecker metric δ ij , and we will write them indifferently as upper or lower indices. The potentials are defined by means of the usual retarded (R) inverse flat d'Alembertian operator, where the sources f are built with the following matter quantities, as well as with products of derivatives of lower-order potentials. The potentials required for the 2.5PN equations of motion were already used in Paper I and read: while the potentials required for the 3.5PN order are given the following definitions: Notice the difference of structure between the sources of these potentials: some sources are proportional to one of the σ quantities and are compact-supported, while others are only proportional to metric potentials. We will call the latter non-compact supported, since their source extends in all space. Next, in keeping with notations of Paper I, we rewrite the covariant equation of motion (2.6b) in the following 3 + 1 form: The leading order of the left-hand side of the force law (3.5) is simply the ordinary acceleration a i = dv i /dt. The Newtonian-like linear momentum and forces are With transparent meaning we shall often call F i the "geodesic part" of the force law (or rather the acceleration), while F i will be referred to as the "Papapetrou part", deviating from geodesic motion. Expanding the metric in terms of the elementary potentials (3.4), one gets the following 3.5PN expressions for P i and F i , which would correspond to the coordinate acceleration of a particle following a geodesic motion: Specializing now to our problem, when computing for instance the equations of motion for the body 1, the velocity v i is to be replaced by v i 1 and the right-hand sides are to be evaluated at y 1 . However, beware that the metric potentials are generically singular at the location of the two particles 1 and 2 where they are meant to be evaluated. This evaluation is thus given a sense through the Hadamard "partie finie" regularization procedure explained in Ref. [54], which is of course to be performed after computing the derivatives of the potentials. We adopt the so-called pure Hadamard-Schwartz prescription [45] for the practical implementation of this regularization. In particular, we use the distributive rule for computing regularizations of products of potentials in Eqs. In this work, we restrict ourselves to contributions to the coordinate acceleration that are of linear order in spins. Considering Eq. (3.5), these linear-in-spin contributions to the geodesic part of the acceleration can have two origins at this level: first, from the spin contributions to the elementary metric potentials (3.4), which we will have to compute at the required order, and secondly, from the time derivative dP i /dt, since it is understood that order-by-order replacements of the accelerations are to be performed, which include in turn spin contributions starting at O(3). Table I indicates which spin parts of metric potentials are needed at which order, for the present computation of the next-to-next-to-leading order spin-orbit contributions to the equations of motion. As explained in Sec. IV, some potentials were computed in all space, but for the most non-linear ones, only their regularized values at the location of one of the particles could be computed. In the latter case, the calculation is different for each different derivative structure of the potential, and Table I gives the list of these needed derivatives.
On the other hand, the effective Papapetrou part of the acceleration or force F i , which corresponds to non-geodesic motion, can be found in Paper I up to 2.5PN order in terms of the spin variable S i FBB . We give here its complete expression up to 3.5PN order, in terms of the spin tensor S ij and the various metric potentials. Defining and expanding the Riemann tensor in Eq. (3.6c) in terms of metric components, themselves expanded in terms of the elementary potentials (3.4), we get: (3.10c) If we are computing the Papapetrou acceleration of body 1, F i 1 , we have to replace in the latter expressions S ij by S ij 1 and v i by v i 1 . Because of the explicit spin factor, at linear order in spins only the non-spin parts of the potentials are needed, and S ij 2 does not appear in this part of the acceleration. The factor m ≡ m 1 in Eq. (3.9) is at the origin of all terms depending on the spins by unit mass, i.e. through S ij 1 /m 1 or S ij 2 /m 2 in the final results. Note that terms with two velocities appear in Eqs. (3.10) because of the replacement of S 0i 1,2 by S ij 1,2 according to (2.9). Again all products of singular potentials have to be regularized following the rule of the pure Hadamard-Schwartz regularization, notably the distributivity rule (3.8). 3 However we shall find in Sec. V that the computation of the value at 1 of one particular potential, namely ∂ jkŶi , which is especially singular, a priori requires the use of dimensional regularization.
As we said, since the Papapetrou part of the acceleration already includes an explicit spin factor, only the non-spin parts of the potentials are required there, and most of these have already been computed in previous works such as [53,55]. Notable exceptions are the second spatial derivatives of the potentialsR i andŶ i , which are defined in Eqs. (3.4); we address the computation of these new regularized potentials in Secs. IV and V.

A. Notation and general points
Throughout the rest of this paper, we use the following notation. The trajectories of the two bodies, with masses m 1 and m 2 and spin tensors S ij 1 and S ij 2 , are denoted y 1 (t) and y 2 (t), and their coordinate velocities v 1 (t) and v 2 (t). The inter-body distance is r 12 = |y 1 − y 2 |, and we set n i 12 = (y i 1 − y i 2 )/|y 1 − y 2 |. For a generic field point x, we pose r 1 = |x − y 1 |, n i 1 = (x i − y i 1 )/|x − y 1 | and similarly for 2. We denote by ∂ i , ∂ 1 i and ∂ 2 i the partial derivative with respect to x i , y i 1 and y i 2 . The symbol 1 ↔ 2 means the expression obtained by the exchange of the two particles. Lengthy calculations are performed with the help of the scientific software Mathematica R , supplemented by the package xTensor [56] dedicated to tensor calculus.
We systematically expand the retardations inside the retarded integral (3.2), truncated to the appropriate order. For instance, truncating at O(3), we have: In this expression, time derivatives may be pulled out of the integrals, and it is understood that accelerations and time derivatives of the spin are to be replaced by the already computed lower order equations of motion and spin precession equations. Note that the second term in this expansion has no dependence to the field point x.
Next, notice that we are using Dirac delta functions in the stress-energy momentum tensor (2.7), while it makes no sense as a distributionà la Schwartz when acting on the class of general functions in the problem, which are generically singular at the locations of the two particles. We refer to [53] for a complete discussion of this issue, which was dealt with by defining a class of "pseudo-functions" and using Hadamard regularization to define the value of a singular function at one of its singular points. In our case, we will give a sense to Dirac delta functions inside integrals, following the rule where δ 1 = δ(x − y 1 ) and (F ) 1 means the Hadamard "partie finie" of F at the point y 1 [54]. This rule is extended in an obvious way to derivatives of delta functions, for instance Other issues are the problem of distributivity of the partie finie already mentionned, and the treatment of the distributional parts of derivatives. Here we shall not use the "extended Hadamard regularization" introduced in Ref. [54]; instead we will apply primarily the "pure Hadamard-Schwartz" prescription described in Ref. [45] which is sufficient for most of our computations. The pure Hadamard-Schwartz regularization constitutes the core of the most powerful and fundamental regularization procedure which is dimensional regularization and has been successfully applied to the problem of the equations of motion in Refs. [45,57]. Thus, in principle, the result of the pure Hadamard-Schwartz regularization of Poisson-type integrals is to be supplemented by a contribution from dimensional regularization which appears when the integral develops a pole ∝ 1/(d−3) in the dimension. The pole corresponds to the appearance of logarithmic divergences in the Hadamard regularization (see Sec. IV D). Here, we stick to the pure Hadamard-Schwartz prescription for the computation of the spin parts of potentials, where no problematic logarithms appear and where we expect that there is no difference between the extended Hadamard or pure Hadamard-Schwartz regularizations and the dimensional regularization (as it was the case up to 2.5PN order in the non-spin equations of motion [55]). On the other hand, as already mentioned we do use dimensional regularization for one particular potential in the Papapetrou part of the equations of motion, as discussed in Sec. V.
For the distributional part of derivatives of singular functions, we employ a particular form of the Gel'Fand-Shilov formula valid for homogeneous functions. Denoting by ∂ i the full derivative, including the distributional part, and by ∂ ord i the ordinary part, which also acts on the Dirac delta functions, we get, in the relevant cases of double derivatives and a simple function of type n L /r m (which is homogeneous of degree −m), where the distributional part is given when ℓ + m is an odd integer by and is zero when ℓ + m is even. 4 In the present context we apply the formulae (4.3)-(4.4) to the expansion of a singular function f around the two singularities y 1 or y 2 .

B. Leading-order potentials
We first deal with the lowest-order potentials, namely the computation of the leading spin-orbit part of the potentials V , V i andŴ ij . We give all the results for completeness, translating results from Paper I in terms of the spin tensor S ij . By expliciting the expression of the dipolar part (2.8) of the stress-energy tensor, one may check that the leading order for the spin part of the sources σ, σ i , and σ ij is: where (ij) indicates symmetrization, δ 1 stands for δ(x − y 1 ), and the gradients are taken with respect to the field point x. Given these expressions for the sources of potentials, and using ∆(1/r 1 ) = −4πδ 1 (with r 1 = |x − y 1 |), Eqs. (4.5) then yield: Notice that the leading-order of the spin part of these potentials has the typical dipolar structure ∂(1/r), while it is simply 1/r (monopolar) for the non-spin part. Note also the important fact that V S starts at O(3) and has no O(1) contribution. This makes the structure of the spin part of potentials different from the one of the non-spin part, and shifts a number of contributions to higher PN order. In particular, the non-compact support part . This in turn implies that the leading order ofŴ S ij has the simpler structure of a compact-supported term, which makes a number of sources of higher potentials less non-linear, such as the non-compact-

C. Higher order compact-support terms
We turn now to the higher-order PN computation of the compact-supported parts of potentials, i.e. whose sources include σ, σ i or σ ij . We will take the potentials V S as an example and give only a schematic view of these computations, since their implementation is relatively straightforward. By expliciting the expression (2.8) of the stress-energy tensor, we get the following structure for the compact-supported sources: indifferently for bodies 1 and 2, where the subscripts M and D refer to the monopolar and dipolar parts of the stress-energy tensor. The quantityμ 1M already intervenes in the study of the equations of motion without spin, and its 3.5PN expression in terms of the metric potentials is given in Eq. (4.3) of Ref. [58]. Next, we Taylor expand the retardations (4.1) for each term. The time derivatives are pulled out of the integrals, which are then evaluated using the rules (4.2) for the delta functions and their derivatives, keeping in mind that the factors 1/ √ −g in (4.7) depend on the field point x. As a check of the method, we also computed all the compact-support parts of the spin potentials using a fully non-distributive prescription, i.e. keeping the dependence in the field point x as long as possible (also when evaluating the quantities such as ν 1D ), and taking the Hadamard partie finie at the very end of the calculation. We obtained no difference with the pure Hadamard distributive prescription (3.8).
Notice that V S , for instance, gets contributions from both the dipolar and monopolar parts of the stress-energy tensor. For the dipolar part, the quantitiesμ 1D , ν 1D and ν i 1D all display an explicit spin factor. For the monopolar part, spin terms appear indirectly in two ways: through the spin contributions to the metric potentials in the expression ofμ 1M , and through the acceleration replacements when evaluating time derivatives.

D. Non-compact support terms
The non-compact supported potentials are more complicated, since their source is itself made of potentials and does not allow the simple integration procedure with delta functions. However, if the source has a simple structure, some integrations can be computed explicitly to obtain a solution valid for a generic field point x. One encounters often the source structure ∂(1/r)∂ 2 (1/r), which arises at leading-order in all ∂A N S ∂B S terms, where A and B are chosen among {V, V i } and {V, V i,Ŵij } respectively, as can be checked from the leading-order structure (4.6). Using the same techniques as in previous works (e.g. [53,55]), in particular using the function g such that ∆g = 1 r 1 r 2 , (4.8a) g ≡ ln(r 1 + r 2 + r 12 ) , (4.8b) we find the following relations: where for instance ∂ 1 i = ∂/∂y i 1 . With these relations in hand, completing the results for the compact-support parts, we are able to straightforwardly compute the leading order O(1) of X S ,R S i ,Ẑ S ij , as well as the 1PN relative order O(3) ofŴ S ij , as indicated in Table I. We now turn to the computation of the partie finie of Poisson integrals, of the type of the first and third term in (4.1). 5 When direct integration is not possible (at least using known results) for the non-compact support part of non-linear potentials, one by-passes the difficulty by directly computing the Hadamard regularized value of the integral at the points y 1 or y 2 . By doing so, one will not have access to the full information about the potentials (hence the metric) everywhere in space, but only the information relevant for the computation of the equations of motion. Since the computation is different for each derivative of the potential, one has to figure out which derivatives of which potential will be needed in the final computation. Table I already gave the required spin parts of potentials. Furthermore we find that two new non-spin (NS ) potentials, featuring two derivatives, have to be computed in addition to known previous results [53]: (4.10) Following Ref. [54] (to which we refer for details and definitions), we define generic Poisson integrals and twice-iterated Poisson integrals: where the source F is singular at y 1 and y 2 . The symbol Pf stands for the Hadamard partie finie for integrals diverging near the particle locations. This regularization procedure consists in removing from the integration domain two spherical balls of radius s surrounding the two singularities, substracting the purely divergent part, and taking the limit s → 0.
The results for the extended Hadamard partie finie (in the sense of Eqs. (5.3)-(5.4) in [54]) of the Poisson potentials P and Q are: In these formulae, the first term would correspond to the naive prescription of taking directly x = y 1 inside the integrals (4.11), and the second term features a regularization constant ln(r ′ 1 ) coming from the singular limit x ′ → y 1 . Notice that the constant s 1 automatically cancels between the two terms, while the Pf term may induce a dependence on s 2 , so that the result depends on two constants r ′ 1 and s 2 . For the present work we had to extend these formulae to the case of a double gradient, in order to compute the potentials (4.10). Using the same method as in Ref. [54], we obtained the corresponding required expressions for the double derivative: One readily checks that these expressions yield the correct results when contracted with δ ij : namely (∆P ) 1 = (F ) 1 and (∆Q) 1 = 2(P ) 1 . Crucial for this check, notice the last term in Eq. (4.13a) which stems from a distributional contribution obtained when evaluating the two derivatives under the integral (4.11a). It can be absorbed by replacing, in the first term of (4.13a), the ordinary factor (3n ′ ij 1 − δ ij )/r ′ 3 1 by the distributional derivative ∂ ′ ij (1/r ′ 1 ). With these tools in hand, we were able to compute all the needed non-compact support parts of potentials. Importantly, we get no contribution of the singular parts proportional to ln(r ′ 1 ) in Eqs. (4.12), in any of the spin parts of potentials listed in Table I. This gives a strong indication that the Hadamard regularization (actually the pure Hadamard-Schwartz version of it) is sufficient to deal with all spin parts of potentials. This is not surprising because the calculation we are doing is only of 2PN relative order, while past experience in the non-spin case [53,55] says that Hadamard's regularization fails no earlier than at 3PN relative order.

V. DIMENSIONAL REGULARIZATION FOR ONE POTENTIAL
By contrast, we have found that in the Papapetrou part of the acceleration, which is specific to the spin-case (and outside past experience), we do get a contribution proportional to ln(r ′ 1 ) in one of the new evaluations we had to perform, namely when computing the double-derivative potential (∂ j ∂ kŶ N S i ) 1 at Newtonian order using Eq. (4.13a). All the other NS potentials, including the other new one (∂ j ∂ kR N S i ) 1 we had to compute, could be safely obtained with the Hadamard regularization of Sec. IV.
The appearance of this ln(r ′ 1 ) tells us, on the contrary, that Hadamard's regularization is insufficient for the potential (∂ j ∂ kŶ N S i ) 1 , and that we a priori need dimensional regularization in order to get rid of possible ambiguities in the final equations of motion, as was the case of the non-spin 3PN equations of motion in [45,53]. On the other hand, when considering directly the Papapetrou contribution to the acceleration, we see that the dangerous contribution exists in only one term, with index structure S jk 1 (∂ i ∂ jŶ N S k ) 1 , and we find that the problematic logarithms ln(r ′ 1 ) cancel out in the final result because of the antisymmetry of S jk 1 . Therefore, we expect beforehand the pole ∝ (d − 3) −1 we shall obtain in dimensional regularization to actually vanish in the final acceleration. Nevertheless, even though the result will be pole-free, we know that the finite part ∝ (d − 3) 0 should play a crucial role, and therefore we must a priori apply the powerful but tedious procedure of dimensional regularization to the problematic potential (∂ j ∂ kŶ N S i ) 1 . We refer to [45] for precise definitions and technical details about the method, which consists in studying the problem in d spatial dimensions, treating d as a complex variable, and taking the analytical continuation of the obtained results when d → 3. We will associate a superscript (d) to quantites which are defined in this d dimensional setting. One considers a class of functions regular everywhere except at the points y 1 and y 2 , and admitting expansions of the type (∀N ∈ N): and similarly around y 2 , with ε ≡ d − 3, and q 0 , q 1 , p 0 ∈ Z. As long as there are no terms in these expansions with both p < 0 and q = 0, and that there is no angular dependence of the p = 0, q = 0 term in the expansion, which is always the case in practice, the presence of the factor r qε 1 allows one to write, by application of the analytical continuation on ε: which justifies the use of distributivity in our previous Hadamard three-dimensional computations. Indeed, it was checked in previous work [45] on the 3PN non-spin equations of motion that the pure Hadamard prescription agreed with the d → 3 limit of the dimensional prescription for all the terms except those developing a pole. As explained in Sec. IV C, although we did not perform a full dimensional regularization computation of all quantities to be evaluated at y 1 or y 2 , we checked that distributive and fully non-distributive prescriptions for the computation of compact-supported potentials yielded the same final result, which gives us confidence in this distributive prescription. 6 The treatment of Poisson integrals such as (4.11) goes as follows. Defining P (d) (x) as: wherek ≡ Γ d−2 2 /π d−2 2 , analytical continuation allows us to take derivatives under the integral, to set apart the distributional contribution, and to make directly the replacement x → y 1 under the remaining integral. We obtain: where the second term on the right-hand side corresponds to the distributional contribution.
0,0 , according to the expansion (5.1), and must not have any angular dependence (since dimensional regularization does not include angular averaging in its definition, unlike Hadamard's partie finie). We set here all the terms r qε 1 1 f (ε) 0,q (n 1 ) with q = 0 to 0 since they are cancelled by analytical continuation (and have in general an angular dependence). We also assume that 1 f (ε) 0,0 has no pole in ε which is true, since no poles are generated at the level of the sources of potentials for this calculation.
When d → 3 we have F (d) → F , where F is the corresponding function in the Hadamard 3-dimensional context. That function admits the expansion The Hadamard partie finie of the function is defined by the angular average over the unit direction n 1 of the term p = 0, say The 3-dimensional coefficients 1 f p (n 1 ) are related to the ε → 0 limits of the d-dimensional p,q (n 1 ) appearing in Eq. (5.1) by which in turn gives us the link between the analytical continuation when d → 3 of F (d) (y 1 ) and the Hadamard partie finie (F ) 1 : Next, we address the link between the Hadamard-regularized Poisson integral (4.11) and its dimensional regularization version (5.3). Defining: and following the same steps as in Ref. [45], we obtain, working at the zeroth order in ε, The result of the PN iteration of the metric, starting with the d + 1-dimensional Einstein equations, and the corresponding definitions for the metric potentials are given in Sec. II of Ref. [45]. In particular, the d-dimensional definition of the dangerous potentialŶ i is given by (2.12f) there. 7 We applied the previous method to the computation of the quantity D(∂ jkŶi )(1) defined above. This required computing the singular expansion of the needed non-compact support sources in d dimensions up to the order p = 0, which was done following the methods of Ref. [45]. Our result for the pole part is quite compact, Since this expression is symmetric by exchange of the two indices i and j (or k), it doesn't contribute to the final equations of motion where only the contraction S jk 1 (∂ ijŶ N S k ) 1 is involved. This confirms that the final equations of motion in dimensional regularization are directly pole-free. Furthermore, we also checked that the finite part ∝ ε 0 coming from Eq. (5.10) also cancels out in the final result (i.e. after contraction with the spin tensor S jk 1 ). Although the latter fact could not a priori be guessed beforehand, it implies that finally the dimensional regularization was superfluous for the computation of the dangerous term S jk 1 (∂ ijŶ N S k ) 1 . This result gives us further confidence that the pure Hadamard-Schwartz regularization is sufficient for all our calculations.

VI. RESULTS FOR THE EQUATIONS OF MOTION AND PRECESSION
Since higher-order results take the form of quite long formulae with many similar terms, we tried to adopt a systematic presentation. First, we split the expressions in terms of powers of G and according to the powers of the masses m 1 and m 2 . Each spin-dependent term is then denoted using the convention that α i p,q gathers all the terms featuring G p+q m p 1 m q 2 . For any two vectors a and b, we use the notation (ab) for the scalar product, i.e. (ab) ≡ a · b = a i b i , and we define (Sab) ≡ S ij a i b j .
Before presenting the results, let us address an important point: every time an expression has at least one free index, it might admit several equivalent writings, in terms of the vectors appearing in the problem and of the spin tensors. Indeed, any of the final results, which are 7 Notice that because we have already seen that the pole 1/ε will necessarily cancel out in the final result, we do not have to worry about extending the Papapetrou part of the force (3.10) to d dimensions. I.e., we can evaluate the term ∂ jkŶi assuming it carries the 3-dimensional coefficient, which is +8 in this case, see one of the terms in Eq. (3.10c).
functions of time only (such as the acceleration of one of the bodies or the time derivative of one of the spins), can be written with the three vectors n i 12 , v i 1 , v i 2 , and with at most one occurence of either one of the spin tensors S ij 1 and S ij 2 since we are working at linear order in spins. As the spatial indices run on three different values, we have the two identities: where the brackets indicate antisymmetrization over the indices, where S ij is one of the spin tensors and where (a, b, c) is a permutation of the set of three vectors (n 12 , v 1 , v 2 ). The number of identities of this kind that one must take into account depends on the number of vectors and tensors at disposal, and on the number of free indices. As an example, (6.1a) with S = S 1 , a = n 12 , b = v 1 and c = v 2 gives, once expanded: Hence, one must keep in mind that there is no unique writing for the results we are going to present, and take this into account when comparing two expressions. This becomes particularly important when using the method of undetermined coefficients (for instance when looking for a contact transformation between harmonic and ADM variables, see Sec. VII D): the system of independent equations that these coefficients have to solve is to be determined only after taking into account the complete list of these identities. Notice also that the use of an antisymmetric spin tensor reduces the number of these identities compared to the use of a vector spin variable and a Levi-Civita symbol, which was one of our motivations for changing the spin variable with respect to Paper I. A straightforward but cumbersome work around of this problem when comparing results is to project everything in an arbitrary orthonormal basis.

A. Spin evolution equation
First, we give the spin evolution equation, i.e. the equation giving the time derivative of the spatial components of the spin tensor S ij , up to 2PN order. This result was in fact already contained in Paper I, under the form of a precession equation for the spin vector S FBB , and it is a mere matter of traduction between the spin variables. Beware that our spin tensor is not of conserved norm, S ij S ij = const beyond leading order. This equation is needed each time we perform order-by-order reduction when evaluating time derivatives.
Since we are working at 2PN relative order, it is sufficient to know the 2PN or O(4) spin evolution equation. This also means that the amount of non-linearity in this computation is less that in the computation of the acceleration. Indeed, expliciting Eq. (2.6a) in terms of potentials, we see that we need only 2PN potentials, and we can ignore the spin contributions in these potentials, working at linear order in spins. However, notice that, since the leadingorder spin contributions to the total angular momentum of the system is of the form S 1 /c + S 2 /c, with S 1,2 spin vectors, the order of this spin evolution equation required for finding a conserved total angular momentum at 3.5PN order is not 2PN but 3PN.
Defining dS ij we obtain for the 1PN order: and for the 2PN order, splitting the result in terms of powers of G and occurence of the masses as explained above: where Note that, with the series of potentials already computed for the 3PN equations of motion without spins [45,59], we are able to control the precession equations up to the next 3PN order; this will be investigated in future work.

B. Acceleration
The spin contributions in the acceleration have the following structure, with our PN counting valid for maximally spinning-objects: where the S subscript indicates the spin-orbit contributions, and the SS subscript indicates contributions that are quadratic in the spins and which we neglect in this work. 8 For the leading-order spin contributions, we get m 2 3S ij 1 n j 12 (n 12 v 12 ) + 6(S 1 n 12 v 12 )n i 12 − 3S ij 1 v j 12 8 Notice that, when doing the calculation in the original harmonic coordinates, one actually obtains some 3PN spin-orbit contributions. However, these are pure gauge, as was already explained in Appendix A of Ref. [60], and are eliminated by the gauge transformation x µ → x µ + δX µ with δX 0 = 0 and whose effect on the accelerations is δa i 1 = δa i 2 = d 2 δX i /dt 2 + O(8). This gauge transformation obviously respects the harmonicity condition in a perturbative sense, since δX µ = O(8). +m 1 6S ij 2 n j 12 (n 12 v 12 ) + 6(S 2 n 12 v 12 )n i 12 − 4S ij 2 v j 12 .

A. Conserved Energy
An important feature of all the spin-orbit contributions we have computed, is that they are associated with the conservative part of the dynamics, i.e. obtained when neglecting the radiation-reaction dissipative terms associated with gravitational radiation. Therefore these contributions should allow for the existence of a set of conserved quantities, namely energy, angular momentum, linear momentum, and center-of-mass integral. Here we have checked that the new terms computed in the equations of motion in Sec. VI B admit corresponding contributions in the conserved energy, which take the following structure: The leading order spin-orbit contribution to the energy reads while the 2.5PN result (1PN relative) was given in Paper I, and reads, once translated into the spin tensor variable S ij : with γ = (1 − V 2 /c 2 ) −1/2 the Lorentz factor. We denote the trajectories of the two bodies in the frame (F ) by y µ 1 = (ct, y 1 ) and y µ 2 = (ct, y 2 ), and by y ′µ 1 = (ct ′ , y ′ 1 ) and y ′µ 2 = (ct ′ , y ′ 2 ) in the boosted frame (F ′ ). The point is that we cannot compare them directly, because the simultaneity surfaces are different in the two different frames. To give a sense to simultaneity between the two frames, it is convenient to define an auxiliary event Ω(ct, x), which can be choosen at will -it could be for instance the point of coordinates (ct, x = 0) in (F ), and whose coordinates in (F ′ ) will be (ct ′ , x ′ ). We define the two events P 1 (ct, y 1 (t)) and P 2 (ct, y 2 (t)) on the two worldines, simultaneous to Ω in the frame (F ), and similarly two events Q 1 (ct ′ , y ′ 1 (t ′ )) and Q 2 (ct ′ , y ′ 2 (t ′ )), simultaneous to Ω in the frame (F ′ ). Ω plays the role of an observer, for which the equations of motion evaluated on the simultaneity surfaces (Ω, P 1 , P 2 ) and (Ω, Q 1 , Q 2 ) must take the same fonctional form, as explained below. Next, we define the times τ 1 and τ 2 , such that the coordinates of Q 1 and Q 2 in (F ) are (cτ 1 , y 1 (τ 1 )) and (cτ 2 , y 2 (τ 2 )). We have, by construction, y ′ 1 µ (t ′ ) = Λ µ ν y ν 1 (τ 1 ) and similarly for 2. The link between y ′ 1 (t ′ ) and y 1 (t) is provided in Ref. [61]. One obtains successively where the right-hand sides are evaluated at t, and (V r 1 ) means V · (x − y 1 ). The expressions for the velocity and acceleration are obtained by taking the time derivative according to ∂ ′ t = γ∂ t + γV i ∂ i . Similar expressions hold for 2, and transformations for quantities such as r 12 or n i 12 are deduced in a perturbative sense from the first of these formulae. We turn now to the transformation rules for the spin tensor. As explained in Paper I, we have for any function f (t): which can be applied to the components S 0i 1 and S ij 1 of the spin tensor. Now, since these are the components of a contravariant tensor, we have Combining these expressions, we arrive at the desired transformation rule for S ij 1 : where S 0i 1 can be further eliminated using the spin supplementary condition (2.9). This transformation is somewhat simpler than (A13)-(A15) in Appendix A of Paper I because the transformation for the spin vector S i FBB is not as simple as in Eqs. (7.10). Lorentz invariance means that the acceleration in the boosted frame, obtained by transforming all variables according to Eqs. (7.8) and (7.11), when systematically truncated at the requested PN order, must take the same functional form as in the original non-boosted frame. We have verified that our result for the acceleration (6.11) passes this test. Note however that the test leaves a lot of freedom in the acceleration. In particular, one might add at the highest level 1/c 7 in the acceleration any quantity depending only on the relative velocity v i 12 = v i 1 − v i 2 and still pass this test. Finally let us remark that the pure Hadamard Schwartz regularization which has been used for all the terms in this computation, finally yields equations of motion which are manifestly Lorentz invariant. This is another indication that there is no need for using a more sophisticated regularization such as dimensional regularization in this problem.

C. Test-mass limit
Another important check of our result is to take the test mass limit and show that we recover the equations of motion of a test particle orbiting in a black hole background. More precisely, we can check that: 1. In the limit where one of the two bodies (say body 1) is a test particle with m 1 → 0 and is spinless (S ij 1 = 0), its acceleration (6.7) reduces to that of a spinless test particle orbiting around a Kerr black hole; 2. In the limit where body 1 is a test particle with spin (i.e. m 1 → 0 with constant ratio S ij 1 /m 1 ), Eq. (6.7) reduces to the equations of motion of a massless spinning particle orbiting around a Schwarzschild black hole.
Note that in case 2. we need to work with a Schwarzschild black hole because we restricted ourselves to spin-orbit effects. In this section, it will be more convenient to rewrite the equation of motion (3.5)-(3.6) so as to make apparent the coordinate acceleration: where we recall that F µ = − 1 2mc R µ νρσ v ν S ρσ and u 0 = 1/ −g ρσ v ρ v σ /c 2 .

Spinless test particle: equivalence with Kerr geodesics
In the limit where m 1 → 0 and S ij 1 → 0, then v i 2 = 0 and S ij 2 = const is a trivial solution of the equations of motion and the (spin part of the) acceleration of body 1 given by Eq. (6.11), simply reduces to (a i 1 ) S = G c 3 r 3 12 6(n 12 v 1 )A i 2 − 4B i 2 + 6(S 2 n 12 v 1 )n i 12 − 6 G c 5 r 3 12 (n 12 v 1 )(S 2 n 12 v 1 )v i 1 which compared to the harmonic-coordinate results in Paper I at the next-to-leading order O (5). In this formalism, the canonical structure is defined for the variables x a (positions of the two bodies), p a (canonical momenta) and S a (canonical spins). Notice that these canonical spins are of conserved norm. The Poisson brackets for these variables read: = c δ ab ε ijk S k , (7.21b) and all the other Poisson brackets vanish. Beware of the additional c factor appearing in the second bracket: our PN counting is different from that in [38], and we have set S i = cŜ i HS to keep close to the convention we use, see Eq. (1.1). One may as well employ a conserved-norm antisymmetric spin tensor defined (exactly) by: The dynamics is entirely contained in the Hamiltonian of the two-body system, expressed in terms of these canonical variables. Besides the recent result for the spin-orbit 3.5PN spinorbit contributions [38], the non-spin part of the Hamiltonian up to 2PN order can be found for instance in [63], and the spin-orbit part up to 2.5PN order is given in [32]. As explained before, since we are working at 2PN relative order with respect to the leading order spinorbit contribution, we will need to know the non-spin part of the dynamics only up to 2PN order. The time derivative of any quantity f is obtained as: df dt = {f, H} + ∂f ∂t , (7.23) where the second term accounts for a possible explicit time dependence of f , absent in our problem. Beware again that, because of the additional c factor in the Poisson brackets for the spin, the PN truncation of the Hamiltonian must be handled carefully. It is well known that, when comparing the results of the two formalisms, it is necessary to perform a contact transformation of the worldlines of the particles, i.e. a time-dependent shift of the worldlines that is not a global coordinate transformation. Thus, one must keep in mind that the positions y a and x a differ at higher-order. We write this symbolically as y = Y(x, p, S) . On the other hand, we may use (7.24) and (7.25) and the link between the spin variables (7.30) to translate the harmonic-coordinates result (6.11) for the acceleration in terms of the ADM variables as: a = a(y, v, S) = a(Y(x, p, S), V(x, p, S), S(x, p, S)) , (7.27) and this expression must be identical to (7.26). Thus, one is led to look for an extension of the contact transformation (7.24) that realizes this identity. Focusing on the non-spin and spin-orbit parts, the structure of this contact transformation is the following: The lower-order already known expressions for this contact transformation can be found in [63] and [32]. Notice that, as when using the Hamiltonian to compute time derivatives in the ADM setting, we need the non-spin part of this contact transformation up to 2PN only. To present these formulae, we use the canonical spin tensor variable (7.22) to get the same index structure as in the other results of the present paper, and we adopt the convenient notation π a ≡ p a /m a to shorten the expressions. We find: Y 2PN 1 = Gm 2 1 2 π i 1 (n 12 π 2 ) − 7 4 π i 2 (n 12 π 2 ) + 1 8 n i 12 5π 2 2 − (n 12 π 2 ) 2 + G 2 m 2 4r 12 n i 12 (7m 1 + m 2 ) , (7.29a) ij 2 π j 2 + (n 12 π 2 )S ij 2 n j 12 + 1 2 (S 2 n 12 π 2 )n i 12 . (7.29c) We need also the conversion rule between the harmonic spin tensor and the ADM canonical spin, say S = S(x, p, S). The conversion rule from the conserved norm variable S c FBB to S is given in Eqs. (6.4) and (6.5) of [32], and the conversion rule between this conserved norm spin and S FBB is given in Eq. (7.4) of [31]. Finally, knowing the link between S FBB and our spin tensor, Eq. (A1) in Appendix A, we get: