We revisit the dynamics of the postNewtonian PN twobody problem for two inspiraling compact bodies. Starting from a matteronly reduced Hamiltonian we present an adapted framework based on the Lie series approach enabling the derivation of complete perturbative solutions within the conservative sector. Our framework supports both circular and eccentric orbits and is applicable to any perturbation respecting rotational invariance and timeindependence. In the context of the ArnowittDeserMisner ADM canonical formalism this includes up to at least 3PN order and local terms beyond. We provide an example application at 2PN recovering classical periapsis advance and orbital period corrections alongside the full orbital evolution in time coordinates. We discuss eventual extension to spinning and timedependent systems.
['General Relativity and Quantum Cosmology']
Hamiltonian normal forms for the postNewtonian binary problem
{'Hamiltonian normal forms for the post-Newtonian binary problem': "C. Aykroyd, 1, ∗ A. Bourgoin, 1 and C. Le Poncin-Lafitte 1 1 SYRTE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universit´es, UPMC Univ. Paris 06, LNE, 61 avenue de l'Observatoire, 75014 Paris, France (Dated: September 20, 2024)  \nWe revisit the dynamics of the post-Newtonian (PN) two-body problem for two inspiraling compact bodies. Starting from a matter-only reduced Hamiltonian, we present an adapted framework based on the Lie series approach, enabling the derivation of complete perturbative solutions within the conservative sector. Our framework supports both circular and eccentric orbits, and is applicable to any perturbation respecting rotational invariance and time-independence. In the context of the Arnowitt-Deser-Misner (ADM) canonical formalism this includes up to at least 3PN order and local terms beyond. We provide an example application at 2PN, recovering classical periapsis advance and orbital period corrections, alongside the full orbital evolution in time coordinates. We discuss eventual extensions to spinning and time-dependent systems.", 'I. INTRODUCTION': "The Einstein field equations form the cornerstone of the general theory of relativity (GR), offering an elegantly profound description of interacting matter, energy, and geometry. Yet these equations are notoriously complex due to their nonlinear, tensorial and interdependent nature. As a consequence, the discovery of exact solutions is chiefly limited to highly idealized symmetrical scenarios. Real astrophysical phenomena, however, are rarely encapsulated by such pristine conditions. Gravitationally radiating systems are a prime example of inherently asymmetric systems, particularly evident in the dynamics of compact binaries. Their study typically segments the process into distinct regimes, each requiring its own set of numerical or analytical approximation tools. The PN framework, popularized from the earliest days of GR [1-4], excels in addressing the inspiralling two-body problem. Initially, PN approximations assume wide, bound orbits with weakly or moderately relativistic speeds. Yet in practice, these solutions have been 'unreasonably effective' [5] in describing a range of compact systems, achieving remarkable accuracy up to the innermost stable circular orbit. This success is partly attributed to GR's 'effacing' property, which permits the omission of self-gravitation and internal structure at lower orders.  \nThe evolution of the PN framework has closely followed technological advancements, notably with the discovery of the Hulse-Taylor pulsar [6], which underscored the necessity for a formalism compatible with high-order calculations. The era of gravitational wave (GW) astronomy, inaugurated by the detection of GW150914 [7], has reinforced this need, aligning analytical methods with the capabilities of second-generation detectors such as Advanced LIGO [8], VIRGO [9], KAGRA [10] (LVK), and the International Pulsar Timing Array [11]. As we approach the advent of third-generation detectors, includ-  \ning LISA [12], the Einstein Telescope [13], and the Cosmic Explorer [14], further refinement of the PN framework remains imperative.  \nPN schemes typically rely on decoupling the field and matter degrees-of-freedom. This can be performed by 'integrating out' the metric either at the level of the field equations or directly at the action-via a combination of perturbative techniques such as multipolar expansions [15], asymptotic matching [16, 17], and Effective Field Theory [18-20]. Eventually, the aim is to establish a matter-only Fokker Lagrangian [21], or to directly compute the particles' equations of motion [16]. Canonical formulations like ADM instead search for a reduced matter-only Hamiltonian or Routhian functional [22, 23]. Accordingly, different harmonic and ADM-TT gauge formulations have been found to be in agreement [16, 24]. But while there has been a notable focus in this fieldmatter decoupling step of PN schemes, there remains a substantial gap in furnishing comprehensive analytical solutions of the reduced functionals. The quasi-Keplerian parametrisation [25-27] has been sucessfully determined up to 4PN to describe the orbits of non-spinning binaries. However, a systematic extension of this approach to broader perturbative situations including spinning or dissipative systems remains a challenge. In many other instances, orbital evolutions are only determined approximately, restricted to secular contributions or assumptions of (quasi-)circular orbits. A common approach employs the Lagrange Planetary Equations [28, 29], which, fundamentally, entails solving by hand the equations of motion for the Keplerian elements, and can hardly be achieved exactly. In the context of canonical formalisms, Hamilton-Jacobi methods [23, 30, 31]-typically represented in action-angle coordinates-provide a straighter path to determining long-term dynamical invariants, but often lack explicit time-coordinate solutions.  \nDeparting from the reduced matter-only canonical formalism, the Lie series approach can instead allow us to achieve an analytic portrayal of the dynamics that captures both long- and the short-term phenomena. Notably, the Lie framework provides a very structural path  \nfor extending 'secular averaging' procedures to high orders and to a broad range of perturbative scenarios. In [32], Lie methods have been employed for determining symplectic ADM centre-of-mass coordinate transformations. They have also been applied to solve the binary problem at 1PN order, in the spinless case [33] and in the spinning case [34, 35]. To our knowledge, a general two-body Lie framework that is directly applicable to post-Newtonian systems (at higher than first order) is not yet established.  \nIn this work, we aim to describe a canonical framework for orbital determination that can be readily applied to a wide range of perturbed two-body systems. We accordingly furnish parametric solutions for time-independent and rotational-invariant perturbations. These assumptions are valid in the ADM (spinless) conservative sector, up to at least 3PN order and to local terms beyond. Our provided solutions shall include the complete analytical orbital behaviour at some PN truncation order, including secular and oscillatory terms, and valid for fully eccentric systems. We shall prescribe a treatment which does not require explicit usage of action-angle coordinates, but which is instead established in terms of Keplerian functions or other parametrisations that can be easily adapted to the problem at hand. This freedom facilitates eventual applications of the framework to spinning systems and perturbations of non-gravitational nature. In Sec. III, we shall apply our framework to the 2PN ADMTT gauge Hamiltonian, not only re-deriving standard periapsis advance and orbital period corrections, but also thoroughly exploring the full long- and short-term orbital evolution. We then verify our results in a numerical example.  \nNotation and conventions.  \nWe presently introduce the notation used throughout the paper. Bold symbols are reserved for 3-vectors, decomposed into norm u = | u | (light-typeface) and direction ˆ u = u /u (overhat). Complex conjugation shall be denoted by the dagger operator ( · ) † . We shall adopt the following Poisson bracket convention:  \n{ f, g } := ∂f ∂ r · ∂g ∂ p -∂f ∂ p · ∂g ∂ r , (1)  \nfor phase space coordinates ( r , p ).", 'II. TWO-BODY PERTURBATIVE FRAMEWORK': "We dedicate this section to recalling the general Lie perturbative method and formulating it in the context of the post-Newtonian two-body problem. Consider a matter-only, local, autonomous Hamiltonian for two point-particles of the form:  \nH ( r , p ) = K -1 ∑ ℓ =0 c -2 ℓ H ℓ ( r , p ) + O ( c -2 K ), (2)  \nexpressed in a center-of-mass coordinate system C . The zeroth-order Newtonian term H 0 is assumed integrable, which is not the case for the higher-order ℓ -th PN perturbations c -2 ℓ H ℓ ( r , p ). Fortunately, not all is lost, and a complete time-integrated solution can often be obtained via the appropriate choice of coordinates.  \nThe essence of the Lie series method is to establish a coordinate system C ∗ in which the Hamiltonian becomes integrable up to an optimal order Q ≤ K -1:  \nH ∗ ( r ∗ , p ∗ ) = Q ∑ ℓ =0 c -2 ℓ H ∗ ℓ ( r ∗ , p ∗ ) + K -1 ∑ ℓ =0 c -2 ℓ R ℓ ( r ∗ , p ∗ ) + O ( c -2 K ), (3)  \nwhere the starred expressions refer to the new set of canonical variables in C ∗ . The new integrable Hamiltonians H ∗ ℓ are in mutual involution. The terms R ℓ are nonintegrable remainders which include eventual resonances from the low orders ℓ < Q ; in general, such resonances must be carefully handled (see e.g. [36] for details). The above expression is known as a Birkhoff normal form.  \nFor the class of systems we shall be considering in this paper, each Hamiltonian term shall be a rotationalinvariant scalar, depending only on r = ∥ r ∥ , p r := p · r /r and p = ∥ p ∥ . We shall demonstrate a posteriori that, remarkably, resonances are absent in our systems-namely, R ℓ = 0 , ∀ ℓ ≤ Q . Furthermore, we shall henceforth restrict H to be of PN order sufficiently low that Q = K -1. In this case, Eq. (3) is reduced to the integrable terms H ∗ ℓ only.  \nThe general strategy to obtain (3) is to consider a family of 'near-identity' coordinate transformations which are canonical by construction. Namely, we consider the group of Hamiltonian symplectomorphisms, coordinate transformations T g generated by the action of smooth Hamiltonian functions g . Equipped with the Poisson bracket these functions form a Lie algebra with elements L g = {· , g } ; each vector element may act infinitesimally on phase-space functions f :  \nδf = L g ( f ) = { f, g } . (4)  \nThe near-identity condition requires the generator to be of PN order O ( c -2 ), with potentially higher-order terms:  \ng ( r , p ) = K -1 ∑ ℓ =1 c -2 ℓ g ℓ ( r , p ) + O ( c -2 K ), (5)  \nwhere each g ℓ is explicitly independent of c . To recover the global coordinate transformation one then exponentiates L g , giving rise to the Lie series :  \nT g ( f ) = e L g ( f ) := ∞ ∑ ℓ =0 1 ℓ ! L ℓ g ( f ) = f + { f, g } + 1 2 { { f, g } , g } + . . . , (6)  \nwhere L ℓ g denotes the ℓ -th repeated composition of L g . By construction, these transformations preserve Poisson brackets. Note that the converse is not true, not all symplectic transformations can be expressed as the application of a Hamiltonian flow. Under this formalism, the phase-space coordinates transform as:  \nr = T g ( r ∗ ), p = T g ( p ∗ ), (7)  \nwith the inverse given by negating the sign of g .", 'A. Homologic equations': 'Time-independent canonical transformations preserve the value of the Hamiltonian, in which case it holds that H ∗ ( r ∗ , p ∗ ) = H ( r , p ) = T g ( H )( r ∗ , p ∗ ). Expanding the Lie operator on H and equating similar-order terms yields the homologic functional relations:  \nH 0 = H ∗ 0 (8a)  \n{ g 1 , H 0 } = H 1 -H ∗ 1 , (8b)  \n{ g 2 , H 0 } = H 2 + 1 2 { H 1 + H ∗ 1 , g 1 } -H ∗ 2 , (8c)  \n.  \n.  \n.  \nwhere the unknowns are the generator sequence ( g k ) and the Hamiltonian sequence ( H ∗ k ). The homologic equations can be solved sequentially from the top, plugging the solutions from the previous iterations; generically, they have the form:  \n{ g ℓ , H 0 } = P ℓ -H ∗ ℓ , (9)  \nwhere the perturbation P ℓ is some function of the solutions up to order ℓ -1. Taking the average along the Newtonian flow Φ H 0 t directly yields a solution for H ∗ ℓ :  \nH ∗ ℓ = ⟨P ℓ ⟩ , (10)  \nwhere  \n⟨ f ⟩ = lim T →∞ 1 T ∫ T/ 2 -T/ 2 f · Φ H 0 t d t . (11)  \nNext, we determine the generator by leveraging the integrability of the Keplerian system. Namely, H 0 is a maximal superintegrable system of 6-dimensional phasespace with 5 independent integrals of motion. These integrals can be specified in many forms; we shall supply some of them in the following section. In this manner, the system can be expressed in terms of five freely-chosen (independent) constants I = ( I 1 , . . . I 5 ), and a remaining Keplerian periodic degree-of-freedom, which we shall call θ . We stress that although these six parameters must be one-to-one with ( r , p ), they are not required to form a  \ncanonical set, allowing more flexibility in the specification of the results 1 . We may subsequently develop the left-side of the homologic equation via chain rule:  \n{ g ℓ , H 0 } = ∂g ℓ ∂θ { θ, H 0 } + ∂g ℓ ∂ I · { I , H 0 } , (12)  \nwhere the five brackets { I , H 0 } identically vanish. A formal solution is thus obtained from Eq. (9) by integrating over θ :  \ng ℓ = ∫ ( P ℓ -H ∗ ℓ ) ( d θ d t ) -1 kep d θ , (13)  \nwhere we have defined the Keplerian time-derivative  \n( d d t ) kep ≡ {· , H 0 } . (14)  \nThe primitive is taken to be average-free, which avoids the emergence of Poisson terms.  \nThe construction of an explicit parametric expansion to these integral solutions requires a convenient description of the system coordinates. For this we propose the Keplerian parametrisation of the following section.', 'B. The Keplerian parametrisation': "We shall now recall the Keplerian osculating parameters, which are typically employed in celestial mechanics (see e.g. [37]) and are convenient for the Lie formalism. We stress that this differs from the post-Keplerian formalism [25, 26], where Newtonian relationships are adjusted with PN corrections. Instead, one adopts phasespace relations which are exact, independently of the system perturbation. Namely, every given point ( r , p ) in phase-space C (or likewise in C ∗ ) is uniquely characterised by six orbital elements; for instance: the semimajor axis a , the eccentricity e , the true anomaly v , the argument of the periapsis ω , the inclination ι , and the longitude of the node Ω. These do not form a canonical set of variables, yet their Poisson brackets do not require PN corrections, facilitating the integration of the Newtonian-level derivatives in the homologic equation. We shall also introduce the eccentric and mean anomalies E and M , the longitude of the periapsis ϖ = ω +Ω, and the total orbital action L = √ a , all of which can be related to the aforementioned six elements. To specify all these parameters, we henceforth assume that the coordinates are appropriately rescaled, so that the Keplerian potential adopts the form:  \nH 0 ( p , r ) ≡ p 2 2 -1 r = -1 2 L 2 . (15)  \nAt each point in C , we introduce an ellipse tangent to the motion, lying within a coordinate plane perpendicular to the orbital angular momentum J ( r , p ) = r × p . Accordingly, we decompose the motion onto a right-handed orbital basis ( ˆ A , ˆ B , ˆ J ):  \nr = r (cos v ˆ A +sin v ˆ B ), (16a)  \np = 1 J ( -sin v ˆ A +( e +cos v ) ˆ B ) , (16b)  \nwith r = J 2 / (1 + e cos v ) and e = √ 1 -J 2 /L 2 . The vectors A = e ˆ A and B = e ˆ B are the periapsis and binormal vectors. Together with J = J ˆ J they form nine (dependent) first integrals of the Newtonian motion, of which A and B are only conserved for perfect inverse-square central forces. The last two can be readily expressed in functional form:  \nA = p × J -r r , B = J p -J × r Jr . (17)  \nThe orbital basis can be further related to a non-rotating reference orthonormal basis ( ˆ X , ˆ Y , ˆ Z ) via a rotation operator R :  \nˆ A = R ( ˆ X ), ˆ B = R ( ˆ Y ), ˆ J = R ( ˆ Z ). (18)  \nIn celestial mechanics, R is typically specified via z -x -z Euler angles (Ω , ι, ω ). Finally, the three anomalies v , E and M are defined to satisfy the following relationships:  \nM = E -e sin E (19)  \nand  \ncos v = 1 e ( J 2 r -1 ) , sin v = Jp r e , (20a)  \ncos E = 1 e ( 1 -r L 2 ) , sin E = rp r eL . (20b)  \nThese angles are chosen on the same Riemann sheet, ensuring that the difference between any two of them does not exceed π .  \nWe emphasize that the Keplerian parameters are to be seen as functions of phase-space and applicable to both the C and C ∗ coordinate systems. Accordingly, for each parameter q , we introduce their time-evolutions under Hamiltonian H and H ∗ :  \n˜ q ( t ) = q ( r ( t ) , p ( t )), (21a)  \n¯ q ( t ) = q ( r ∗ ( t ) , p ∗ ( t )). (21b)  \nWe shall refer to ¯ q as the secular evolution of q and to ˜ q as the complete time evolution of q . We can immediately deduce the following relationship:  \n˜ q = T g ( q ). (22)  \nThe evolution equation of each parameter can be extracted from the generalized Hamilton's equations,  \n˙ ˜ q = ˜ { q, H} , (23a)  \n˙ ¯ q = { q, H ∗ } , (23b)  \nvalid for non-canonical variables. In these equations, the application of tilde (resp. bar ) denote that expressions are evaluated in C (resp. C ∗ ) coordinates [cf. Eq. (21)], after computation of the brackets. The advantage of the Lie approach lies in that solving Eq. (23b) and then applying (22) is in general much easier than directly solving (23a). The Poisson bracket structure of each parameter q can be deduced from the provided definitions, and is already known from celestial mechanics 2 . We remark that these brackets remain unchanged under the symplectic transformation T g .", 'C. Solutions to the homologic equation': 'The integral solutions to the homologic equation can be explicitly solved for a large class of rotationally-invariant problems. Namely, we consider the asymptoticallyvanishing polynomial expansion  \nP ℓ = A ℓ, 0 + A ℓ, 1 r -1 + . . . + A ℓ,s r -s , (24)  \nexpressed such that the coefficients A ℓ, 0 , . . . , A ℓ,s depend solely on the Keplerian integrals of motion I (via use of the relationships from Sec. II B). In this case, we develop Eqs. (10, 13) with:  \ng ℓ = s ∑ k =1 A ℓ,k ( ξ k -M n 0 h k ) , H ∗ ℓ = s ∑ k =0 A ℓ,k h k , (25)  \nwith n 0 = ( -2 H 0 ) 3 / 2 .The basis elements turn out to be:  \nξ k = ∫ 1 r k ( d θ d t ) -1 kep d θ , h k = 〈 1 r k 〉 , (26)  \nwhere the primitive is chosen ensuring that ⟨ ξ k ⟩ = ( ⟨ M ⟩ /n 0 ) h k = 0. Explicit integrated expressions for ξ k and h k are provided in Appendix A.  \nWith this the Hamiltonian H ∗ = ∑ ℓ c -2 ℓ H ∗ ℓ and the generator g = ∑ ℓ c -2 ℓ g ℓ are now fully determined. The explicit coordinate-time evolutions r ∗ ( t ) and p ∗ ( t ) can hence be obtained from the integrable equations of motion; for instance, via the Kepler elements ¯ q ( t ) outlined in the previous section. Initial conditions in C ∗ -space are determined through application of the inverse Lie operator T -g . Eventually, the phase-space evolution in the original coordinate system C can be recovered via application of T g .  \nAs an immediate example, we shall apply the framework to the conservative sector of the ADM matter Hamiltonian at 2PN order.', 'III. APPLICATION TO THE CANONICAL ADM HAMILTONIAN AT 2PN': "Asymptotically flat coordinates in GR admit the Poincar'e group as a global symmetry [24, 31]. The generators of this group are realised as 10 surface integrals, conserved on-shell as dictated by Noether's theorem: the total energy E and total linear momentum Π , generators of space-time translations; the angular momentum ϑ , generator of rotations; and the generator of boosts Κ = Γ -t Π associated to the centre-of-mass constant Γ . In canonical formalisms, each of these quantities is seen as a function of phase-space. The ADM-TT gauge Hamiltonian is not manifestly Poincar'e invariant since most quantities are realised in nonlinear fashion. Instead, this invariance must be explicitly shown as done in [23, 25, 38-40]. However, for non-spinning particles, the ADM-TT gauge does manifestly respect the Euclidean group, and thus the momenta are simply Π = ∑ a p a and ϑ = ∑ a r a × p a , valid at all post-Newtonian orders (where ( r a , p a ) are the coordinates of particle a ). The energy, as usual, refers to the conserved Hamiltonian, namely E = Hmc 2 , where m = ∑ a m a represents the combined mass parameter of the system. When including the radiative sector to account for losses via the gravitational radiative reaction, the Noether quantities are no longer conserved; at orders 2 . 5PN and above, one must consider flux-balance-type equations (see e.g. [21, 24, 41]) or time-dependent Hamiltonians [23, 42]. In our case we are concerned in studying the conservative sector of the two-body problem.  \nFrom the freedom of performing a Poincar'e transformation, we adopt coordinates of the centre-of-mass (CM) where Π = Γ = Κ = 0. Rothe and Schafer [32] explicitly determine a symplectic transformation to a canonical CM frame which is valid for any Π and Γ . In practice, under CM coordinate conditions ( Π = Γ = Κ = 0), the frame transformation reduces to the Newtonian-like relations p = p 1 = -p 2 , r = r 2 -r 1 (see also [22, 30, 43]). As our starting point, we shall adopt the rescaled CM coordinates:  \np ' = p µ , r ' = r Gm , H ' = Hmc 2 µ , t ' = t Gm , (27)  \nwhere µ = m 1 m 2 /m and ν = µ/m . We also rescale E ' = E /µ and ϑ ' = ϑ / ( Gmµ ). For notational simplicity, we shall henceforth omit the primes. In rescaled coordinates, the reduced matter-only ADM-TT Hamilto-  \n[23, 30]:  \nH 1 ( r , p ) = 1 2 r 2 -1 8 (1 -3 ν ) p 4 -1 2 r ( (3 + ν ) p 2 + ν ( n · p ) 2 ) , (28a) H 2 ( r , p ) = 1 16 ( 1 -5 ν +5 ν 2 ) p 6 + 1 8 r ( ( 5 -20 ν -3 ν 2 ) p 4 -2 ν 2 p 2 ( n · p ) 2 -3 ν 2 ( n · p ) 4 ) + 1 2 r 2 ( (5 + 8 ν ) p 2 +3 ν ( n · p ) 2 ) -1 4 r 3 (1 + 3 ν ), (28b)  \nwith n = r /r and leading zeroth-order term H 0 as specified in Eq. (15).", 'A. Solution and observables': 'We employ the methods outlined in Sec. II to transform the rescaled CM Hamiltonian. Specifically, Eq. (25) parametrically furnishes the normal-form and generator as a function of the perturbation at each PN order. We consider iteratively the two perturbations:  \nP 1 = H 1 , P 2 = H 2 + 1 2 {H 1 + H ∗ 1 , g 1 } . (29)  \nEach perturbation is then expanded into the form (24) by appropriately expressing momenta in terms of the first integrals and the separation, with the aid of the following two expressions:  \np 2 = 2 r -1 L 2 , p 2 r = p 2 -J 2 r 2 . (30)  \nThe procedure leads to the following Hamiltonian:  \nH ∗ ( J, L ) = -1 2 L 2 + 1 c 2 { 15 -ν 8 L 4 -3 JL 3 } + 1 c 4 { 5(2 ν -7) 4 J 3 L 3 -27 2 J 2 L 4 + 3(35 -4 ν ) 4 JL 5 -ν 2 -15 ν +145 16 L 6 } . (31)  \nWe recall that the parameters L = L ( r ∗ , p ∗ ) and J = J ( r ∗ , p ∗ ) are treated as pure Newtonian-order functions of phase-space coordinates. This enables calculations to be treated agnostically in terms of the parameter choice and more readily adapted to a new problem. Simultaneously, Hamiltonian (31) can also be regarded as if it were expressed in Delaunay variables, since the involved transformation is canonical. Unsurprisingly, we find that our expression precisely matches the Delaunay form that can be derived via Hamilton-Jacobi methods (cf. [23, 30]). Nevertheless, Lie methods have the advantage of also providing the generator:  \ng = 1 c 2 { p r r ( J (4 -ν ) -6 L ) 2 JL 2 -6 J arctan ( Lrp r JL + r ) + νp r 2 } + 1 c 4 { p r r 2 ( J 2 ( ν 2 -2 ν ) 8 -3 J 3 ν 4( J + L ) ) + p r r ( 9 J 2( J + L ) + 6 -ν 4 ) + p r ( 5(2 ν -7) 4 J 2 + 3( ν +2) 4 J ( J + L ) + 3 JL + -ν 2 +10 ν -32 8 L 2 ) + p r r ( 5(2 ν -7) 4 J 3 L + 9(3 -ν ) 4 JL 3 -ν +16 8 L 4 ) + ( 10 ν -35 2 J 3 + 15 -6 ν 2 JL 2 ) arctan ( Lrp r JL + r ) } . (32)  \nWe have verified the consistency of the obtained generator, checking that our expressions for H ∗ and g identically satisfy the homologic equations-to order O ( c -6 )-with T g correctly transforming H into H ∗ . In the next steps, we shall derive the integrable dynamics from Hamiltonian (31), and eventually incorporate the short-term oscillatory behaviour to obtain the complete dynamics via application of the Lie transform generated by (32).', '1. Integrable dynamics': "Effectively, the parameter space C ∗ encodes the secular or 'average' behaviour of the system. The dynamics can be easily determined from the transformed Hamiltonian H ∗ via the generalised Hamilton's equations (23b), without explicitly resorting to the variable transformations alluded to in the previous section. We presently provide the expression of six independent osculating Keplerian parameters, following the bar-tilde conventions from Eq. (21). The remaining parameters, such as ¯ v ( t ) and ¯ E ( t ), are obtained via the Keplerian relationships (19)-(20), applied at each moment in time. The full set of parameters can then be used to reconstruct the evolution of the phase-space coordinates ( r ∗ , p ∗ ) under the integrable Hamiltonian H ∗ . As expected, the secular PN perturbations act only on the orbital phase and the periapsis angle:  \n¯ a ( t ) = ¯ a 0 , ¯ e ( t ) = ¯ e 0 , ¯ M ( t ) = ¯ M 0 + ∂ H ∗ ∂L t , (33a)  \n¯ ι ( t ) = ¯ ι 0 , ¯ Ω( t ) = ¯ Ω 0 , ¯ ϖ ( t ) = ¯ ϖ 0 + ∂ H ∗ ∂J t . (33b)  \nThe subscripted zeros refer to initial values obtained from the osculating functions applied to ( r ∗ 0 , p ∗ 0 ). Two of the constants, ¯ a 0 and ¯ e 0 , are directly linked to the Noether quantities E and ϑ , which are also the dynamical invariants of the transformation T g . Namely, this symplectomorphism inherits both the time and the rotation symmetries of the system, as is clear from Eq. (32). Rotational invariance in particular ensures the preservation of ϑ , so that J = ϑ = ¯ J . Moreover, the energy relation E = H ∗ ( ¯ L, J ) can be inverted order-by-order, yielding  \na useful expression:  \n¯ L = 1 √ -2 E + 1 c 2 { 3 J + ν -15 8 √ -2 E } + 1 c 4 { 35 -10 ν 4 J 3 -6 ν -15 2 J E + 3 ν 2 +30 ν +35 128 ( -2 E ) 3 / 2 } . (34)  \nThe values for ¯ a 0 and ¯ e 0 can be deduced from these expressions by inputting ¯ a 0 = ¯ L 2 and ¯ e 0 = (1 -¯ J 2 / ¯ L 2 ) 1 / 2 . Next, we extract the precession and mean orbital fre-  \nquencies from the normal-form Hamiltonian:  \n˙ ¯ M ( t ) = 1 ¯ L 3 + 1 c 2 { 9 ¯ J ¯ L 4 + ν -15 2 ¯ L 5 } + 1 c 4 { 15(7 -2 ν ) 4 ¯ J 3 ¯ L 4 + 54 ¯ J 2 ¯ L 5 + 60 ν -525 4 ¯ J ¯ L 6 + 3 ( ν 2 -15 ν +145 ) 8 ¯ L 7 } , (35a)  \n˙ ¯ ϖ ( t ) = 1 c 2 3 ¯ J 2 ¯ L 3 + 1 c 4 { 15(7 -2 ν ) 4 ¯ J 4 ¯ L 3 + 27 ¯ J 3 ¯ L 4 + 12 ν -105 4 ¯ J 2 ¯ L 5 } . (35b)  \nThese frequencies are associated to well-known observables of the integrable PN dynamics. These are the periapsis-to-periapsis period T and the periapsis advance per orbit ∆ ϕ :  \n˙ ¯ M = Gm ( 2 π T ) , ˙ ¯ ϖ = ∆ ϕ T . (36)  \nFor convenience, we provide the observables T and ∆ ϕ in physical units. Explicit expressions for P and ∆ ϕ have also been derived by Damour, Jaranowski & Schafer [30] in the framework of the Hamilton-Jacobi formalism. It is worth mentioning that their expressions perfectly align with our results when we incorporate Eq. (34).", '2. Complete dynamics of the original system': "We may now return to the original coordinate system C . Applying T g to each Keplerian parameter, we extract their full temporal evolution as a function of the secular variables derived in the previous section [cf. Eq. (22)].  \nAfter some lengthy bracket computations, we obtain the expressions that follow, strictly valid at O ( c -6 ):  \n˜ a = ¯ a + 6 ∑ k =0 A ( a ) k cos( k ¯ v ), (37a)  \n˜ e = ¯ e + 6 ∑ k =0 A ( e ) k cos( k ¯ v ), (37b)  \n˜ v = ¯ v + 6 ∑ k =1 A ( v ) k sin( k ¯ v ) + (¯ v -¯ E ) ( 2 ∑ k =0 B ( v ) k cos( k ¯ v ) ) , (37c)  \n˜ ϖ = ¯ ϖ + 6 ∑ k =1 A ( ϖ ) k sin( k ¯ v ) + B ( ϖ ) (¯ v -¯ M ). (37d)  \nThe Kepler parameters ˜ ι = ¯ ι and ˜ Ω = ¯ Ω are identically unchanged by the T g . Each of the above coefficients is time-independent; they are reported in Appendix B.  \nWe remark that certain coefficients [Eqs. (37b-37d)] contain coordinate singularities in the quasi-circular regime (¯ e → 0 or equivalently ¯ J → ¯ L ). Specifically, their Laurent expansions contain single or double poles in ¯ e = 0, and for small enough eccentricity (¯ e ≲ 1 /c 2 ) the residuals are unbounded. This is an artefact of expanding the post-Newtonian osculating parametrisation and can be bypassed by adopting 'Poincar'e-like' elements:  \nz = e exp(i ϖ ), λ = v + ϖ , ζ = sin( ι/ 2) exp(iΩ), (38)  \nWe have determined the full temporal evolution of the aforementioned 'regularised' elements:  \n˜ z = ¯ z +e i ¯ ϖ { 6 ∑ k = -6 ( A ( z ) k +i B ( z ) k ( ¯ v -¯ E ) + C ( z ) k (¯ v -¯ E ) 2 ) e i k ( ¯ ϖ -¯ λ ) + 2 ∑ k = -2 ( D ( z ) k +i E ( z ) k ( ¯ v -¯ E ) ) e i k ( ¯ ϖ -¯ λ E ) } , (39a)  \n˜ λ = ¯ λ + 4 ∑ k =1 i A ( λ ) k ( ¯ z k -(¯ z † ) k ) cos( k ¯ λ ) + B ( λ ) k ( ¯ z k +(¯ z † ) k ) sin( k ¯ λ )  \n+(¯ v -¯ E ) ( 2 ∑ k =0 C ( λ ) k ( ¯ z k +(¯ z † ) k ) cos( k ¯ λ ) + i D ( λ ) k ( ¯ z k -(¯ z † ) k ) sin( k ¯ λ ) ) + B ( ϖ ) (¯ v -¯ M ), (39b)  \n˜ ζ = ¯ ζ , (39c)  \nwith λ E = E + ϖ .  \nFor convenience, the eccentricity can also be extracted via the exact expression instead of the truncated series:  \n˜ e = √ 1 -¯ J 2 / ˜ a , (40)  \nwith ˜ a given by Eq. (37a). 3", '3. Numerical example': "The analytical expressions provided in the previous section solve the ADM Hamiltonian (28) at the O ( c -6 ) level. We shall now illustrate their validity numerically, for an unequal-mass moderately-relativistic binary in both quasi-circular and eccentric orbit configurations. For each scenario, we integrate H by employing an 8thorder implicit Runge-Kutta scheme. Integration is performed in polar coordinates ( r, λ, p r , J ), where Hamilton's equations reduce to (see for instance [31]; rotational invariance is required):  \n˙ r = ∂ H ∂p r , ˙ p r = -∂ H ∂r , (41a)  \n˙ λ = ∂ H ∂J , ˙ J = -∂ H ∂λ = 0, (41b)  \nwhere an implicit change of variables is performed with the aid of Eqs. (30). Figure 1 depicts the phase error by which the perturbative solution (in C coordinate space) deviates from the numerical integrator in each scenario. As we might anticipate, the analytical residue manifests close to 2PN marks below that of a standard Keplerian solution. This residue arises from the truncation of higher-order terms in the Lie expansion of T g , as portrayed in Fig. 2. Specfically, the secular Kepler parameters ¯ q ( t ) are obtained from a trunction of H ∗ to its integrable terms, while a second truncation occurs in the transformation ˜ q ( t ) = T g (¯ q )( t ). Consequently, this residue manifests as terms of order O ( c -6 ) that are linear in time for ˜ ϖ ( t ) and ˜ M ( t ) and periodic for the other variables, which is consistent with a 2PN solution.  \nFIG. 1. Evolution of the phase residue ( λ ) of the truncated analytical solution (in solid red ) when compared to numerical integration evolving under the 2PN ADM Hamiltonian [Eq. (28)], for a moderately-relativistic binary of velocity p/c ∼ 10 -3 and mass-ratio ν = 2 / 9. Two eccentricity configurations are considered: on the left , a quasi-circular orbit ( e = 0 . 01); on the right , a highly eccentric orbit ( e = 0 . 8). The solid blue benchmark denotes the error of a purely Keplerian analytical solution. The numerical accuracy baseline at 10 -15 is represented in dashed black , estimated from the residue of the Keplerian solution when c →∞ . We have also verified that the conservation of the ADM constants E and J respect this accuracy level.  \n<!-- image -->  \nFIG. 2. Schematics of the Lie canonical framework. The leftmost column conceptualises a complete relativistic Hamiltonian system with an infinite series of PN terms. The middle column represents a truncation at O ( c -2 K ). The rightmost column corresponds to the secular system with Hamiltonian H ∗ in normal form, obtained by applying the Lie transformation T g to H . Indeed, the red pathway illustrates the key transformation steps of the Lie approach; meanwhile, the blue arrow represents direct integration of H , a nonlinear process with a potentially infinite amount of terms. In practical computations, each application of the Lie transform T g is truncated to order O ( c -2 K ), leading to a computation of ˜ q ( t ) that is formally accurate to the same order.  \n<!-- image -->  \nℓ  \n=0", 'IV. CONCLUSIONS': 'We have provided a canonical Lie framework for analyzing the perturbed two-body problem under generic conditions of time-independence and rotational invariance. Our framework allows for the computation of the complete coordinate-time evolution in the PN conservative sector, using centre-of-mass coordinates C . It is also applicable to more general perturbative two-body systems under similar conditions. The formalism can theoretically be applied to perturbations of arbitrary order, although the residue is expected to grow quickly at very high orders (see e.g. [36]). We provide a generic family of solutions [Eq. (25)] for systems within the aforementioned conditions. Accordingly, we specify an equivalent Hamiltonian in a new set of phase-space coordinates C ∗ , which is integrable at given order; as well as the corresponding Lie transformation T g : C ∗ → C , enabling the conversion between both coordinate systems. From the boundedness of the obtained Lie generator, it is clear that resonances are absent from the class of systems which  \nhave been studied. We have also applied our framework to the local, conservative ADM sector of the non-spinning point-mass binary problem, recovering classical results.  \nGoing forward, our formalism may be broadened to other contexts, such as spinning pole-dipole models [23, 40, 44, 45] or time-dependent systems, including the ADM dissipative sector and non-local terms at 4PN [23]. Notably, the radiation-reaction term at 2 . 5PN may be treated by either doubling the phase-space variables [23, 46-48] or incorporating explicit time-dependence. In general, the presence of non-autonomous perturbative terms no longer leads to the preservation of the value of the conservative Hamiltonian by the Lie transform T g ; regardlessly, similar homologic relations to Eq. (8) can be derived, but in extended phase-space. Spinning systems may benefit from the phase-space description of the Keplerian variables specified in Sec. II B; when properly normalised, the angular momentum and eccentricity vectors form the algebra of SO (3) × SO (3) ∼ SO (4), potentially expressible in spinor form (see [49] for an example). We conclude by remarking that in instances where  \nresonances are present, order-by-order execution of the normal-form process is recommended (see e.g. [36] for details).', 'SUPPLEMENTARY INFORMATION': "We provide source files in Mathematica language, organised as follows:  \n- · A Wolfram Language package file 'KeplerPert.wl' containing methods for computing the canonical Lie transform and for solving autonomous rotationally-invariant two-body homologic equations;\n- · A Wolfram Notebook file 'ADM 2PN.nb' exemplifying the package's use applied to the 2PN ADM Hamiltonian;\n- · An auxiliary package 'PoissonBracket.wl' implementing some basic properties of Poisson bracket algebra;\n- · An auxiliary package 'Utils.wl' with a collection of generic helper functions.", 'ACKNOWLEDGMENTS': "C.A. acknowledges the joint finantial support of Centre National d' ' Etudes Spatiales (CNES) and ' Ecole Doctorale Astronomie et Astrophysique d'Ile de France (ED127 AAIF). This work was also supported by the Programme National GRAM, by PNPS (CNRS/INSU), by INP and IN2P3 co-funded by CNES, and by CNES LISA grants at CEA/IRFU. We are very grateful to G. Faye for the constructive comments and discussions.", 'Appendix A: Flow integrals': 'In this section we search for explicit expressions to the basis elements ξ k and h k , which parametrise the solutions to the homologic equations. We recall their definitions:  \nξ k = ∫ 1 r k ( d θ d t ) -1 kep d θ , h k = 〈 1 r k 〉 ,  \nwith the integration constant chosen such that ⟨ ξ k ⟩ = 0. It is immediately clear from the definition that h 0 = 1. Constraining θ to a Kepler angle further leads to h k = ( n 0 / 2 π ) ( ξ k ( π ) -ξ k ( -π ) ) , where n 0 = 1 /L 3 . For higher orders of ξ k and h k we shall split the computation into two distinct cases. Namely, using the definitions in Sec. II B, we work out the relationships that follow:  \n- · For k = 1, identifying θ = E , we have:  \n( d E d t ) kep = 1 Lr , E = M + rp r L . (A1)  \n- · For k ≥ 2, identifying θ = v , we have instead:  \n( d v d t ) kep = J r 2 , v = M + rp r L +2arctan ( Lrp r JL + r ) . (A2)  \nThe integrals for k = 1 , 2 are trivial and yield:  \nξ 1 = LE , h 1 = Ln 0 . (A3)  \nξ 2 = v J , h 2 = n 0 J . (A4)  \nThe evaluation of higher values of k involves a little more work. Namely, we consider, for integer k ≥ 3:  \nξ k = 1 J 2 k -3 k -2 ∑ j =0 ( k -2 j ) e j ∫ cos j v d v . (A5)  \nFortunately, the integrals of powers of cosine on the right K α = ∫ cos α v d v can be determined analytically. They follow from the recurrence relationship:  \nK α = ( α -1 α ) K α -2 + 1 α sin v cos α -1 v , K 0 = v . (A6)  \nWe verify via induction that the recurrence is satisfied by:  \nK α = sin v ∑ β ∈ A β !! α !! ( α -1)!! ( β -1)!! 1 β cos β -1 v + α e ( α -1)!! α !! v ,  \n(A7)  \nwhere ( · )!! is the double factorial, α e := 1 -( α mod 2), and the sum is performed over the set A = { α e +1 , α e + 3 , . . . , α -2 , α } . The last term is only present when α is even, as evidenced by the presence of the α e factor. Putting everything together and rearranging:  \nξ k = p r J 2 k -4 k -2 ∑ α =0 ( k -2 α ) × ∑ β ∈ A β !! α !! ( α -1)!! ( β -1)!! e α -β β ( J 2 r -1 ) β -1 + v J 2 k -3 ⌊ k 2 ⌋-1 ∑ α =0 ( k -2 2 α )( 2 α α ) ( e 2 ) 2 α .  \nWe can now extract h k :  \nh k = n 0 J 2 k -3 ⌊ k 2 ⌋-1 ∑ α =0 ( k -2 2 α )( 2 α α ) ( e 2 ) 2 α , (A8)  \nvalid for integer k ≥ 2, and where ⌊ x ⌋ is the largest integer smaller or equal to x . The angular relationships provided [Eqs. (A1, A2)] allow us to express g ℓ = ∑ k A ℓ,k ( ξ k -( M/n 0 ) h k ) solely as a function of r , p r and p . Finally, we can confirm via parity arguments that the average of each term does not contribute to g ℓ :  \n〈 ξ k -M n 0 h k 〉 = 0. (A9)', 'Appendix B: Coefficients of the Keplerian coordinate transformation': 'We provide the coefficients of the 2PN ADM transformations for the Keplerian parameters from C ∗ to C as detailed in Sec. III A 2:', '· coefficients for a:': 'A ( a ) 0 = 1 c 2 { 4 -ν -6 ¯ L ¯ J + ¯ L 2 (2 ν -7) ¯ J 2 + ¯ L 4 (9 -ν ) ¯ J 4 + 1 c 2 ( 7 ( 7 ν 2 -96 ν +360 ) ¯ L 6 16 ¯ J 8 -( 207 ν 2 -1876 ν +4360 ) ¯ L 4 32 ¯ J 6 + ( 17 ν 2 -98 ν +208 ) ¯ L 2 4 ¯ J 4 + 3(7 ν -90) ¯ L 3 4 ¯ J 5 -( ν -32) ¯ L 4 ¯ J 3 -3(2 ν -5) 2 ¯ J ¯ L -27 ν 2 -196 ν +40 32 ¯ J 2 + -ν -16 4 ¯ L 2 ) } (B1)  \nA ( a ) 1 = ¯ e c 2 { ¯ L 2 (7 ν -16) 4 ¯ J 2 + ¯ L 4 (48 -7 ν ) 4 ¯ J 4 + 1 c 2 ( -24 ¯ J ( ¯ J + ¯ L ) + 640 + 124 ν -13 ν 2 32 ¯ J 2 -3(9 ν +40) ¯ L 8 ¯ J 3 (B2) + ( 75 ν 2 -426 ν +1306 ) ¯ L 2 16 ¯ J 4 + 3(35 ν -312) ¯ L 3 8 ¯ J 5 -( 313 ν 2 -2296 ν +3684 ) ¯ L 4 32 ¯ J 6 + ( 11 ν 2 -147 ν +504 ) ¯ L 6 2 ¯ J 8 ) } ,  \nA ( a ) 2 = ¯ e 2 c 2 { ¯ L 4 (3 -ν ) ¯ J 4 + 1 c 2 ( 6( ν -7) ¯ J ( ¯ J + ¯ L ) -6( ν -7) ¯ J 2 + 6( ν -7) ¯ L ¯ J 3 + ( 67 ν 2 -632 ν +2120 ) ¯ L 2 32 ¯ J 4 + 3(5 ν -24) ¯ L 3 ¯ J 5 -( 105 ν 2 -532 ν +468 ) ¯ L 4 16 ¯ J 6 + ( 127 ν 2 -1536 ν +4032 ) ¯ L 6 32 ¯ J 8 ) } (B3)  \nA ( a ) 3 = ¯ e c 2 { ¯ L 2 ν 4 ¯ J 2 -¯ L 4 ν 4 ¯ J 4 + 1 c 2 ( -29 ν 2 +200 ν -256 64 ¯ J 2 -3(11 ν -16) ¯ L 16 ¯ J 3 + ( 121 ν 2 -710 ν +1116 ) ¯ L 2 32 ¯ J 4 + 3(55 ν -144) ¯ L 3 16 ¯ J 5 -( 357 ν 2 -2084 ν +2744 ) ¯ L 4 64 ¯ J 6 + 3 ( 3 ν 2 -29 ν +48 ) ¯ L 6 4 ¯ J 8 ) } (B4)  \nA ( a ) 4 = ¯ e 2 c 4 { (3 ν -4)(11 ν -48) ¯ L 2 32 ¯ J 4 + 3(5 ν -6) ¯ L 3 4 ¯ J 5 -( 63 ν 2 -260 ν +192 ) ¯ L 4 32 ¯ J 6 + 3 ( 5 ν 2 -32 ν +24 ) ¯ L 6 16 ¯ J 8 } (B5)  \nA ( a ) 5 = ¯ e 2 c 4 { 3 ν (16 -3 ν ) 64 ¯ J 2 -9 ν ¯ L 16 ¯ J 3 + ν (17 ν -54) ¯ L 2 32 ¯ J 4 + 9 ν ¯ L 3 16 ¯ J 5 -ν (41 ν -108) ¯ L 4 64 ¯ J 6 + ( ν -3) ν ¯ L 6 4 ¯ J 8 } (B6)  \nA ( a ) 6 = ¯ e 2 c 4 { ν 2 ¯ L 2 32 ¯ J 4 -ν 2 ¯ L 4 16 ¯ J 6 + ν 2 ¯ L 6 32 ¯ J 8 } , (B7)', '· Coefficients for e :': 'A ( e ) 0 = ¯ e c 2 { 3 ¯ J ¯ L -3 ¯ J ( ¯ J + ¯ L ) -ν -9 2 ¯ J 2 + ν -4 2 ¯ L 2 + 1 c 2 ( -181 ν 2 +1144 ν -496 128 ¯ J 2 ¯ L 2 + 2 ν +13 2 ¯ J 3 ( ¯ J + ¯ L ) + 94 -35 ν 8 ¯ J 3 ¯ L + 3(8 ν -29) 4 ¯ J ¯ L 3 + 2 ¯ J 3 ( ¯ L -¯ J ) -9 2 ¯ J 2 ( ¯ J + ¯ L ) 2 + 157 ν 2 -2296 ν +5080 128 ¯ J 4 + 3 ν 2 -23 ν +64 8 ¯ L 4 ) } , (B8)  \nA ( e ) 1 = 1 c 2 { 7 ν -16 8 ¯ L 2 -7 ν -48 8 ¯ J 2 + 1 c 2 ( -89 ν 2 +724 ν -1174 32 ¯ J 2 ¯ L 2 -9(7 ν -24) 16 ¯ J 3 ¯ L + 9(11 ν -24) 16 ¯ J ¯ L 3 + 107 ν 2 -1596 ν +3516 64 ¯ J 4 + 71 ν 2 -404 ν +640 64 ¯ L 4 ) } , (B9)  \nA ( e ) 2 = ¯ e c 2 { -( ν -3) 2 ¯ J 2 + 1 c 2 ( -113 ν 2 +352 ν +32 256 ¯ J 2 ¯ L 2 + 3 ν -29 2 ¯ J 3 ( ¯ J + ¯ L ) -3 2 ¯ J 3 ¯ L -2 ¯ J 3 ( ¯ L -¯ J ) + 49 ν 2 -1632 ν +3744 256 ¯ J 4 ) } , (B10)  \nA ( e ) 3 = ¯ e 2 c 2 { ν 8 ¯ L 2 -ν 8 ¯ J 2 + 1 c 2 ( 39 ν 2 -340 ν +828 64 ¯ J 2 ¯ L 2 + 3( ν +16) 32 ¯ J ¯ L 3 + 9(13 ν -48) 32 ¯ J 3 ¯ L -73 ν 2 -224 ν +248 128 ¯ J 4 + -5 ν 2 +104 ν -256 128 ¯ L 4 ) } , (B11)  \nA ( e ) 4 = ¯ e c 4 { 45 ν 2 -328 ν +384 128 ¯ J 2 ¯ L 2 + 3(5 ν -6) 8 ¯ J 3 ¯ L -53 ν 2 -168 ν +168 128 ¯ J 4 } , (B12)  \nA ( e ) 5 = 1 c 4 { 9 ν 32 ¯ J 3 ¯ L + (11 ν -36) ν 64 ¯ J 2 ¯ L 2 -9 ν 32 ¯ J ¯ L 3 -(13 ν -24) ν 128 ¯ J 4 -3(3 ν -16) ν 128 ¯ L 4 } , (B13)  \nA ( e ) 6 = ¯ e 3 c 4 { -ν 2 256 ¯ J 4 } , (B14)', '· Coefficients for v :': "A ( v ) 1 = 1 c 2 { -2 ¯ J ( ¯ L -¯ J ) -4 ¯ J ( ¯ J + ¯ L ) + 13 ν -32 8 ¯ J 2 + 1 c 2 ( 43 ν 2 -260 ν +1152 64 ¯ J 2 ¯ L 2 + 9(23 ν -72) 16 ¯ J 3 ¯ L + 43 ν -74 16 ¯ J 3 ( ¯ L -¯ J ) + 178 -41 ν 16 ¯ J 3 ( ¯ J + ¯ L ) -3 ¯ J 2 ( ¯ J + ¯ L ) 2 + ν 2 +288 ν -68 64 ¯ J 4 ) } , (B15)  \nA ( v ) 2 = 1 c 2 { -3 2 ¯ J ¯ L + ν -2 ¯ J 2 + 4 -ν 2 ¯ L 2 + 1 c 2 ( -53 ν 2 +688 ν -1568 128 ¯ J 2 ¯ L 2 -3(14 ν -47) 8 ¯ J ¯ L 3 + 79 -8 ν 8 ¯ J 3 ¯ L + 4 ¯ J 3 ( ¯ L -¯ J ) -10 ¯ J 3 ( ¯ J + ¯ L ) + 117 ν 2 -80 ν +608 128 ¯ J 4 + -ν 2 +7 ν -32 4 ¯ L 4 ) } , (B16)  \nA ( v ) 3 = ¯ e c 2 { ν 8 ¯ J 2 + 1 c 2 ( -67 ν 2 +600 ν -1280 128 ¯ J 2 ¯ L 2 + 3 -ν ¯ J 3 ( ¯ L -¯ J ) -3(41 ν -176) 32 ¯ J 3 ¯ L + -19 ν -36 8 ¯ J 3 ( ¯ J + ¯ L ) + 163 ν 2 -580 ν +792 128 ¯ J 4 ) } , (B17)  \nA ( v ) 4 = 1 c 4 { -21 ν 2 +126 ν -92 32 ¯ J 2 ¯ L 2 + 3( ν -4) 4 ¯ J ¯ L 3 -3(7 ν -8) 8 ¯ J 3 ¯ L + ( ν -2)(21 ν -32) 32 ¯ J 4 + ( ν -4) 2 8 ¯ L 4 } , (B18)  \nA ( v ) 5 = ¯ e c 4 { -9 ν 32 ¯ J 3 ¯ L -3(3 ν -16) ν 128 ¯ J 2 ¯ L 2 + (17 ν -36) ν 128 ¯ J 4 } , (B19)  \nA ( v ) 6 = ¯ e 2 c 4 { ν 2 128 ¯ J 4 } , (B20)  \nB ( v ) 0 = 1 c 4 { 9(2 ν -5) 4 ¯ J 4 -3(2 ν -5) 4 ¯ J 2 ¯ L 2 } , (B21)  \nB ( v ) 1 = ¯ e c 4 { 3(2 ν -5) ¯ J 4 } , (B22)  \nB ( v ) 2 = ¯ e 2 c 4 { 3(2 ν -5) 4 ¯ J 4 } . (B23)  \n- · Coefficients for ϖ :  \nA ( ϖ ) 1 = ¯ e c 2 { 2 ¯ J ( ¯ L -¯ J ) -2 ¯ J ( ¯ J + ¯ L ) -5 ν -48 8 ¯ J 2 + 1 c 2 ( -7 ν 2 +28 ν -128 64 ¯ J 2 ¯ L 2 -9(3 ν -8) 16 ¯ J 3 ¯ L + 74 -43 ν 16 ¯ J 3 ( ¯ L -¯ J ) + 25 ν -362 16 ¯ J 3 ( ¯ J + ¯ L ) -6 ¯ J 2 ( ¯ J + ¯ L ) 2 + 11 ν 2 -712 ν +2100 64 ¯ J 4 ) } , (B24)  \nA ( ϖ ) 2 = 1 c 2 { -ν -3 2 ¯ J 2 + 1 c 2 ( 21 ν 2 -304 ν +1296 128 ¯ J 2 ¯ L 2 + 3( ν -14) 4 ¯ J 3 ¯ L -4 ¯ J 3 ( ¯ L -¯ J ) -8 ¯ J 3 ( ¯ J + ¯ L ) -53 ν 2 -176 ν +336 128 ¯ J 4 ) } , (B25)  \nA ( ϖ ) 3 = ¯ e c 2 { -ν 8 ¯ J 2 + 1 c 2 ( 11 ν 2 -104 ν +256 128 ¯ J 2 ¯ L 2 + 3( ν -16) 32 ¯ J 3 ¯ L + ν -3 ¯ J 3 ( ¯ L -¯ J ) + 19 ν -72 8 ¯ J 3 ( ¯ J + ¯ L ) -107 ν 2 -564 ν +824 128 ¯ J 4 ) } , (B26)  \nA ( ϖ ) 4 = 1 c 4 { 13 ν 2 -86 ν +96 32 ¯ J 2 ¯ L 2 + 3(5 ν -6) 8 ¯ J 3 ¯ L -17 ν 2 -66 ν +60 32 ¯ J 4 } , (B27)  \nA ( ϖ ) 5 = ¯ e c 4 { 9 ν 32 ¯ J 3 ¯ L + 3(3 ν -16) ν 128 ¯ J 2 ¯ L 2 -(17 ν -36) ν 128 ¯ J 4 } , (B28)  \nA ( ϖ ) 6 = ¯ e 2 c 4 { -ν 2 128 ¯ J 4 } , (B29)  \nB ( ϖ ) = 1 c 2 { 3 ¯ J 2 + 1 c 2 ( 3(2 ν -5) 4 ¯ J 2 ¯ L 2 -15(2 ν -7) 4 ¯ J 4 ) } . (B30)  \n- · The non-zero coefficients for z are:  \nA ( z ) 0 = ¯ e c 2 { 3 ¯ J ¯ L -3 ¯ J ( ¯ J + ¯ L ) -ν -9 2 ¯ J 2 + ν -4 2 ¯ L 2 + 1 c 2 ( -21 ν 2 +133 ν -10 16 ¯ J 2 ¯ L 2 -29( ν -2) 8 ¯ J 3 ¯ L + 3(8 ν -29) 4 ¯ J ¯ L 3 + 2 ν +17 2 ¯ J 3 ( ¯ J + ¯ L ) -9 2 ¯ J 2 ( ¯ J + ¯ L ) 2 + 17 ν 2 -251 ν +446 16 ¯ J 4 + 3 ν 2 -23 ν +64 8 ¯ L 4 ) } (B31)  \nA ( z ) 1 = 1 c 2 { -ν ¯ e 2 8 ¯ J 2 + 1 c 2 ( -13 ν 2 +63 ν -149 16 ¯ J 2 ¯ L 2 -3( ν -72) 16 ¯ J 3 ¯ L + 3(5 ν -8) 16 ¯ J ¯ L 3 + 15 ν 2 -92 ν -150 32 ¯ J 4 + 11 ν 2 -58 ν +64 32 ¯ L 4 ) } (B32)  \nA ( z ) -1 = 1 c 2 { 3(8 -ν )¯ e 2 4 ¯ J 2 + 1 c 2 ( -3 ( 19 ν 2 -188 ν +276 ) 32 ¯ J 2 ¯ L 2 -9(3 ν -8) 8 ¯ J 3 ¯ L + 3(17 ν -48) 8 ¯ J ¯ L 3 + 65 ν 2 -1280 ν +3528 64 ¯ J 4 + 49 ν 2 -288 ν +512 64 ¯ L 4 ) } , (B33)  \nA ( z ) 2 = ¯ e c 4 { -ν 2 +8 ν -48 16 ¯ J 2 ¯ L 2 + 3( ν +6) 8 ¯ J 3 ¯ L + ν 2 -14 ν +12 16 ¯ J 4 } (B34)  \nA ( z ) -2 = ¯ e c 2 { -ν -3 2 ¯ J 2 + 1 c 2 ( -7 ν 2 +22 ν +34 16 ¯ J 2 ¯ L 2 + 3( ν -11) 2 ¯ J 3 ( ¯ J + ¯ L ) -3(3 ν -2) 8 ¯ J 3 ¯ L + 3 ν 2 -112 ν +366 16 ¯ J 4 ) } (B35)  \nA ( z ) 3 = 1 c 4 { 3 ν 32 ¯ J 3 ¯ L + ( ν -10) ν 64 ¯ J 2 ¯ L 2 -3 ν 32 ¯ J ¯ L 3 -( ν -4) ν 128 ¯ J 4 -( ν -16) ν 128 ¯ L 4 } (B36)  \nA ( z ) -3 = 1 c 2 { ν 8 ¯ L 2 -ν 8 ¯ J 2 + 1 c 2 ( 7 ν 2 -67 ν +183 16 ¯ J 2 ¯ L 2 + 3(3 ν +8) 16 ¯ J ¯ L 3 + 3(17 ν -72) 16 ¯ J 3 ¯ L -13 ν 2 +8 ν -82 32 ¯ J 4 + -ν 2 +22 ν -64 32 ¯ L 4 ) } (B37)  \nA ( z ) -4 = ¯ e c 2 { 5 ν 2 -39 ν +48 16 ¯ J 2 ¯ L 2 + 3(5 ν -6) 8 ¯ J 3 ¯ L -5 ν 2 -9 ν +12 16 ¯ J 4 } (B38)  \nA ( z ) -5 = 1 c 4 { 9 ν 32 ¯ J 3 ¯ L + 3(3 ν -10) ν 64 ¯ J 2 ¯ L 2 -9 ν 32 ¯ J ¯ L 3 -3(3 ν -4) ν 128 ¯ J 4 -3(3 ν -16) ν 128 ¯ L 4 } (B39)  \nB ( z ) 0 = ¯ e c 2 { 3 ¯ J 2 + 1 c 2 ( 3(4 ν -13) 4 ¯ J 2 ¯ L 2 + 9 ¯ J 3 ¯ L -9 ¯ J 3 ( ¯ J + ¯ L ) -3(12 ν -53) 4 ¯ J 4 ) } , (B40)  \nB ( z ) 1 = ¯ e 2 c 4 { -3 ν 8 ¯ J 4 } , (B41)  \nB ( z ) -1 = 1 c 4 { 3(3 ν -8) 4 ¯ J 2 ¯ L 2 -9( ν -8) 4 ¯ J 4 } (B42)  \nB ( z ) -2 = ¯ e c 4 { -3( ν -3) 2 ¯ J 4 } (B43)  \nB ( z ) -3 = ¯ e 2 c 4 { -3 ν 8 ¯ J 4 } (B44)  \nE  \n(  \nz  \n)  \n1  \n=  \nC ( z ) 0 = 1 c 4 { -9 2 ¯ J 4 } (B45)  \nD ( z ) 1 = 1 c 2 { -3¯ e 2 2 ¯ J 2 + 1 c 2 ( -3( ν -10) 4 ¯ J ¯ L 3 -3(3 ν -10) 2 ¯ J 2 ¯ L 2 -9 2 ¯ J 3 ¯ L + 15(2 ν -7) 8 ¯ J 4 + 3(4 ν -13) 8 ¯ L 4 ) } , (B46)  \nD ( z ) -1 = 1 c 2 { 3¯ e 2 2 ¯ J 2 + 1 c 2 ( -3( ν +2) 4 ¯ J ¯ L 3 + 3(3 ν -10) 2 ¯ J 2 ¯ L 2 -9 2 ¯ J 3 ¯ L -15(2 ν -7) 8 ¯ J 4 -3(4 ν -13) 8 ¯ L 4 ) } , (B47)  \nD ( z ) 2 = ¯ e 3 c 4 { 9 8 ¯ J 4 } , (B48)  \nD ( z ) -2 = ¯ e 3 c 4 { 9 8 ¯ J 4 } , (B49)  \n1  \n2  \nc  \n-  \n2  \n9  \n¯  \nJ  \n4  \n}  \n,  \n(B50)  \nE ( z ) -1 = 1 c 2 { 9 2 ¯ J 4 } (B51)  \nThe remaining coefficients are identically zero.  \n- · The non-zero coefficients for λ are:  \nA ( λ ) 1 = 1 c 2 { ( ν +2) 2 J 2 -3 J ( J + L ) + 1 c 2 ( ( 3 ν 2 -106 ν +508 ) 32 J 4 -(2 ν +23) 4 J 3 ( J + L ) + 9(5 ν -16) 8 J 3 L + ( 9 ν 2 -58 ν +256 ) 32 J 2 L 2 -9 2 J 2 ( J + L ) 2 ) } (B52)  \nA ( λ ) 2 = 1 c 2 { ( ν -1) 4 J 2 -3 4 J ( J + L ) + 1 c 2 ( ( 4 ν 2 +6 ν +17 ) 16 J 4 -(11 ν +2) 4 J 3 ( J + L ) + 3(14 ν -47) 16 J 3 L + ( ν 2 -7 ν +32 ) 8 J 2 L 2 + 9 2 J 2 ( J + L ) 2 ) } (B53)  \nA ( λ ) 3 = 1 c 4 { ( 7 ν 2 -2 ν -4 ) 32 J 4 -3(5 ν -2) 8 J 3 ( J + L ) + 27 8 J 2 ( J + L ) 2 } (B54)  \nA ( λ ) 4 = 1 c 4 { ( ν -1) 2 16 J 4 -3( ν -1) 8 J 3 ( J + L ) + 9 16 J 2 ( J + L ) 2 } (B55)  \nB ( λ ) 1 = 1 c 2 { ν +2 2 J 2 -3 J ( J + L ) + 1 c 2 ( 3 ν 2 -106 ν +508 32 J 4 + 9(5 ν -16) 8 J 3 L + -2 ν -23 4 J 3 ( J + L ) + 9 ν 2 -58 ν +256 32 J 2 L 2 -9 2 J 2 ( J + L ) 2 ) } (B56)  \n{  \nB ( λ ) 2 = 1 c 2 { ν -1 4 J 2 -3 4 J ( J + L ) + 1 c 2 ( 4 ν 2 +6 ν +17 16 J 4 + 3(14 ν -47) 16 J 3 L + -11 ν -2 4 J 3 ( J + L ) + ν 2 -7 ν +32 8 J 2 L 2 + 9 2 J 2 ( J + L ) 2 ) } (B57)  \nB ( λ ) 3 = 1 c 4 { 7 ν 2 -2 ν -4 32 J 4 -3(5 ν -2) 8 J 3 ( J + L ) + 27 8 J 2 ( J + L ) 2 } (B58)  \nB ( λ ) 4 = 1 c 4 { ( ν -1) 2 16 J 4 -3( ν -1) 8 J 3 ( J + L ) + 9 16 J 2 ( J + L ) 2 } (B59)  \nC ( λ ) 0 = 1 c 4 { 15(2 ν -7) 4 J 4 -3( ν -5) J 4 -3(2 ν -5) 4 J 2 L 2 } (B60)  \nC ( λ ) 1 = C ( λ ) 2 = 1 c 4 { 3(2 ν -5) 2 J 4 } (B61)  \nD ( λ ) 1 = D ( λ ) 2 = 1 c 4 { -3(2 ν -5) 2 J 4 } (B62)  \nThe remaining coefficients are identically zero.  \n- [1] A. Einstein et al. , Annalen Phys 49 , 769 (1916).\n- [2] H. Lorentz and J. Droste, in Collected Papers: Volume V (Springer, Dordrecht, 1937) pp. 330-355.\n- [3] T. Levi-Civita, American Journal of Mathematics 59 , 225 (1937).\n- [4] A. Einstein, L. Infeld, and B. Hoffmann, Annals of mathematics , 65 (1938).\n- [5] C. M. Will, Proceedings of the National Academy of Sciences 108 , 5938-5945 (2011).\n- [6] R. A. Hulse and J. H. Taylor, The Astrophysical Journal Letters 195 , L51 (1975).\n- [7] B. P. Abbott et al. , Phys. Rev. Lett. 116 , 061102 (2016).\n- [8] T. L. S. Collaboration and J. Aasi, Class. Quantum Gravity 32 , 074001 (2015).\n- [9] D. Bersanetti, B. Patricelli, O. J. Piccinni, F. Piergiovanni, F. Salemi, and V. Sequino, Universe 7 , 322 (2021).\n- [10] T. Akutsu et al. , Nature Astronomy 3 , 35 (2019).\n- [11] R. Manchester, Classical and Quantum Gravity 30 , 224010 (2013).\n- [12] P. Amaro-Seoane et al. , Laser interferometer space antenna (2017).\n- [13] M. Maggiore et al. , Journal of Cosmology and Astroparticle Physics 2020 (03), 050.\n- [14] D. Reitze et al. , Cosmic explorer: The u.s. contribution to gravitational-wave astronomy beyond ligo (2019).\n- [15] A. Ross, Physical Review D 85 , 125033 (2012).\n- [16] L. Blanchet, Living Reviews in Relativity 17 , 10.12942/lrr-2014-2 (2014).\n- [17] E. Poisson and C. M. Will, Gravity: Newtonian, PostNewtonian, Relativistic (Cambridge University Press, Cambridge, 2014).\n- [18] W. D. Goldberger and I. Z. Rothstein, Physical Review D 73 , 104029 (2006).\n- [19] R. A. Porto, Physics Reports 633 , 1 (2016).\n- [20] M. Levi, Reports on Progress in Physics 83 , 075901 (2020).\n- [21] Q. Henry, F. Larrouturou, and C. Le Poncin-Lafitte, Physical Review D 108 , 024020 (2023).\n- [22] P. Jaranowski and G. Schafer, Physical Review D 57 , 7274 (1998).\n- [23] G. Schafer and P. Jaranowski, Living Reviews in Relativity 27 , 10.1007/s41114-024-00048-7 (2024).\n- [24] V. C. de Andrade, L. Blanchet, and G. Faye, Classical and Quantum Gravity 18 , 753 (2001).\n- [25] T. Damour and N. Deruelle, Annales De L Institut Henri Poincare-physique Theorique 43 , 107 (1985).\n- [26] R.-M. Memmesheimer, A. Gopakumar, and G. Schafer, Physical Review D 70 , 104011 (2004).\n- [27] G. Cho, S. Tanay, A. Gopakumar, and H. M. Lee, Phys. Rev. D 105 , 064010 (2022).\n- [28] C. W. Lincoln and C. M. Will, Phys. Rev. D 42 , 1123 (1990).\n- [29] A. Tucker and C. M. Will, Physical Review D 104 , 104023 (2021).\n- [30] T. Damour, P. Jaranowski, and G. Schafer, Physical Review D 62 , 044024 (2000).\n- [31] A. Le Tiec, Physical Review D 92 , 084021 (2015).\n- [32] T. J. Rothe and G. Schafer, Journal of Mathematical Physics 51 , 082501 (2010).\n- [33] D. L. Richardson and T. J. Kelly, Celestial Mechanics 43 , 193 (1988).\n- [34] F. Biscani and S. Carloni, Monthly Notices of the Royal Astronomical Society 428 , 2295-2310 (2012).\n- [35] F. Biscani and S. Carloni, Monthly Notices of the Royal Astronomical Society 446 , 3062 (2014).\n- [36] A. Morbidelli, Modern celestial mechanics : aspects of solar system dynamics (Taylor & Francis, London, 2002).\n- [37] C. D. Murray and S. F. Dermott, Solar System Dynamics (Cambridge University Press, Cambridge, 2000).\n- [38] T. Damour and N. Deruelle, C. R. Seances Acad. Sci. (1981).\n- [39] T. Damour, P. Jaranowski, and G. Schafer, Poincar'e invariance in the adm hamiltonian approach to the general relativistic two-body problem (2000), flux-balance equations: derivation of non-conservation laws dissipative sys-  \n- ems., gr-qc/0003051 [gr-qc].\n- [40] T. Damour, P. Jaranowski, and G. Schafer, Physical Review D 77 , 064032 (2008).\n- [41] G. Comp'ere, R. Oliveri, and A. Seraj, Journal of High Energy Physics 2020 , 10.1007/jhep10(2020)116 (2020).\n- [42] C. Konigsdorffer, G. Faye, and G. Schafer, Phys. Rev. D 68 , 044004 (2003).\n- [43] T. Damour and G. Schafer, Nuovo Cimento B Serie 101B , 127 (1988).\n- [44] E. Barausse, E. Racine, and A. Buonanno, Physical Review D 80 , 104025 (2009).\n- [45] L. Blanchet, A. Buonanno, and A. Le Tiec, Physical Review D 87 , 024030 (2013).\n- [46] A. Buonanno and T. Damour, Physical Review D 59 , 084006 (1999).\n- [47] C. R. Galley, Physical Review Letters 110 , 174301 (2013).\n- [48] C. R. Galley, D. Tsang, and L. C. Stein, The principle of stationary nonconservative action for classical mechanics and field theories (2014), 1412.3082 [math-ph].\n- [49] B. Deme and J.-B. Fouvry, Milankovitch equations with spinors (2022), 2212.05929 [physics.class-ph]."}