In Section 7.3 we have given the equations of motion of non-spinning compact binary systems in the frame of the center-of-mass for general orbits at the 3PN and even 3.5PN order. We shall now investigate (in this section and the next one) the explicit solution to those equations. In particular, let us discuss the general “doubly-periodic” structure of the post-Newtonian solution, closely following Refs. [142, 143*, 149*].
The 3PN equations of motion admit, when neglecting the radiation reaction terms at 2.5PN order, ten first integrals of the motion corresponding to the conservation of energy, angular momentum, linear momentum, and center of mass position. When restricted to the frame of the center of mass, the equations admit four first integrals associated with the energy and the angular momentum vector , given in harmonic coordinates at 3PN order by Eqs. (4.8) – (4.9) of Ref. .
The motion takes place in the plane orthogonal to . Denoting by the binary’s orbital separation in that plane, and by the relative velocity, we find that and are functions of , , and . We adopt polar coordinates in the orbital plane, and express and the norm , thanks to , as some explicit functions of , and . The latter functions can be inverted by means of a straightforward post-Newtonian iteration to give and in terms of and the constants of motion and . Hence,
where and denote certain polynomials in , the degree of which depends on the post-Newtonian approximation in question; for instance it is seventh degree for both and at 3PN order [312*]. The various coefficients of the powers of are themselves polynomials in and , and also, of course, depend on the total mass and symmetric mass ratio . In the case of bounded elliptic-like motion, one can prove  that the function admits two real roots, say and such that , which admit some non-zero finite Newtonian limits when , and represent respectively the radii of the orbit’s periastron (p) and apastron (a). The other roots are complex and tend to zero when .
Let us consider a given binary’s orbital configuration, fully specified by some values of the integrals of motion and corresponding to quasi-elliptic motion.70 The binary’s orbital period, or time of return to the periastron, is obtained by integrating the radial motion asper orbital revolution, per period is given by . As tends to one in the limit (as is easily checked from the usual Newtonian solution), it is often convenient to pose , which will then entirely describe the relativistic precession.
Let us then define the mean anomaly and the mean motion by
Here denotes the instant of passage to the periastron. For a given value of the mean anomaly , the orbital separation is obtained by inversion of the integral equation336*) by writing it in the form
Here denotes a certain function of the mean anomaly which is periodic in with period , hence periodic in time with period . According to Eq. (336*) this function is given in terms of the mean anomaly by338) exhibits clearly the nature of the compact binary motion, which may be called doubly periodic in that the mean anomaly is periodic with period , and the periastron advance is periodic with period . Notice however that, though standard, the term “doubly periodic” is misleading since the motion in physical space is not periodic in general. The radial motion is periodic with period while the angular motion is periodic [modulo ] with period where . Only when the two periods are commensurable, i.e., when where , is the motion periodic in physical space (with period ).
The quasi-Keplerian (QK) representation of the motion of compact binaries is an elegant formulation of the solution of the 1PN equations of motion parametrized by the eccentric anomaly (entering a specific generalization of Kepler’s equation) and depending on various orbital elements, such as three types of eccentricities. It was introduced by Damour & Deruelle [149*, 150] to study the problem of binary pulsar timing data including relativistic corrections at the 1PN order, where the relativistic periastron precession complicates the simpler Keplerian solution.
In the QK representation the radial motion is given in standard parametric form as334). The “angular” eccentricity enters the equation for the angular motion at 1PN order which is written in the form 71 per orbital revolution defined by Eq. (333*); it may be written as where is the angle of return to the periastron.
Crucial to the formalism are the explicit expressions for the orbital elements , , , , and in terms of the conserved energy and angular momentum of the orbit. For convenience we introduce two dimensionless parameters directly linked to and by
where is the reduced mass with the total mass (recall that for bound orbits) and we have used the intermediate standard notation . The equations to follow will then appear as expansions in powers of the small post-Newtonian parameter ,72 with coefficients depending on and the dimensionless reduced mass ratio ; notice that the parameter is at Newtonian order, . We have [149*]
The dependence of such relations on the coordinate system in use will be discussed later. Notice the interesting point that there is no dependence of the mean motion and the radial semi-major axis on the angular momentum up to the 1PN order; such dependence will start only at 2PN order, see e.g., Eq. (347a).
The above QK representation of the compact binary motion at 1PN order has been generalized at the 2PN order in Refs. [170*, 379*, 420*], and at the 3PN order by Memmesheimer, Gopakumar & Schäfer [312*]. The construction of a generalized QK representation at 3PN order exploits the fact that the radial equation given by Eq. (331a) is a polynomial in (of seventh degree at 3PN order). However, this is true only in coordinate systems avoiding the appearance of terms with the logarithm ; the presence of logarithms in the standard harmonic (SH) coordinates at the 3PN order will obstruct the construction of the QK parametrization. Therefore Ref. [312*] obtained it in the ADM coordinate system and also in the modified harmonic (MH) coordinates, obtained by applying the gauge transformation given in Eq. (204*) on the SH coordinates. The equations of motion in the center-of-mass frame in MH coordinates have been given in Eqs. (222); see also Ref. [9*] for details about the transformation between SH and MH coordinates.
in which the true anomaly is still given by Eq. (343*). The new orbital elements , , , , , , and parametrize the 2PN and 3PN relativistic corrections.73 All the orbital elements are now to be related, similarly to Eqs. (345), to the constants and with 3PN accuracy in a given coordinate system. Let us make clear that in different coordinate systems such as MH and ADM coordinates, the QK representation takes exactly the same form as given by Eqs. (340*) and (346). But, the relations linking the various orbital elements , , , , , , to and or and , are different, with the notable exceptions of and .
Indeed, an important point related to the use of gauge invariant variables in the elliptical orbit case is that the functional forms of the mean motion and periastron advance in terms of the gauge invariant variables and are identical in different coordinate systems like the MH and ADM coordinates [170*]. Their explicit expressions at 3PN order read
Because of their gauge invariant meaning, it is natural to use and as two independent gauge-invariant variables in the general orbit case. Actually, instead of working with the mean motion it is often preferable to use the orbital frequency which has been defined for general quasi-elliptic orbits in Eq. (337*). Moreover we can pose
Besides the very useful gauge-invariant quantities , and , the other orbital elements , , , , , , , , , , , parametrizing Eqs. (340*) and (346) are not gauge invariant; their expressions in terms of and depend on the coordinate system in use. We refer to Refs. [312*, 9*] for the full expressions of all the orbital elements at 3PN order in both MH and ADM coordinate systems. Here, for future use, we only give the expression of the time eccentricity (squared) in MH coordinates:
Again, with our notation (344), this appears as a post-Newtonian expansion in the small parameter with fixed “Newtonian” parameter .
In the case of a circular orbit, the angular momentum variable, say , is related to the constant of energy by the 3PN gauge-invariant expansion
This permits to reduce various quantities to circular orbits, for instance, the periastron advance is found to be well defined in the limiting case of a circular orbit, and is given at 3PN order in terms of the PN parameter (230*) [or (348*)] by
See Ref.  for a comparison between the PN prediction for the periastron advance of circular orbits and numerical calculations based on self-force theory in the small mass ratio limit.
The gravitational wave energy and angular momentum fluxes from a system of two point masses in elliptic motion was first computed by Peters & Mathews [340*, 339*] at Newtonian level. The 1PN and 1.5PN corrections to the fluxes were provided in Refs. [416, 86*, 267*, 87*, 366*] and used to study the associated secular evolution of orbital elements under gravitational radiation reaction using the QK representation of the binary’s orbit at 1PN order . These results were extended to 2PN order in Refs. [224*, 225] for the instantaneous terms (leaving aside the tails) using the generalized QK representation [170, 379, 420]; the energy flux and waveform were in agreement with those of Ref.  obtained using a different method. Arun et al. [10*, 9*, 12*] have fully generalized the results at 3PN order, including all tails and related hereditary contributions, by computing the averaged energy and angular momentum fluxes for quasi-elliptical orbits using the QK representation at 3PN order , and deriving the secular evolution of the orbital elements under 3PN gravitational radiation reaction.74
The secular evolution of orbital elements under gravitational radiation reaction is in principle only the starting point for constructing templates for eccentric binary orbits. To go beyond the secular evolution one needs to include in the evolution of the orbital elements, besides the averaged contributions in the fluxes, the terms rapidly oscillating at the orbital period. An analytic approach, based on an improved method of variation of constants, has been discussed in Ref. [153*] for dealing with this issue at the leading 2.5PN radiation reaction order.
The generalized QK representation of the motion discussed in Section 10.2 plays a crucial role in the procedure of averaging the energy and angular momentum fluxes and over one orbit.75 Actually the averaging procedure applies to the “instantaneous” parts of the fluxes, while the “hereditary” parts are treated separately for technical reasons [10*, 9*, 12*]. Following the decomposition (308*) we pose where the hereditary part of the energy flux is composed of tails and tail-of-tails. For the angular momentum flux one needs also to include a contribution from the memory effect [12*]. We thus have to compute for the instantaneous part
Thanks to the QK representation, we can express , which is initially a function of the natural variables , and , as a function of the varying eccentric anomaly , and depending on two constants: The frequency-related parameter defined by (348*), and the “time” eccentricity given by (350). To do so one must select a particular coordinate system – the MH coordinates for instance. The choice of rather than (say) is a matter of convenience; since appears in the Kepler-like equation (346a) at leading order, it will directly be dealt with when averaging over one orbit. We note that in the expression of the energy flux at the 3PN order there are some logarithmic terms of the type even in MH coordinates. Indeed, as we have seen in Section 7.3, the MH coordinates permit the removal of the logarithms in the equations of motion, where is the UV scale associated with Hadamard’s self-field regularization [see Eq. (221*)]; however there are still some logarithms which involve the IR constant entering the definition of the multipole moments for general sources, see Theorem 6 where the finite part contains the regularization factor (42*). As a result we find that the general structure of (and similarly for , the norm of the angular momentum flux) consists of a finite sum of terms of the type351*). The coefficients , and are straightforwardly computed using the QK parametrization as functions of and the time eccentricity . The ’s correspond to 2.5PN radiation-reaction terms and will play no role, while the ’s correspond to the logarithmic terms arising at the 3PN order. For convenience the dependence on the constant has been included into the coefficients ’s. To compute the average we dispose of the following integration formulas ()76
In the right-hand sides of Eqs. (353b) and (353c) we have to differentiate times with respect to the intermediate variable before applying . The equation (353c), necessary for dealing with the logarithmic terms, contains the not so trivial function353a) we see that there will be no radiation-reaction terms at 2.5PN order in the final result; the 2.5PN contribution is proportional to and vanishes after averaging since it involves only odd functions of .
Finally, after implementing all the above integrations, the averaged instantaneous energy flux in MH coordinates at the 3PN order is obtained in the form 348*). The various instantaneous post-Newtonian pieces depend on the symmetric mass ratio and the time eccentricity in MH coordinates as
The Newtonian coefficient is nothing but the Peters & Mathews  enhancement function of eccentricity that enters in the orbital gravitational radiation decay of the binary pulsar; see Eq. (11*). For ease of presentation we did not add a label on to indicate that it is the time eccentricity in MH coordinates; such MH-coordinates is given by Eq. (350). Recall that on the contrary is gauge invariant, so no such label is required on it.
The last term in the 3PN coefficient is proportional to some logarithm which directly arises from the integration formula (353c). Inside the logarithm we have posed42*). Only after adding the hereditary contribution to the 3PN energy flux can we check the required cancellation of the constant . The hereditary part is made of tails and tails-of-tails, and is of the form [10*]
where , , , and are certain “enhancement” functions of the eccentricity.
Among them the four functions , , and appearing in Eqs. (359) do not admit analytic closed-form expressions. They have been obtained in Refs. [10*] (extending Ref. [87*]) in the form of infinite series made out of quadratic products of Bessel functions. Numerical plots of these four enhancement factors as functions of eccentricity have been provided in Ref. [10*]; we give in Figure 3* the graph of the function which enters the dominant 1.5PN tail term in Eq. (358*).
Furthermore their leading correction term in the limit of small eccentricity can be obtained analytically as 
On the other hand the function in factor of the logarithm in the 3PN piece does admit some closed analytic form:
The latter analytical result is very important for checking that the arbitrary constant disappears from the final result. Indeed we immediately verify from comparing the last term in Eq. (356d) with Eq. (359c) that cancels out from the sum of the instantaneous and hereditary contributions in the 3PN energy flux. This fact was already observed for the circular orbit case in Ref. ; see also the discussions around Eqs. (93*) – (94*) and at the end of Section 4.2.
Finally we can check that the correct circular orbit limit, which is given by Eq. (314), is recovered from the sum . The next correction of order when can be deduced from Eqs. (360) – (361*) in analytic form; having the flux in analytic form may be useful for studying the gravitational waves from binary black hole systems with moderately high eccentricities, such as those formed in globular clusters .
Previously the averaged energy flux was represented using – the gauge invariant variable (348*) – and the time eccentricity which however is gauge dependent. Of course it is possible to provide a fully gauge invariant formulation of the energy flux. The most natural choice is to express the result in terms of the conserved energy and angular momentum , or, rather, in terms of the pair of rescaled variables (, ) defined by Eqs. (344). To this end it suffices to replace by its MH-coordinate expression (350) and to use Eq. (349) to re-express in terms of and . However, there are other possible choices for a couple of gauge invariant quantities. As we have seen the mean motion and the periastron precession are separately gauge invariant so we may define the pair of variables (, ), where is given by (348*) and we pose
As we are interested in the phasing of binaries moving in quasi-eccentric orbits in the adiabatic approximation, we require the orbital averages not only of the energy flux but also of the angular momentum flux . Since the quasi-Keplerian orbit is planar, we only need to average the magnitude of the angular momentum flux. The complete computation thus becomes a generalisation of the previous computation of the averaged energy flux requiring similar steps (see Ref. [12*]): The angular momentum flux is split into instantaneous and hereditary contributions; the instantaneous part is averaged using the QK representation in either MH or ADM coordinates; the hereditary part is evaluated separately and defined by means of several types of enhancement functions of the time eccentricity ; finally these are obtained numerically as well as analytically to next-to-leading order . At this stage we dispose of both the averaged energy and angular momentum fluxes and .
The procedure to compute the secular evolution of the orbital elements under gravitational radiation-reaction is straightforward. Differentiating the orbital elements with respect to time, and using the heuristic balance equations, we equate the decreases of energy and angular momentum to the corresponding averaged fluxes and at 3PN order . This extends earlier analyses at previous orders: Newtonian  as we have reviewed in Section 1.2; 1PN [86, 267]; 1.5PN [87, 366] and 2PN [224, 153*]. Let us take the example of the mean motion . From Eq. (347a) together with the definitions (344) we know the function at 3PN order, where and are the orbit’s constant energy and angular momentum. Thus,
have already been used at Newtonian order in Eqs. (9). Although heuristically assumed at 3PN order, they have been proved through 1.5PN order in Section 5.4. With the averaged fluxes known through 3PN order, we obtain the 3PN averaged evolution equation assecular evolution under gravitational radiation reaction for eccentric orbits. The complete evolution includes also, superimposed on the averaged adiabatic evolution, some fast but smaller post-adiabatic oscillations at the orbital time scale [153, 279].