{'fraction_non_alphanumeric': 0.04608770120768303, 'fraction_numerical': 0.01586940944124998, 'mean_word_length': 4.124305928712162, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Post-Newtonian diagnostic of quasi-equilibrium binary configurations of compact objects

Using equations of motion accurate to the third post-Newtonian (3PN) order (O(v/c)^6 beyond Newtonian gravity), we derive expressions for the total energy E and angular momentum J of the orbits of compact binary systems (black holes or neutron stars) for arbitrary orbital eccentricity. We also incorporate finite-size contributions such as spin-orbit and spin-spin coupling, and rotational and tidal distortions, calculated to the lowest order of approximation, but we exclude the effects of gravitational radiation damping. We describe how these formulae may be used as an accurate diagnostic of the physical content of quasi-equilibrium configurations of compact binary systems of black holes and neutron stars generated using numerical relativity. As an example, we show that quasi-equilibrium configurations of corotating neutron stars recently reported by Miller et al. can be fit by our diagnostic to better than one per cent with a circular orbit and with physically reasonable tidal coefficients.

an issue, because neutron stars are much larger objects, so the numerical simulations necessarily commence at larger separations, where PN methods are presumably more reliable.
Numerical simulations of compact binary inspiral start with a solution of the initial value equations of Einstein's theory; these provide the initial data for the evolution equations (some initial-data models [18] solve in addition one of the six dynamical field equations). The initial state is assumed to consist of two compact objects (neutron stars or black holes) in an initially circular orbit. For stellar-mass systems that have evolved in isolation for eons, gravitational radiation is expected to leave the orbit in an accurately circular state, apart from the adiabatic inspiral induced by the loss of orbital energy; that inspiral is ignored in the initial-data models. (Miller has analysed the consequences of this particular assumption [19]).
The circular-orbit condition is imposed by demanding that dr/dt = 0 initially, where r is a measure of the orbital separation. One way to achieve this is to require that the system have an initial "helical Killing vector" (HKV), which corresponds to a kind of rigid rotation of the binary system. Some initial-data models assume that the objects are co-rotating, a condition which is astrophysically unlikely, albeit computationally advantageous, while others assume that the bodies are irrotational, i.e. non-rotating in an inertial frame. To simplify the problem, an approximation for the spatial metric is generally made; one is the assumption of conformal flatness, an approximation that is known to be invalid in full general relativity. This approximation is usually justified by the neglect of radiation reaction in the initial state. Other approximations, derived from post-Newtonian theory, or from sums of Kerr geometries, have also been used. For black hole binaries, suitable horizon boundary conditions must be imposed, while for neutron star binaries, equations of hydrostationary equilibrium and an equation of state must be provided.
One important product of these initial value solutions is a relationship between the energy E and angular momentum J of the system as measured at infinity, and the orbital frequency Ω. The energy could be the total energy as measured at infinity, consisting of the masses of the two stars plus the orbital energy, or it could be the total energy less the energy of the same two stars in isolation. The latter quantity would be a measure of the orbital binding energy. As all quantities are well-defined and gauge invariant, they are useful variables for making comparisons with PN methods.
We have developed formulae for E(Ω) and J(Ω) using PN methods. Our analytic formulae include point-mass terms through 3PN order, but ignore radiation reaction. They also include rotational energy and spin-orbit and spin-spin terms for the case in which the bodies are rotating. They further include a Newtonian calculation of the effects of tidal and rotational distortions, applicable to stars of arbitrary density distribution, expressed in terms of so-called "apsidal constants" (i.e. we do not restrict attention to homogeneous ellipsoids [20]), and including effects at quadrupole and octupole order. We verify that, for black holes, tidal effects can be ignored, while for neutron-star binaries, they must be included. In contrast to previous work [16,21,22,23], our formulae apply to general eccentric orbits, not just to circular orbits.
In an earlier paper [24], we compared this formula with HKV numerical solutions for corotating binary black holes obtained by Grandclément et al. [18], for the regime where the black holes are separated from the location of the innermost circular orbit by a factor of around two, where PN results might be expected to work well (Gm/rc 2 ∼ 0.1). We found that when we assumed circular PN orbits, our 3PN formulae for E(Ω) and J(Ω) agreed to within 0.5 % with other PN methods, including our own formulae truncated at 2PN order, and 3PN formulae derived using resummation or Padé techniques. However all PN methods consistently and systematically underestimated the binding energy and overestimated the angular momentum, compared to the values derived from the numerical HKV initial-data models, by amounts that were up to 10 times larger than the spread among the PN methods. But when we relaxed the assumption of a circular orbit and demanded only that dr/dt = 0, our PN formula could be made to agree extremely well with the numerical data by assuming that the system being simulated is initially at the apocenter of a slightly eccentric orbit. For values of GmΩ/c 3 ranging from 0.03 to 0.06, corresponding to orbital v/c between 0.3 and 0.4, or orbital separation between 10 and 6 Gm/c 2 , nearly perfect agreement with the binding energy and the angular momentum could be obtained with eccentricities that range from 0.03 to 0.05.
The concordance within fractions of a percent between the various 2PN, 3PN and resummation PN results matches expectation, since (Gm/rc 2 ) 3 ∼ 10 −3 . Presuming that all relevant physical effects have been included, we argued that the PN results in this range of GmΩ/c 3 are robust. We suggested the possibility that the approximations made in most numerical initial-data models could lead to an apparent eccentricity in what was expected to be a quasi-circular orbit. At present, however, the discrepancy between the two approaches can only be considered a hint of possible eccentricity, because the results of [18] did not include quantitative error bars for the variables E(Ω) and J(Ω).
These results motivate us to propose a "post-Newtonian diagnostic", a tool that can be used to extract physical information from numerical simulations, and that may also be an aid to guide some of the assumptions and approximations inherent in numerical initial data computations toward those that lead to the desired physical configuration, such as a true quasi-circular orbit.
In this paper we provide the physical assumptions, mathematical details, and justifications for the approximations that underly this proposed diagnostic tool. We give the detailed foundations for the analysis carried out in [24] for black-hole binary systems, and also extend that work to the case of neutron-star systems by including tidal effects.
As an application of our diagnostic to neutron-star systems, we analyse recent numerical models of quasi-equilibrium orbits of neutron stars by Miller et al. [25]. In contrast to the black-hole case, we find that the orbital energy in the neutron-star initial-data models of [25] can be fit to better than one percent, and importantly, within the error bars provided in [25], using circular orbits with physically reasonable tidal parameters appropriate to the "Γ = 2" equation of state used in that numerical work. The results illustrate the robustness of the PN approximation well into the strongly relativistic regime of compact binaries, especially when augmented with physically movitated finite-size effects. Application of this PN diagnostic to other numerical models will be subject of future papers.
The remainder of this paper provides the details underlying these conclusions. In Sec. II, we solve the post-Newtonian equations of motion calculated to third post-Newtonian (3PN) order, for general eccentric orbits. Neglecting radiation reaction effects, we then express the total conserved orbital energy and angular momentum in terms of a pair of "covariant" orbit elements e (eccentricity) and ζ (related to the semi-latus rectum). In Sec. III, we calculate the effects of finite size in binary systems with bodies whose spin axes are perpendicular to the orbital plane. These include tidal and rotational distortions, spin-orbit terms and spin-spin terms. In Sec. IV, we analyse our diagnostic quantitatively, and apply it to co-rotating, equal-mass binaries of black holes and of neutron stars. Two appendices provide the detailed derivations of the expressions for the tidal and rotational distortion included in our diagnostic: Appendix A uses Newtonian gravity to solve the general problem of the equilibrium configurations of gravitating bodies disturbed by an external force, paralleling the treatment in the classic monographs of Kopal [26,27], and Appendix B specializes the results to linear perturbations caused by rotational and tidal disturbances.

II. ENERGY AND ANGULAR MOMENTUM FOR "POINT" MASSES TO 3PN ORDER
A. Orbits at the turning point in post-Newtonian gravity Since our ultimate focus will be on orbits that are possibly eccentric, but that momentarily haveṙ = 0, it will be useful to review the characteristics of orbits at turning points in Newtonian theory. In Newtonian gravity, the orbit of a pair of point masses may be described by the set of equations where p = a(1 − e 2 ) is the semi-latus rectum (a is the semi-major axis), ω is the angle of pericenter, m = m 1 + m 2 is the total mass, µ = m 1 m 2 /m is the reduced mass, and E and J are the total orbital energy and angular momentum, respectively (henceforth we use units in which G = c = 1). A circular orbit corresponds to e = 0, with r = a = constant, Ω 2 = m/a 3 , E/µ = a 2 Ω 2 /2 − m/a = −(mΩ) 2/3 /2, and J/µ = √ ma = (m/Ω) 1/3 . However, if we demand only that the orbit be at apocenter, so thatṙ = 0 only, we have φ = ω + π, r = p/(1 − e), Ω 2 = (m/p 3 )(1 − e) 4 , so that, in terms of Ω a , the angular velocity at apocenter, To obtain expressions in terms of Ω p , the angular velocity at pericenter, one makes the replacements Ω a → Ω p and e → −e in Eqs. (2.2). However, at higher PN orders, neither the orbital eccentricity e nor the semi-latus rectum p is uniquely or invariantly defined. One definition of eccentricity used by Lincoln and Will [29] in their analysis of orbits at 2.5PN order was that of a Newtonian orbit momentarily tangent to the true orbit (the "osculating" eccentricity); it had the unusual property that it did not tend to zero for a circular PN orbit, but tended toward a constant value of order m/p, while the rate of pericenter advance approached the same rate of rotation as the orbit itself. In this language, the true orbit was a non-circular orbit at perpetual periastron, thereby maintaining a constant separation r. In an effort to avoid this anomaly, other authors [30] adopted a "quasi-Keplerian" parametrization, which defined multiple "eccentricities" to encapsulate different aspects of non-circular orbits at PN order.
In an effort to find a parametrization of non-circular PN orbits that will be useful in comparing with numerical models, we [24] proposed an alternative measure of eccentricity and semi-latus rectum according to: where Ω p is the value of Ω where it passes through a local maximum (pericenter), and Ω a is the value of Ω where it passes through the next local minimum (apocenter). These definitions have the following virtues: (1) they reduce precisely to the normal eccentricity e and semi-latus rectum p in the Newtonian limit, as can be verified from Eqs. (2.1); (2) they are constant in the absence of radiation reaction; (3) they are somewhat more directly connected to measurable quantities, since Ω is the angular velocity as seen from infinity (eg. as measured in the gravitational-wave signal) and one calculates only maximum and minimum values, without concern for the coordinate location in the orbit; and (4) they are straightforward to calculate in a numerical model of orbits without resorting to complicated definitions of "distance" between bodies.
They have the defect that, when radiation reaction is included, they are not local, continuously evolving variables, but rather are some kind of orbit-averaged quantities (for this reason, they may not be as "covariant" as they seem -see Sec. II E below). Nevertheless, when an eccentric orbit decays and circularizes under radiation reaction the definition of e has the virtue that it tends naturally to zero when the orbital frequency turns from ocillatory behavior to monotonically increasing behavior (i.e. the maxima and minima merge).
By virtue of these definitions, ζ has the further property that (2.4) We will derive expressions for orbital energy and angular momentum in terms of these parameters e and ζ; for comparision with numerical models of quasi-equilibrium parametrized in terms of Ω atṙ = 0 (Ω a or Ω p ), one can simply substitute for ζ from Eq. (2.4). In this section we will focus on 3PN expressions for point masses; in the next section, we will incorporate effects due to rotation and finite size.

B. 3PN equations of motion
We use the standard form of the equations of motion, written in a "Newtonian-like" manner. The acceleration of body 1 is given schematically by where x a and m a denote the position and the mass of the body a, r is the separation between the two bodies, n = (x 1 − x 2 )/r is the unit vector from 2 to 1, and v = v 1 − v 2 the relative velocity. The equation for body 2 is obtained by making the replacement 1 ⇋ 2. The notation nP N represents the n th post-Newtonian correction to Newtonian gravity. These equations are valid only for point-like, non-spinning bodies. Post-Newtonian terms nP N include even (integer) and odd (half-odd integer, such as 2.5PN, or 5/2 PN) orders. Even terms are conservative, in the sense that the equations of motion admit conserved quantities such as energy and angular momentum. Odd terms correspond to gravitational radiation reaction, and therefore are not conservative. In particular, they will cause the orbit to shrink, and the eccentricity to decrease.
We convert the two-body problem to an effective one-body problem. For this purpose we choose the origin to be at the center of mass of the system, which is defined by an integral of the motion (a conserved quantity to the 3PN order of approximation to which we will be working). We then change all variables to the relative coordinates x = x 1 − x 2 using relations of the type where η = µ/m = m 1 m 2 /m 2 is the reduced mass parameter (0 < η ≤ 1/4), and δm = m 1 − m 2 . The result is a set of equations of motion in terms of relative coordinates: where A and B represent post-Newtonian terms. To date, the two-body equations of motion have been computed up to and including 3.5PN order. In an appropriate harmonic gauge, writing A = A 1 + A 2 + · · · and B = B 1 + B 2 + · · · , the expressions for A and B read [31]:    At 3PN order, the computation implemented by Blanchet et al. [8,9] produced logarithmic terms, proportional to ln(r/r ′ 1 ) and ln(r/r ′ 2 ), where r ′ 1 and r ′ 2 are constants related to a scale of radius for each body. In obtaining Eqs. (2.8) and (2.9), we removed these logarithms using a 3PN coordinate transformation x µ → x µ + δx µ , with [9]: where y a = |x − x a | denotes the coordinate separation between the considered point and the body a. We note that we have δx µ = 0, except at the location of the two bodies. This ensures that the harmonic condition is still respected in the new gauge to the required order. In addition, the parameter λ, which was initially undetermined in [8,9,31] has now been fixed to be λ = −1987/3080 by different techniques [11,13,32]; that value has been incorporated into all equations.
In the absence of the 2.5PN and 3.5PN terms, these equations of motion admit conserved total energy E and total angular momentum J. Writing E = E 0 + E 1 + E 2 + E 3 and J = J 0 + J 1 + J 2 + J 3 , we have: (2.12d)

C. Solution of the 3PN equations of motion
In order to solve these equations, we shall initially adopt the method of osculating orbital elements, which is welladapted to the perturbed two-body Kepler problem. The osculating orbit elements are defined by the Keplerian orbit that is tangent to the actual trajectory at a particular moment of time. In the Newtonian case, the osculating elements are constants of the motion; in a perturbed Newtonian problem, they change smoothly with time (see [29] for more details about the method of osculating elements applied to the post-Newtonian problem).
From the equations of motion we can easily deduce that the trajectory is planar, which allows us to reduce the number of variables from six to four. If we assume that the plane of the motion is perpendicular toẑ (x, y, z being a standard cartesian coordinate system), our new set of variables (α, β, p, φ) is related to the old set (x, y, v x , v y ) by the definitions (some of which are redundant): Reciprocally, we can deduce the osculating elements from the orbital variables by using the following relations: (2.14) One additional expression will be useful:ṙ We note that the vector (α, β) has as its norm the ordinary Keplerian osculating eccentricity e and as its phase angle the direction ω of the Keplerian osculating periastron, so that we have α = e cos ω and β = e sin ω.
In what follows, we will use the parameter u = m/p rather than p. Note that u is of order ǫ ∼ m/r. In the Newtonian case, u, α and β are constants of the motion; in the post-Newtonian problem, these parameters vary according to the following "Lagrange planetary equations" (so-called from their extensive use in solar-system studies): , we get a set of coupled first-order differential equations in the variables α(φ), β(φ) and u(φ). The planetary equations derived from Eqs. (2.16) are too long to be reproduced here (they can be found through 2.5PN order in [29]). However we can schematically write them in the general form: where Du n , Dα n and Dβ n (n ∈ {1, 2, 5/2, 3, 7/2}) are polynomials in α and β and simple trigonometric functions of φ. We quote, for illustration, the first post-Newtonian expressions for these polynomials: We want to solve these equations iteratively. At zeroth (Newtonian) order u, α and β are constants of the motioñ u,α andβ, and can be related to the initial state of the orbit. Post-Newtonian effects cause them to vary slowly over a post-Newtonian timescale or a radiation-reaction timescale, related to the orbital phase φ by φ/ǫ and φ/ǫ 5/2 , respectively. Superimposed upon this will be variations on an orbital timescale. To take these two effects into account, we use a two-scale approach [33]. We define a variable θ = ǫφ, and we assume that the osculating elements can be written as functions of θ and φ in the generic form u = u(ũ(θ),α(θ),β(θ), φ), with θ and φ now treated as independent variables. We then expand the elements in powers of ǫ: Notice that, by its very nature, u begins at order ǫ. We write the derivative with respect to φ in the form (2.20) We also expand the derivatives with respect to θ in powers of ǫ: Now we have reduced our study to the search for α i , β i , u i on the one hand, which will give the dependence oñ α,β,ũ, and φ, and dα i , dβ i , dũ i on the other hand, which will give differential equations allowing solution for the θ-dependence, or long-term variation of the parameters. Note that this is not the only way to decompose the problem, but is a natural way, given the split into orbital and secular evolution of the variables. We now define the average and the average-free part of a function f (φ) by: where the "independent" variable θ is held fixed. (An equivalent procedure would be to convert all functions of φ into 2π-periodic functions and constants.) We rewrite Eqs. (2.16) with our new variables, and we collect terms of common powers of ǫ. At first order in ǫ we get where the expressions on the right-hand-side are given by Eqs. (2.18), withα replacing α, and so on. Reading off the average parts of Eqs. (2.23), we find dũ 1 = 0, dα 1 = −3ũβ, and dβ 1 = 3ũα. Definingα ≡ẽ cosω andβ ≡ẽ sinω we find, to first PN order that These results express the well-known fact that the orbital eccentricity and semi-latus rectum do not evolve secularly to 1PN order; in fact, this holds true at 2PN and 3PN order; they only evolve secularly as a result of radiation reaction. The angle of pericenterω evolves secularly at 1PN order via the standard advance; there are also 2PN and 3PN contributions, but no radiation-reaction contributions to the advance ofω, through 3.5PN order. Then, integrating the average-free parts of Eqs. (2.23), we obtain, for example, The role of the second AF is to get rid of the constant of integration. The same method yields similar results for β 1 and u 1 . At second order in ǫ, we obtain equations of the form where α 1 , β 1 , u 1 , dα 1 , dβ 1 and dũ 1 are known from the first order solution. For the same reasons as previously we have: Using this procedure systematically up to 3.5PN order, we completely determine α(α,β,ũ, φ), β(α,β,ũ, φ) and u(α,β,ũ, φ), as well as dα dφ (α,β,ũ), dβ dφ (α,β,ũ) and dũ dφ (α,β,ũ). From this and Eqs. (2.13) we can deduce the explicit expressions for x, v, r, etc.
To 3.5PN order, the secular evolution ofũ andẽ is governed by radiation reaction, and is given by the coupled equations (we now set ǫ = 1) We note that the eccentricity decreases as the orbit shrinks. The periastron advance is driven by the conservative part of the equations: (2.29)

D. Energy and angular momentum in terms of new orbit elements
We now wish to convert from the osculating orbit elementsũ andẽ to our alternative quantities defined in Eqs.
we can easily show that the maxima and minima of Ω occur at φ =ω (pericenter) and φ =ω + π, (apocenter) respectively. We then express Ω p and Ω a , and thence our new orbit elements e and ζ as functions ofẽ andũ. To 2PN order, the relationships are given by Notice that a circular orbit corresponds toẽ = e = 0. We invert these relations and substitute the expressions forẽ(e, ζ) andũ(e, ζ) into the solution of the equations of motion. The results for m/r and r 2φ to 3PN order are too long to be reproduced here. However, in order to give an idea of what they look like, we quote them to 2PN order, expressed in terms of our new orbit elements.
The leading term corresponds to the Newtonian solution. Note that ζ, e and ω are now our post-Newtonian orbital elements, and should not be mistaken for the Newtonian u, e and ω introduced in Sec. II A. An alternative method for integrating the post-Newtonian equations of motion was developed by Wagoner and Will [28]. In that method, perturbations of the velocity and angular momentum were defined by the equations where p,ê and ω are constants. Taking a time derivative of both equations, substituting the 3PN equations of motion (ignoring radiation reaction terms), converting to derivatives with respect to φ and integrating, one obtains expressions for the perturbed r 2 dφ/dt and v. One then integrates the identity (d/dφ)r −1 = −(r 2 dφ/dt) −1 (n · v) with respect to φ, setting the constants of integration at each PN order so that the identity d(rn)/dt = v is reproduced.
In terms of the bare orbit elementsê,û = m/p and ω, the orbit equations look different at each PN order from those derived above in terms ofẽ,ũ andω. But when the solution so derived is used to identify Ω a and Ω p and thence to define new orbit elements from Eqs. (2.3), the resulting orbit solution in terms of our new orbit elements is identical to Eqs. (2.32), through 3PN order. The equations describing the evolution with time of our new orbit elements then become: Now the problem is entirely solved. Equations (2.32) pushed to 3PN order, characterize the motion, while Eqs. (2.34) give the pericenter advance and effect of radiation reaction on the orbital elements. We now ignore the effects of radiation reaction, express all the orbital variables r, x, v,ṙ to 3PN order in terms of our new orbit elements and the angle φ, and substitute into the expressions (2.11) and (2.12). As expected E and J are constant (independent of φ) through 3PN order. DefiningẼ = E/µ andJ = |J|/µm, with J = |J|ẑ, we find for a general eccentric orbit: Notice that E Harm is proportional to (1 − e 2 ) through 3PN order, indicating that E Harm = 0 for the limiting unbound orbit e = 1; this is another appropriate feature of our "covariant" eccentricity. The energy and angular momentum are well-defined, physically observable quantities, so one can alternatively express our orbit elements ζ and e as functions ofẼ andJ. Here we give the results to 1PN order, but the calculation can be done to 3PN order: The foregoing results are valid in harmonic gauge. That gauge is characterized by the condition ∂ ν h µν = 0, where h µν = √ −gg µν − η µν , g µν and g are the physical metric and its determinant, and η µν is a background Minkowski metric. In this gauge, Einstein's equations take the form where = η µν ∂ µ ∂ ν is the flat d'Alembertian operator, and the source term τ µν depends both on the matter stressenergy tensor T µν and on non-linear contributions of the gravitational field. The local equation of motion ∇ ν T µν = 0 is equivalent to ∂ ν τ µν = 0, which follows from the harmonic gauge condition. There actually is an infinity of distinct harmonic gauges, and the equations of motion will generally depend on the choice of a particular gauge. We already saw an example of this in the choice of eliminating logarithmic terms from the 3PN contributions [Eq. (2.10) above]. A different approach to the two-body problem, implemented through 3PN order by Damour, Jaranowski and Schäfer [4,5,6], is to compute the Hamiltonian of the system rather than the equations of motion. Unlike other methods, this does not use a harmonic coordinate system, but a so-called ADM (Arnowitt-Deser-Misner), or "Hamiltonian" gauge, or coordinate system. It has been proven to be equivalent to the harmonic formulation [10].
The Hamiltonian has been computed up to 3PN order; because it is a Hamiltionian approach, it explicitly suppresses the 2.5PN and 3.5PN contributions of radiation reaction. We quote it here only to 1PN order (see [6] for a complete expression): (2.38) We convert this two-body problem into an effective one-body problem by using the simple relation p = p 1 = −p 2 , valid in the center-of-mass frame. Thus we get a new expression for H ADM (x, p). From Hamilton's equations: we iteratively extract the equations of motion and write them in the same form as equation (2.7), but with different A and B. Substituting the expression of p as a function of v and x into the Hamiltonian, we obtain the total conserved energy E ADM . Similarly, we get J ADM by calculating x × p. For both the equations of motion and the expressions for energy and angular momentum, the harmonic and ADM-Hamiltonian terms coincide at 1PN order, but they differ at 2PN and 3PN orders. We apply the method described in section II C to find solutions to the ADM equations of motion and expressions for E ADM and J ADM in terms of the osculating orbit elements. In this case,ẽ andũ are strictly constant because radiation reaction is not present in the Hamiltonian approach. We then find expressions for our new orbit elements e and ζ in terms ofẽ andũ and write E ADM and J ADM in terms of these elements. The results are: (2.40b) We observe two features of the harmonic and the ADM versions of these expressions: (i) the "circular" parts (e = 0) of the formulae coincide. In that case the angular velocity Ω = Ω a = Ω p is the same as that observed from infinity for both harmonic and ADM coordinates; (ii) the expressions also coincide for η → 0, i.e. in the test-mass limit. As mentioned before, the differences between the formulae only occur at 2PN and 3PN orders. It is actually possible to relate the coordinate positions and velocities in the two gauges. In particular, the relation betweenφ ADM andφ Harm , r Harm , etc. allows us to find a relation between (e ADM , ζ ADM ) and (e Harm , ζ Harm ), and thus account for the differences in the coefficients of E and J. We found that a transformation of the typė This is reflected in Eq. (2.41) where the difference between the two angular velocities is of 2PN order. Furthermore, for the small eccentricity orbits that we wish to consider, the corrections are proportional to e, and are thus further suppressed. Thus we argue that our definitions of e and ζ are "almost" covariant.

A. Estimates for compact binaries
In reality, the bodies in our binary system cannot be treated as purely point masses. They may be rotating, and thus subject to a number of effects, including rotational kinetic energy, rotational flattening, and spin-orbit and spin-spin interactions. Furthermore, there will be tidal deformations. These effects will not only make direct contributions to the energy and angular momentum of the system, they may also modify the equations of motion, and thereby modify the expressions for our alternative eccentricity and semi-latus rectum. However because they depend on the size of the bodies, which, for neutron stars and black holes, are of order m, we expect these effects to be "effectively" of high PN order, even if they are Newtonian in origin, such as tidal effects. To see this, we estimate each finite-size effect in turn and compare it with the Newtonian orbital energy E N ∼ m 2 /r. We assume that the rotational angular velocity ω of each body ranges from zero to the orbital angular velocity, given by Ω ∼ (m/r 3 ) 1/2 , and we let the radius of each body be of the form R a ∼ qm a , where q ∼ 1 for black holes (in harmonic coordinates), and q ∼ 5 for neutron stars.
• Rotational kinetic energy: E Rot ∼ Iω 2 /2 ≤ mR 2 (m/r 3 ) ∼ E N q 2 (m/r) 2 . This is effectively 2PN order. There will be PN corrections to the kinetic energy, given by E Rot−PN ∼ E Rot (Rω) 2 ∼ mR 4 ω 4 ∼ E N q 4 (m/r) 5 . These are effectively 5PN order, but, because of the q 4 dependence, could be important for neutron stars.
• Rotational flattening: E Flat ∼ δIω 2 /2, where δ is a measure of the deformation of the body, given by the ratio of rotational to gravitational energy, δ ∼ (Iω 2 )/(m 2 /R), so that E Flat ∼ ω 4 R 5 ≤ E N q 5 (m/r) 5 . There is an equivalent contribution of rotational flattening to the gravitational internal energy. These are effectively 5PN order, but because of the q 5 dependence, could be important for neutron stars.
• Tidal deformations: E Tidal ∼ (δ ′ m) 2 /R, where δ ′ is the ratio of gravitational energy due to the tidal force of the companion to the internal gravitational energy of the body, δ ′ ∼ (mR 2 /r 3 )/(m/R) ∼ (R/r) 3 . Thus E Tidal ∼ m 2 (R/r) 6 /R ∼ E N q 5 (m/r) 5 . There is also a contribution from the rotational kinetic energy of the tidal bulge, given by E KE−bulge ∼ δ ′ Iω 2 ∼ mR 5 ω 2 /r 3 ≤ E N q 5 (m/r) 5 . These are effectively 5PN order, but could be significant for neutron stars.
A parallel heirarchy of finite-size effects applies to the total angular momentum of the system. The largest effect in principle is that due to the rotational kinetic energy of the bodies and thus requires some care. For black holes, we can apply the general formulas for mass and angular momentum of isolated Kerr black holes, in terms of the irreducible mass and angular velocity. For neutron stars, no such general formula exists, so it may be necessary to rely upon numerical results for energy and angular momentum of isolated rotating neutron star models in order to take accurate account of this effect. On the other hand, it does not directly affect the equations of motion.
Because the remaining effects are effectively of 3PN order and higher, our strategy will be to evaluate them analytically to the lowest non-trivial order. For tidal and rotational flattening terms, this will mean using Newtonian theory. For spin-orbit and spin-spin terms, we will use the well-known 1PN formulae. We will ignore any coupling among these effects, or between these effects and the point-mass PN effects described in the previous section. Accordingly, we will calculate the separate contribution of each effect to the energy and angular momentum and simply add them all up.

B. Newtonian Tidal and Rotational Effects
In Appendix A we derived the general form of the equations of motion and the conserved energy and angular momentum for a binary system of tidally and rotationally deformed bodies, and in Appendix B we specialized to linear perturbations and multipole indices l = 2 and l = 3. We now specialize further to systems more relevant to the initial configurations in numerical relativity which we wish to study, namely binary systems in which the spin axes of both stars are perpendicular to the orbital plane. where the three perturbing terms correspond repectively to the effects of rotational distortions, quadrupole tidal distortions (l = 2) and octupole tidal distortions (l = 3), with the coefficients given by 3 m 1 /m 2 ) .

(3.2)
For each body, R a denotes its radius, k denote the "apsidal constants" for angular harmonics l = 2 and l = 3, respectively, andω a denotes the body's angular velocity at a chosen point in the orbit (see Appendix B for details). Apsidal constants are dimensionless coefficients that depend on the degree of central condensation of the star, and that determine the size of distortion of a given angular degree l produced by a given external perturbation. Note that A < R 5 /m 2 r 3 ∼ q 5 (m/r) 3 , so that, despite appearances, this term, like the purely tidal term from B, is effectively 5PN order. The energy and angular momentum that are conserved by virtue of the full fluid equations of motion are given by where, for each body, I a denotes the moment of inertia, W a denotes the self-gravitational energy of the undistorted configuration, andr denotes the orbital separation at the point at which the star's angular velocity isω. The chosen point in our case will be the pericenter or apocenter. In Eq. (3.3), the split among the intrinsic energy and spins of the bodies E Self and S, the constant distortion terms E Distort and J Distort , and the orbital terms is clear. The angular momentum components are all referred to the axis perpendicular to the orbital plane.
We now repeat the method of subsections II C -II D to obtain the general solution to the equations of motion to first order in the tidal and rotational perturbations. We then obtain our new orbit elements e and ζ in terms of the bare elementsẽ andũ; for example, e is given by

C. Spin-orbit and Spin-Spin effects
Spin-orbit and spin-spin interactions produce corrections in the equations of motion that are formally of 1PN order. For systems with the spins perpendicular to the orbital plane they are given by where S = S 1 + S 2 and ∆ = m(S 2 /m 2 − S 1 /m 1 ). The individual spins are constants of the motion when they are both aligned perpendicular to the orbital plane. The conserved energy and total angular momentum are given by where L N = µx × v, and Eq. (3.9b) denotes the component perpendicular to the orbital plane (for the complete equations of motion, see, for example, [35,36]). We define the dimensionless quantities where D ∼ F ∼ (R/r) 2 ∼ q 2 (m/r) 2 , and G ∼ (R/r) 4 ∼ q 4 (m/r) 4 , making the spin-orbit and spin-spin terms effectively 3PN and 5PN order respectively [34]. With these definitions, the equation of motion takes the form of Eq.
(2.7), with (3.12) Again we solve the equations of motion using the method of subsections II C -II D and define our new orbit elements. In this case, for example, the eccentricity is given by In terms of our new elements, the pericenter advance is given by (3.14) while e and ζ undergo no secular changes. When expressed in terms of our new orbit elements, the spin-orbit and spin-spin contributions to the total energy and angular momentum have the form In deriving the "point-mass" equations of motion, the underlying assumption was that the masses that enter the equations are the total mass of each body, comprised of baryonic mass, gravitational binding energy and rotational kinetic energy, if any. Thus, each m a should be written m a = m B a − W a + E Rot a . In many numerical approaches, sequences of models are constructed in which the total (or ADM) mass of each corresponding non-rotating star is held fixed along the sequence. Thus, for making comparisons with such sequences, we should replace each m a in Eqs. (2.35) or (2.40) with m 0 a + I a ω 2 a /2 (or, in the case of black holes, with a suitable formula in terms of the irreducible mass and ω a ). But because I a ω 2 a ∼ q 2 (m/r) 3 , the main contribution, at effectively 3PN order, comes from making this replacement in the Newtonian expressions. Expressing E N and J N in terms of Ω as E N = − 1 2 ηm(1 − e 2 )(mΩ a ) 2/3 /(1 − e) 4/3 , and J N = ηm 2 (1 − e) 2/3 /(mΩ a ) 1/3 , and making the above replacement, we find the corrections to the Newtonian energy and angular momentum where all masses now are those of the equivalent non-rotating body. For neutron stars, this would be that of the same baryonic mass; for black holes, it would be that of the same irreducible mass.

A. Estimates of effects
We now have all the ingredients to formulate a post-Newtonian diagnostic for quasi-equilibrium configurations of compact bodies. The ingredients are the various contributions to the total energy and angular momentum of the system in terms of the "covariant" orbit elements e and ζ, together with the relationships connecting the value of ζ with the orbital angular velocity at a turning point of the orbit, namely ζ = (mΩ p ) 2/3 /(1 + e) 4/3 or ζ = (mΩ a ) 2/3 /(1 − e) 4/3 , corresponding to pericenter and apocenter, respectively. The ingredients are: • Point-mass orbital contributions through 3PN order. Eqs. (2.35) or (2.40). It is straightforward to show that, because the harmonic and ADM versions differ by 2PN terms proportional to ηe 2 and higher, the differences between the two versions are negligible for all cases of interest. Henceforth we will adopt the harmonic version of Eqs. (2.35).
• Self Terms. Eqs. (3.7a) and (3.7b). We add a suitably defined total "rest" mass for the bodies to the definition of E Self . Because the rotational kinetic energy and the spin angular momentum are effectively of 2PN order, they will have to be treated with some care.
In order to assess the applicability of this diagnostic, we first study the sizes of various effects for systems of interest. In general we will consider systems of solar-mass scale neutron stars or black holes, in circular or small-eccentricity orbits, in the vicinity of the onset of an unstable plunge and merger. This corresponds to ζ ∼ m/r < 1/5 for black holes, or to ζ < m/(2R) ∼ 1/q for neutron stars. For q between 4 and 6, the two ranges are comparable. Both correspond to mΩ < 0.1. We will generally choose a range 0.01 < mΩ < 0.1. First we look at the relative contributions of point-mass PN corrections. Figure 1 shows the contribution, relative to the Newtonian orbital energy, of the 1PN, 2PN and 3PN terms in the energy, for η = 0 and η = 1/4, as a function of mΩ. Results for the angular momentum are similar. While the 1PN terms are essentially insensitive to η, and the 2PN terms are only 15 % smaller for equal masses than for the test-mass limit, the 3PN terms are suppressed for equal masses by more than a factor of 10 compared to the test-mass limit. As Blanchet [16,17] has argued, this suggests that the 3PN approximation may be quite accurate for comparable-mass systems, without the need for sophisticated resummation techniques. At the largest angular velocity considered, 3PN terms contribute less than one per cent of the total binding energy and angular momentum of the orbit.
Next we consider the effects of tidal and rotational distortions. We consider systems of identical bodies (m 1 = m 2 ) which are corotating (ω 1 =ω 2 = Ω). For neutron stars, we adopt the maximum values of the apsidal constants (k 2 = 3/4 and k 3 = 3/8, see Appendix B 3), and choose two representative values of q = R a /m a for neutron star models with reasonable equations of state, namely q = 4 and q = 6. The results, plotted as a fraction of the Newtonian orbital terms, are shown in Figures 2 and 3, along with the PN contributions for comparison. As expected, tidal effects are very sensitive to the stellar radii. For q = 4, the l = 2 tidal terms become comparable to the 2PN and 1PN terms only around mΩ ∼ 0.09, while the l = 3 terms are an order of magnitude smaller. For q = 6, the l = 2 tidal terms exceed the 1PN terms already by mΩ ∼ 0.05, while the l = 3 terms are small, approaching the 2PN terms only at the largest allowed mΩ ∼ 0.07, corresponding to the point at which these larger stars are touching. For irrotational stars (ω 1 =ω 2 = 0), the tidal effects are very similar.
These curves illustrate that tidal effects need to be taken into account carefully in an accurate diagnostic for neutron star binaries, but are not so large that they invalidate our approximation scheme. Their modest size also supports our use of Newtonian theory to calculate them. They only become problematical for the largest neutron stars near the very endpoint of their inspiral. It should also be pointed out that, in making these estimates, we have adopted the largest values of the apsidal constants, corresponding to uniform-density stars. While neutron stars are not as centrally condensed as, say, non-degenerate stars, they are also not uniform density, so the k l may well be smaller than their maximum values. For example, for a Newtonian polytrope, p = kρ Γ , with Γ = 2, k 2 = 0.26, so the q = 6 tidal terms in Fig. 3 are reduced by a factor of three, bringing them to a level at or below the 1PN terms over the whole range of mΩ. On the other hand, very little, if anything, is known about the values of k l for general relativistic neutron stars over a range of equations of state. This is a subject that we are currently investigating. Figure 4 shows the effects of tides for co-rotating black-hole binaries. There we choose q = 1 (R = m in harmonic coordinates), k 2 = 3/4 and k 3 = 3/8 (for slowly rotating black holes, k 2 from rotational distortions happens to be precisely 3/4; see, eg. [37]). We see, not surprisingly, that tidal effects are utterly negligible over the entire range of mΩ.
Finally, we examine spin effects. Again we consider identical, co-rotating bodies. For neutron stars, we assume that S a = I a Ω, with the moment of inertia given by that for a uniform density body, I a = (2/5)m a R 2 a = (2/5)q 2 m 3 a . The results are shown also in Figures 2 and 3. For q = 4, spin-orbit effects are small but significant, just below the 2PN terms, while spin-spin effects are negligible. For q = 6, spin-orbit terms exceed 2PN terms by mΩ ∼ 0.04 and become comparable to the 1PN terms by the maximum angular velocity, while spin-spin terms barely exceed the 3PN effects.
For black holes, we use the fact that S a = 4m 3 a Ω. Figure 4 shows that the spin-orbit terms lie between the 2PN and 3PN contributions and thus must be included, while spin-spin terms are negligible (though larger than the tidal terms).

B. Corotating, identical black holes
For black hole binaries, we ignore tidal and spin-spin effects. We set m 1 = m 2 , η = 1/4, andω 1 =ω 2 = Ω. We exploit the fact that there exist exact formulae for the energy and spin of isolated Kerr black holes in terms of the irreducible mass, irr ω/[1 − 4(M irr ω) 2 ] 1/2 . The total energy and angular momentum of the system are then given by where where m irr is the total irreducible mass of the system, given by (m irr ) 1 + (m irr ) 2 . In Eqs. (4.2a) and (4.2b), we have expanded the Kerr formulae for M and S in powers of m irr Ω, assumed to be small compared to unity, keeping as many higher-order terms as needed to reach a precision comparable to our 3PN formulae. To obtain E Tot and J Tot at a turning point as functions of Ω, we substitute ζ = (m irr Ω a ) 2/3 /(1 − e) 4/3 or ζ = (m irr Ω p ) 2/3 /(1 + e) 4/3 for apocenter or pericenter, respectively (in calculating E N,Corr and J N,Corr , we have already changed the dependence in ζ from the total mass of the rotating bodies to the total irreducible mass of the non-rotating counterparts). These are the formulas used in [24] to compare with the numerical HKV quasi-equilibrium solutions of Grandclément et al. [18]. When E Tot and J Tot are scaled by m irr and m 2 irr respectively, there remains only one free parameter, the eccentricity of the orbit, and we found [24] that a substantially better fit to the numerical data was obtained for non-zero values of e, of the order of 0.03, with the system at apocenter, than for e = 0. We suggested that such apparent eccentricity could be a result of the inevitable approximations (such as the conformally flat approximation) and numerical errors in such initial-data models, but, in the absence of detailed estimates of the sizes of those errors, it was difficult to draw firm conclusions. On the other hand, those engaged in numerical models of black hole binaries could use our diagnostic as a guide to know when, say, a suitable circular orbit has been achieved, or whether further numerical experiments with different grid sizes or larger computational domains are necessary to reach the desired physically meaningful state.

C. Corotating, identical neutron stars
For neutron stars, we must include tidal effects. We set m 1 = m 2 , η = 1/4, andω 1 =ω 2 = Ω; we let the apsidal constants and radius factors be common for both stars, given by k 2 , k 3 , and q, respectively, and express all quantities in terms of the total mass m 0 = (m 0 ) 1 + (m 0 ) 2 of two non-rotating stars with the same equation of state. We also define for each star the coefficient α a = I a /m a R 2 a , and also assume it to be common for both stars. The result is where the 3PN point-mass expressions E Harm and J Harm are given in Eqs. (4.2c) and (4.2d), and where (1 + 7e 2 + 7e 4 + e 6 )q 7 k 3 ζ 13/2 , (4.4h) We also need to fix the coefficients q and α. From data provided by Miller (private communication), the radius of each isolated non-rotating star in isotropic coordinates is given by R I = 6.77M ADM ,while the baryonic moment of inertia, calculated using isotropic coordinates, is given by I 0 = 9.412M 3 0 . We work in harmonic coordinates, but since R H = R I (1 + M 2 ADM /4R 2 I ), the difference between the two coordinates is only of order 1/2 %, so we read off q = 6.77. The ADM moment of inertia can be identified as Also shown is a model with the same Γ = 2 apsidal constants, but with a non-zero eccentricity e = 0.02 and with the system at apocenter. This marginally fits the numerical data within the error bars, but consistently gives lower (more negative) energies. We conclude that these quasi-equilibrium neutron-star configurations are fit to better than one per cent by our PN diagnostic with a circular orbit, and with physically reasonable tidal terms.
In future work we plan to compare this diagnostic with results of other numerical models of quasi-equilibrium black hole and neutron star binaries. Our 3PN equations of motion, together with tidal and spin terms, augmented by radiation reaction terms, can also be used to develop a "dynamical" diagnostic, to compare with numerical simulations of evolutions from the quasi-equilibrium initial data [25,38]. quantities (such as the external tidal potential, seen either from the global inertial frame, or from the rotating frame of a given body) is sufficiently long that each body can be assumed to be in hydrostatic equilibrium. In other words, we will ignore dynamical tides [39]. This is a reasonable assumption as long as we are focusing on quasi-equilibrium initial data. Consider one of the bodies in the binary system. From the equation of hydrostatic equilibrium, ∇p = ρ∇Ψ, where p, ρ and Ψ are the pressure, density and total gravitational potential, respectively, we conclude that ∇ρ × ∇Ψ = 0, and thus that surfaces of constant ρ and Ψ coincide. We label surfaces of constant ρ by the radial parameter a, and let the equation of those surfaces have the form where Y lm (Ω) are spherical harmonics corresponding to the directionΩ, and where the dimensionless distortion functions f lm have the property f * lm = (−1) m f l,−m . On general grounds we expect f lm ∼ (R/r) l+1 ∼ q l+1 (m/r) l+1 for tidal effects, and, for l = 2, f 2m ∼ ω 2 /ρ ∼ (R/r) 3 ∼ q 3 (m/r) 3 from rotational effects. The effect of these distortions on the external potential of a body is of order f lm (R/r) l ∼ q 2l+1 (m/r) 2l+1 . For l = 2, this means effectively 5PN order; l = 3 effects would be effectively 7PN order, and so on. However, for neutron stars, with q ∼ 4 and m/r ∼ 0.1, an l = 2 distortion effect becomes numerically comparable to a 2PN term, while l = 3 is comparable to a 3PN term. For black holes, with q < 1, the effects are much smaller. Thus, in the end, we will keep only l = 2 and l = 3 distortion terms. Also, non linear corrections to f lm would be of order (R/r) l+1 ∼ q l+1 (m/r) l+1 smaller than the dominant linear effects, and thus, effectively of 8PN order for l = 2 (for neutron stars, these non-linear corrections would be numerically smaller than 3PN). The exception to this is in the internal gravitational energy of each body, where a quadratic contribution yields (m 2 /R)f lm 2 ∼ (m 2 /r)(R/r) 2l+1 , which is comparable to the other effectively 5PN contributions for l = 2. We begin, however, with a general analysis, keeping lm arbitrary, and working to second order in the small quantities f lm . Later (Appendix B) we will specialize to l = 2 and l = 3 linear perturbations. To second order, it is straightforward to show that, for any n, where and C lm αβ;γδ is defined in terms of Clebsch-Gordan coefficients, C lm αβ;γδ = (2α + 1)(2γ + 1) 4π(2l + 1) Note that the various angular momentum quantum numbers are connected by the constraints l = α + γ , α + γ − 2 , . . . , |α − γ|, and m = β + δ; the C lm αβ;γδ are symmetric under (αβ) ⇋ (γδ). Also note that X 00 = (16π) −1/2 αβ f αβ f * αβ . We expand the gravitational potential U of the body and the disturbing potential V in the form where the subscript > (<) corresponds to the larger (smaller) of r and r ′ . The disturbing potential consists of a part, with disturbing coefficients d lm , that corresponds to a potential with ∇ 2 V = 0, such as the gravitational potential from another body, or the Laplacian-free part of a centrifugal potential, plus the spherical part of a centrifugal potential, with coefficientd. We now substitute Eqs. (A1) and (A2) into (A5), convert all expressions from r to a, and demand that, for l = 0, the external gravitational potential of our body have the form U = m/r (i.e. the perturbation does not change the mass of the body), and that, for l = 0, the total potential U + V be constant at a given a. The first can be satisfied if f 00 + 2X 00 = 0, while the second holds if, for l = 0, where m(a) = a 0 4πρ(a)a 2 da , and A denotes the surface of the body. The left-hand-side of Eq. (A6) is first order in f lm , while the right-hand-side is second-order. Dividing the first-order terms by a l , differentiating with respect to a and multiplying by a 2l+2 , we obtain the first-order result where prime denotes differentiation with respect to a. Substituting this and the first-order solution of Eq. (A6) back into the right-hand-side of Eq. (A6), it is straightforward to show that the term involving ′ αβ;γδ reduces to m(a)(2X lm − aX ′ lm )/4πa, to second order. The basic equation for the distortion functions f lm can then be written in the form of an integral equation where R lm contains all contributions quadratic in small quantities: R lm (a) = 2l + 1 4π am(a)X ′ lm + 2(2l + 1)da 2 f lm Combining Eq. (A9) with various derivatives of it, one obtains the following useful equations, evaluated at the surface a = A of the star: where m = A 0 4πρa 2 da, and Another combination of first and second derivatives of Eq. (A9) yields a second-order differential equation for f lm , sometimes called Clairaut's equation: For a given density distribution ρ(a), this equation can be solved, subject to the boundary conditions that f lm be regular at a = 0, and that, at the surface, f lm satisfy Eq. (A11b).

Energy and angular momentum of the system
Given a solution for the distortion functions f lm (a), we can calculate all the quantities needed for the equations of motion and the energy and angular momentum of the orbit. The external potential of our body, for example, is given by where ′ denotes summation for l = 0, and where we define the "apsidal constant" k lm by The total gravitational energy of the system is given by The self-energy W 11 of body 1 can be written Substituting Eqs. (A1) and (A2) we find no contribution linear in f lm . To second order, we obtain Since the second term is already second order, we can integrate by parts and use the first-order versions of Eqs. (A11) and (A13) to obtain the alternative form where we define the self-gravitational binding energy W of the undistorted configuration by In W 12 , we substitute the external potential of body 2 evaluated inside body 1, to obtain where y 2 = x − x 2 , k (2) lm is the apsidal constant of body 2, and d (2) lm is the coefficient of the disturbing potential acting on body 2. Since the interaction energy is smaller than the self energy by a factor of R/r, we only need to keep terms linear in the deformations f lm or the disturbing coefficients d lm ; consequently we carry out a multipole expansion of 1/y 2 in the first term in Eq. (A21), then convert from r to a using Eq. (A1), but we evaluate the second term at the center of mass of body 1 and do the lowest-order spherical integral. Effectively, we are ignoring multipole-multipole coupling between the bodies, which can be shown to lead to effects of order (A/r) 2l+1 (A/r) 2n+1 ∼ (m/r) 10 , or 10PN order for l = n = 2. The result is where now r = |x 1 − x 2 | and n = (x 1 − x 2 )/r. Combining Eqs. (A19) and (A22), the final result is where, under the interchange, n → −n. The kinetic energy of the system is given by T = ρv 2 d 3 x. Splitting the velocity of an element of fluid into center-of-mass, rotational, and random parts, and noting that m Y * lm (n)Y lm (n ′ ) = (2l + 1)!! 4πl! n <L> n ′ <L> = 2l + 1 4π P l (n · n ′ ) , wheren <L> denotes an STF product of l unit vectors (a capitalized superscript denotes a multi-index), the product n <L> n ′ <L> denotes contraction on all indices, and P l is a Legendre polynomial, we may write whereλ 1 = ω 1 /ω 1 . Converting from r to a using Eq. (A1), recalling that f 00 = −2X 00 , and noting that ω 2 is already of first order in disturbing quantities, we obtain, to second order in small quantities, where I 1 = (2/3) A1 0 4πρa 4 da. The angular momentum of the system is given by J i = ǫ ijk ρx j v k . Using the same split of the velocities, we obtain, to the analogous order of precision, where we define the symmetric trace free (STF) tensor with the following properties with the coefficient for body 2 obtained by interchange. Working to linearized order, we can factor out the azimuthal m dependence by defining f T lm = f T l Y * lm (n), then the Clairaut equation and the outer boundary condition for body 1 take the form Note that the apsidal constant depends only on f T l , and is thus independent of m, so with a corresponding expression for body 2, where .
Note that the overall scale of f T l has no effect on the apsidal constant, to linear order. For a uniformly rotating body, the disturbing potential at a point x ′ is the centrifugal potential Thus we read off the rotational coefficient for body 1 with the coefficient for body 2 obtained by interchange (the spherical coefficientd only contributes at second order in small quantities). Similarly defining f R 2m = f R 2 Y * 2m (λ), we find that f R 2 also satisfies Eq. (B3a) for l = 2, but with the boundary condition The rotational apsidal constant is also independent of m, and, since the overall scale of f R 2 (A 1 ) is irrelevant, it is equal to the l = 2 tidal apsidal constant: k 2 . This equality will not hold when non-linear corrections are included (it will also not hold in general relativity, when frame-dragging and other relativistic effects are included).