Secular dynamics of binaries in stellar clusters-III. Doubly-averaged dynamics in the presence of general relativistic precession

Secular evolution of binaries driven by an external (tidal) potential is a classic astrophysical problem. Tidal perturbations can arise due to an external point mass, as in the Lidov-Kozai (LK) theory of hierarchical triples, or due to an extended stellar system (e.g. galaxy or globular cluster) in which the binary resides. For many applications, general-relativistic (GR) apsidal precession is important, and has been accounted for in some LK calculations. Here we generalise and extend these studies by exploring in detail the effect of GR precession on (quadrupole-level) tidal evolution of binaries orbiting in arbitrary axisymmetric potentials (which includes LK theory as a special case). We study the (doubly-averaged) orbital dynamics for arbitrary strengths of GR and binary initial conditions and uncover entirely new phase space morphologies with important implications for the binary orbital evolution. We also explore how GR precession affects secular evolution of binary orbital elements when the binary reaches high eccentricity (e → 1) and delineate several different dynamical regimes. Our results are applicable to a variety of astrophysical systems. In particular, they can be used to understand the high-eccentricity behaviour of (cluster) tide-driven compact object mergers — i.e. LIGO/Virgo gravitational wave sources — for which GR effects are crucial.


INTRODUCTION
The problem of the relative motion of two bound point masses has formed the basis of celestial mechanics since it was first successfully tackled mathematically by Newton in 1687. In 1915 Einstein updated the solution, showing that the lowest order correction to Newton's elliptical orbit (in the small parameter G(m1 + m2)/ac 2 , with mi the constituent masses and a the binary semimajor axis) was simply an extra prograde apsidal precession at a ratė ωGR =ω where e is the orbital eccentricity, andωGR|e=0 is the GR precession rate for a circular orbit. Einstein's solution is now known as the first post-Newtonian (1PN) approximation to the two-body problem. A century on, the LIGO/Virgo Collaboration has detected, and continues to detect, dozens of merging compact object (black hole or neutron star) binaries (The LIGO Scientific Collaboration et al. 2018, 2020. These discoveries certainly warrant an astrophysical explanation, which is complicated by the fact that the timescale for an isolated compact object binary to merge via gravitational wave (GW) emission is often much longer than the age of the Universe. For instance, an isolated circular black hole binary with m1 = m2 = 30M⊙ will only merge within 10 10 yr if its initial ⋆ E-mail: ch783@cam.ac.uk † John N. Bahcall Fellow at the Institute for Advanced Study semimajor axis a is 0.2AU. Thus nature must have a way of forcing these relativistic binaries to such small separations.
One way to achieve this outcome is by driving an initially wide binary to a very high eccentricity, e → 1, so that, for a given semimajor axis, a binary's pericentre distance p ≡ a(1 − e) is greatly diminished. In this case the repeated close approaches of the binary components allow significant energy and angular momentum to be dissipated in bursts of GWs, efficiently shrinking the binary orbit and accelerating the merger. Thus, in recent years much effort has gone into searching for mechanisms by which binaries might achieve very high eccentricity. One broad category of proposed mechanisms consists of secular eccentricity excitation of binaries by some perturbing tidal 1 potential. This could be the tidal potential due to a tertiary point mass (e.g. a star) that is gravitationally bound to the binary, in which case the dynamics are described by the Lidov-Kozai (LK) theory (Lidov 1962;Kozai 1962), or simply the mean field potential of the star cluster or galaxy in which the binary resides (Heisler & Tremaine 1986;Brasser et al. 2006;Hamilton & Rafikov 2019a,b;Bub & Petrovich 2019).
However, when investigating such merger channels it is almost always necessary to account for the effect of (1PN) GR precession of the binary's pericentre angle. This is because GW emission primar-ily occurs during close pericentre passages when the binary is highly eccentric, and this is precisely the regime in which GR precession is most important (equation (1)). For similar reasons it is often necessary to include GR precession (as well as other short-range precession effects such as rotational or tidal bulges, see e.g. Liu et al. 2015;Muñoz et al. 2016) in studies of LK secular evolution, which rely on tidal dissipation (inside one or both binary components) to shrink the binary orbit (Fabrycky & Tremaine 2007;Antonini et al. 2016). Tidal dissipation is strongest when e → 1, meaning that GR precession is important as well.
GR precession is routinely accounted for in LK population synthesis calculations (e.g. Antonini et al. 2014;Rodriguez et al. 2015;Liu et al. 2015;Liu & Lai 2018;Hamers et al. 2018;Samsing et al. 2019). Its effect is understood as increasing a binary's prograde apsidal precession rate as e → 1, preventing the perturber from coherently torquing the binary and effectively stopping the reduction of the binary's angular momentum. A result of this so-called 'relativistic quenching' effect is a reduction in the maximum eccentricity that the binary can reach, even if the initial inclination between binary and perturber orbits is favourable (Fabrycky & Tremaine 2007). Some authors have derived approximations to this maximum eccentricity in the limit where GR precession can be treated as a small perturbation to the LK evolution (Miller & Hamilton 2002;Blaes et al. 2002;Wen 2003;Veras & Ford 2010;Liu et al. 2015;Anderson et al. 2017;Grishin et al. 2018). Also, Iwasa & Seto (2016) looked at the modification of the phase space portrait of the LK problem in the presence of GR precession, although their study was far from exhaustive. However, so far nobody has studied carefully the impact of the GR precession for binaries perturbed by general tidal potentials (Brasser et al. 2006;Hamilton & Rafikov 2019c;Bub & Petrovich 2019) where similar considerations should apply.
The main purpose of this paper is to systematically explore the effect of GR precession on the underlying phase space dynamics and eccentricity evolution of a tidally perturbed binary. We will focus exclusively upon the 'doubly-averaged' (hereafter DA) 2 dynamics of binaries perturbed by quadrupole-order tidal potentials. We will also make the test-particle approximation, i.e. assume that the binary's outer orbital motion relative to its perturber contains much more angular momentum than its internal Keplerian orbit. The quadrupolar and test-particle approximations are very good ones for the applications we have in mind (such as compact object binaries perturbed by globular cluster tides), but they can be relaxed, see §5.3. The DA approximation will be relaxed in an upcoming paper (Hamilton & Rafikov, in prep.). The present study complements our investigation of DA cluster tide-driven dynamics of binaries begun in Hamilton & Rafikov (2019a,b), hereafter 'Paper I' and 'Paper II' respectively. The DA theory developed in Papers I and II includes the test-particle quadrupole LK problem as a limiting case, but is more general and dynamically richer, particularly when GR precession is included, as we show here.
In §2 we write down the doubly-averaged perturbing Hamiltonian and establish the notation that we will use for the rest of the paper. In particular we introduce the key parameter ǫGR which measures the strength of GR precession relative to tidal torques. In §3 we explore the phase space behaviour as ǫGR is varied; the quantitative results that we quote in this section are derived in Appendix A. In §4 we investigate very high eccentricity behaviour in the presence of weak or moderate GR precession. In particular we explore how finite ǫGR modifies both the maximum eccentricity reached and the timescale of high eccentricity episodes. In §5 we discuss our results in the light of the existing literature, and comment on the limitations of our study. We summarise in §6. Lastly, in Appendix C we provide for the first time an explicit, analytical solution to the DA equations of motion for all orbital elements in the high eccentricity limit. We also check the accuracy of this solution against direct numerical integration of the DA equations of motion.

DYNAMICAL FRAMEWORK
We consider a binary orbiting in an arbitrary time-independent, axisymmetric external potential Φ. For the remainder of the paper we will refer to this as the 'cluster' potential, though it can in reality be due to an axisymmetric galaxy, point mass, or whatever. We refer to the binary's barycentric motion around the cluster -which is assumed to follow the trajectory of a test particle in the potential Φ, and need not be circular -as the 'outer' orbit. The binary's 'inner' orbit is described by the usual orbital elements: semimajor axis a, eccentricity e, inclination i, argument of pericentre ω, longitude of the ascending node Ω and mean anomaly M . Inclination is measured relative to the plane perpendicular to the cluster's symmetry (Z) axis or, in the case of a spherical potential, relative to the outer orbital plane. The angle to the line of nodes is measured relative to an arbitrary fixed axis in this plane.
Here A is a constant with units of (frequency) 2 which measures the strength of the tidal torque and sets the timescale of secular evolution. It is completely determined by stipulating the form of the cluster potential Φ and the outer orbit of the binary; to order of magnitude, A −1/2 is comparable to the period of the binary's outer orbit. For reference, we provide the value of A for the LK problem, i.e. when Φ is the Keplerian potential (Paper I, Appendix B): where M is the tertiary perturber's mass and ag, eg are respectively the semimajor axis and eccentricity of the binary's outer orbit relative to the tertiary perturber. Next, H * 1 and H * GR are the dimensionless Hamiltonians accounting for quadrupole-order cluster tides and GR pericentre precession, respectively: H * 1 = (2 + 3e 2 )(1 − 3Γ cos 2 i) − 15Γe 2 sin 2 i cos 2ω, The crucial quantity Γ in (4) is a dimensionless parameter which is fully determined (like A) by stipulating Φ and the outer orbitsee §2.1 for discussion. The relative strength of GR precession is measured in equation (5)  (7) Here nK = G(m1 + m2)/a 3 , T b = 2π/nK, v ∼ G(m1 + m2)/a are the Keplerian mean motion, period, and typical orbital speed of the inner orbit of the binary, respectively, while tsec ∼ (T b A) −1 is the timescale of secular eccentricity oscillations (Paper II, equations (33)-(34)). In the numerical estimate (7) we have assumed that the binary is orbiting a spherical cluster with scale radius b and total mass M. Typical values of the dimensionless parameter A * ≡ A/(GM/b 3 ) are mapped out in §6 of Paper I.
As in Paper I we introduce Delaunay variables (actions) L = G(m1 + m2)a, J = L √ 1 − e 2 , and Jz = J cos i, their corresponding angles being M , ω and Ω respectively. Since (2) is independent of the mean anomaly M , the action L is conserved, so we can choose to work with the following dimensionless variables: Obviously j is just the dimensionless angular momentum. The definitions (8) imply that e and j must obey to be physically meaningful for a given Θ. With this notation established, we can rewrite the dimensionless Hamiltonians (4)-(5) in the form Since both Hamiltonians are independent of Ω, the dimensionless quantity Θ is an integral of motion. The total dimensionless Hamiltonian H * = H * 1 + H * GR can be taken as the other integral of motion. The equations of motion fully describing the evolution of the dimensionless variables ω, j are Since ω, j are decoupled from Ω, the evolution of the nodal angle Ω can be explored separately using the equation of motion Obviously, the equation of motion for Jz is trivial, dJz/dt = −∂H/∂Ω = 0.
Given that H * (ω, j) is a constant we can use equations (10)-(11) to eliminate ω from equation (13). Following a derivation analogous to that of equation (30) in Paper II, and without making any approximations, we find where with Note that the definitions of j 2 ± , j 2 0 , Σ, and D are equivalent to those given in Paper II (equations (18), (17), (19) and (15) respectively) except that we have replaced H * 1 in equation (19) by we have used the value of the Hamiltonian that includes GR precession. Therefore in the limit ǫGR = 0, equation (15) reduces to equation (30) of Paper II. Note also that j 2 ± , j 2 0 are not necessarily positive. We will use equation (15) extensively when we study high-eccentricity behaviour in §4.
Some shorthand notation will be necessary as we proceed. In particular, several different values of e and j will come with distinct subscripts. We provide a summary of our notation in Table 1.

A note on Γ
The dimensionless quantity Γ is crucial for our investigation because it determines the morphology of the phase space, and therefore sets fundamental constraints on the allowed dynamical behaviour (Paper II). Its value is computed by time-averaging the outer orbit in the potential Φ (Paper I). Typically it has to be computed numerically but in special cases (e.g. spherically symmetric potentials) it can be calculated (semi-)analytically. For example, the test-particle quadrupole LK problem (e.g. Antognini 2015; Naoz 2016) corresponds exactly to the limit Γ = 1, while the problem of binaries orbiting in the midplane of a thin disk (e.g. Heisler & Tremaine 1986) corresponds to Γ = 1/3. Binaries orbiting inside the potential generated by a homogenous sphere (similar to the inner regions of a globular cluster) effectively have Γ → 0. For a binary on a circular outer orbit of radius R in a Plummer potential with a scale radius b, we get Γ = (1 + 4ζ 2 ) −1 , where ζ ≡ b/R (this special case was considered by Brasser et al. 2006; compare their disturbing function (A.5) with our equation (4)).
In general, Γ can take any value. Moreover, we showed in Paper II (for ǫGR = 0) that at critical Γ values of ±1/5, 0, bifurcations occur in dynamics, meaning that we need to explore separately four (8) e f , j f e, j values of fixed points at ω = ±π/2 when ǫ GR = 0.
(34) emax, j min Maximum e, minimum j values. e lim Upper limit on possible values of eccentricity, e lim = √ 1 − Θ.
(16), (17) regimes: However, one can show that for sensible spherical potentials only 0 < Γ ≤ 1 is possible (Paper I, Appendix D). If one lets the binary population in a cluster simply trace the underlying stellar density profile then it turns out that cored clusters (those with a flat density profile ρ(r) → const. as r → 0, like the Plummer profile) have a significant fraction of their binaries in the 0 < Γ ≤ 1/5 regime (21), and the rest in the Γ > 1/5 regime (20). On the contrary, cusped clusters (those with ρ(r) ∝ r −p as r → 0 for some p > 0, like the Hernquist profile) host a much larger fraction of their binaries in the regime Γ > 1/5 and relatively few in the regime 0 < Γ ≤ 1/5. This has strong implications for the dynamical evolution of binaries in these various types of cluster ( §3).
To get negative Γ values typically requires a highly inclined outer orbit in a sufficiently non-spherical potential (Paper I). Since our applications are mostly concerned with spherical or near-spherical potentials such as those of globular clusters, for which negative Γ values are very rare, we concentrate on the regimes (20)-(21) in the main body of the paper. Discussion of the regimes (22)-(23) can be found in Appendix D .

PHASE SPACE BEHAVIOUR
To gain a qualitative understanding of the dynamics driven by the Hamiltonian equations (12)-(13), one can fix the values of Γ, Θ, ǫGR and then plot (ω, e) phase space trajectories, i.e. contours of constant H * in the (ω, e) plane. We call such a plot a 'phase portrait'.
For Γ > 0 -i.e. in the important regimes (20)-(21) -fixed points always exist in the phase portrait provided Θ is small enough (which can be achieved e.g. by starting with sufficiently large inclination). The precise requirement for that is (Paper II, §2.2) In §4 we will be interested exclusively in situations where some fraction of binaries can reach eccentricities very close to unity (i.e. 1−e ≪ 1). Given that Θ = (1−e 2 ) cos 2 i is conserved, a necessary condition for this is that Θ ≪ 1. For Γ > 0, fixed points always exist for sufficiently small Θ when ǫGR = 0. Neither this nor the existence of fixed points exclusively at ω = π/2 are generally true for ǫGR = 0, as we will see below.
Fixed points are crucial because they may force an initially low-e binary to very high maximum eccentricity emax. Indeed, a key result of Paper II (again with ǫGR = 0) was that for Γ > 1/5, whenever fixed points exist, e f provides a lower bound on emax. Equation (24) implies that e f is close to unity whenever Θ ≪ 1; high eccentricity excitation is then ubiquitous. On the other hand, for 0 < Γ ≤ 1/5 the fixed points no longer provide a lower bound on circulating trajectories' emax and so high-e excitation is much rarer. One consequence of this result is that cored clusters, which have a significant fraction of binaries in the 0 < Γ ≤ 1/5 regime, produce few tidally-driven compact object mergers compared to cusped clusters -see Hamilton & Rafikov (2019c).
In the rest of this section we explore how the phase space behaviour uncovered in Paper II (i.e. for ǫGR = 0) is modified in the case of finite ǫGR, which we do separately for Γ > 1/5 ( §3.1) and for 0 < Γ ≤ 1/5 ( §3.2). We describe some properties of the fixed points in §3.3, and then show how to calculate the maximum eccentricity of a given binary in §3.4. Details of the mathematical results that we quote throughout § §3.1-3.4 are given in Appendix A. Note that the Γ ≤ 0 regimes are treated in Appendix D.
3.1 Phase space behaviour in the case Γ > 1/5 Figure 1 shows phase portraits for the case Γ = 0.5 > 1/5. In the top (bottom) row we set Θ = 0.1 (0.5). From left to right we vary ǫGR taking ǫGR = 0, 2, 5, 10, 30. In each panel a black horizontal dashed line shows the limiting possible eccentricity e lim = √ 1 − Θ (equation (9)). Contours are spaced linearly from the minimum (blue) to maximum (red) value of H * indicated by the colour bar at the top of each panel. Since the linearly sampled contours In the top (bottom) row we fix Θ = 0.1 (Θ = 0.5). We increase ǫ GR from left to right, using the values ǫ GR = 0, 2, 5, 10, 30 indicated in each panel. Contours are spaced linearly from the minimum (blue) to maximum (red) value of H * -see the colour bar at the top of each panel. We have also added by hand dashed contours passing through (ω, e) = (0, 0.01) and (ω, e) = (±π/2, 0.1) in each panel. Dashed black horizontal lines indicate the limiting possible eccentricity e lim = √ 1 − Θ, while fixed points are shown with grey crosses should they exist. Black separatrices illustrate the boundary between families of librating and circulating phase space trajectories. become too widely separated at low eccentricity, to illustrate the low-e behaviour we have added dashed contours passing through (ω, e) = (0, 0.01) and (ω, e) = (±π/2, 0.1) in each panel. Just like in Paper II, trajectories are split into librating and circulating families. We plot the separatrices between these families with solid black lines. Fixed points are denoted with grey crosses.
In panels (a) and (f) we encounter the usual ǫGR = 0 behaviour familiar from the LK problem: (I) there are fixed points at ω = ±π/2, each of which is surrounded by a region of librating orbits, (II) these librating islands are connected to e = 0 line, (III) the family of circulating orbits runs 'over the top' of the librating regions 4 , and (IV) all phase space trajectories reach maximum eccentricity at ω = ±π/2.
Inspecting the other panels, we see that the effect of increasing ǫGR from zero is simply to push the fixed points at ω = ±π/2 to lower eccentricity. As a result, large amplitude eccentricity oscillations along a given secular trajectory are noticeably quenched as ǫGR is increased, and the region of librating orbits is diminished in both area and vertical extent. Eventually the eccentricity of the fixed points reaches zero and so they vanish altogether, leaving only circulating orbits (panels (e), (i) and (j)).
The phase space evolution for non-zero ǫGR exhibited in Figure  1 is characteristic of all systems with Γ > 1/5, including the LK case Γ = 1, which has already been discussed to some degree by Iwasa & Seto (2016) -see §5.2. Of course, the precise characteristics, such as the eccentricity of the fixed points, the critical ǫGR for fixed points to vanish, etc., do depend on the value of Γ, as we detail in §3.3.
3.2 Phase space behaviour in the case 0 < Γ ≤ 1/5 For 0 < Γ ≤ 1/5 (a regime typical of binaries orbiting the inner regions of a cored cluster), a similar but slightly more complex picture emerges. In Figure 2 we show phase portraits similar to Figure 1 except that we now take Γ = 0.1 < 1/5, and pick some new values of ǫGR to better demonstrate the modified phase space behaviour. The strength of GR still increases from left to right. As in Figure 1 we have added in dashed contours that take the values H * (ω = 0, e = 0.01) and H * (ω = ±π/2, e = 0.1).
Starting with the non-GR case ǫGR = 0, we immediately notice a qualitative difference between the phase space morphologies for 0 < Γ ≤ 1/5 (Figures 2a,f) and Γ > 1/5 (Figures 1a,f), discussed at length in Paper II. Although there are again fixed points at ω = ±π/2, the librating islands that surround them are now connected to e = e lim (and not to e = 0, like in the Γ > 1/5 case). As a result, circulating orbits run 'underneath' librating orbits (rather than 'over the top' as for Γ > 1/5) and the maximum eccentricity of circulating orbits is found at ω = 0 (rather than at ω = ±π/2). Crucially, unlike for Γ > 1/5, a binary that starts at low eccentricity does not necessarily reach a high eccentricity even if there are fixed points located near e = 1. This fact is responsible for the dearth of cluster-tide driven mergers in cored clusters, which have many binaries with 0 < Γ ≤ 1/5 (Hamilton & Rafikov 2019c). As we increase ǫGR from zero, the fixed points again get pushed to lower eccentricity (panels (b) and (g)). However, the effect of this for Γ = 0.1 is to initially increase, rather than decrease, the fraction of the phase space area that is encompassed by the librating islands. Additionally, as the fixed points get pushed to lower eccentricity, a new family of high-eccentricity circulating orbits emerges once ǫGR exceeds a threshold value which we determine in §3.3.2. These phase space trajectories run 'over the top' of the fixed points and have their eccentricity maxima at ω = ±π/2 (panels (b), (c) and (h)). Physically the emergence of a new circulating family is expected. This is because at high binary eccentricity (near e = e lim ) GR precession starts to dominate the dynamics, causing the binary's pericentre angle ω to precess rapidly, which prevents the tidal torque from significantly affecting binary eccentricity. The qualitative change from ǫGR = 0 is reflected in the fact that the librating island is now truly an island, disconnected from both e = 0 and e = e lim . This is different from the case Γ > 1/5, in which the lower portion of the librating regions always stretches down to e = 0 until ǫGR becomes so large that fixed points cease to exist.
Simultaneously with the new high-e family of orbits, two saddle points (i.e. fixed points that are not local extrema of H * ) emerge at ω = 0, π. Passing through them are separatrices that isolate the distinct phase space orbital families in panels (b), (c) and (h). This is an entirely new phase space feature that is not found in LK theory, as it is only possible for Γ ≤ 1/5 and only when GR is present, as we show in §3.3.2. The 'two-eyed' phase space structure of panels (b), (c) and (h) has therefore not been uncovered before. Note that for a system exhibiting this structure, a circulating trajectory 'above' the saddle point can have the same H * value as a circulating trajectory 'below' the saddle point. In other words a single value of the Hamiltonian can correspond to two entirely different phase space trajectories. This can be seen in Figures 2c,h where the dashed contours circulating above the librating islands appear because they have the same H * values as the manually added dashed low-e contours passing through (ω, e) = (±π/2, 0.1) and (0, 0.01).
As ǫGR is increased further, the eccentricity of the saddle points diminishes, similar to the fixed points at ω = π/2. It is interesting to note that these various types of fixed points move at different 'speeds' down the phase portrait as ǫGR grows. In particular, panels (d) and (i) of Figure 2 demonstrate that for 0 < Γ ≤ 1/5 there is a range of ǫGR values where the saddle point at ω = 0 has gone below e = 0 and so no longer exists, but the ω = ±π/2 fixed points still do exist.
Even these remaining fixed points get pushed to (and past) e = 0 as ǫGR is increased ever further, leaving the entire phase space filled with circulating trajectories that have their eccentricity maxima at ω = ±π/2, just as for Γ > 1/5 (Figure 2e,j). The amplitude of eccentricity oscillations decreases correspondingly until cluster tides are completely negligible and only GR apsidal precession remains.

Fixed points
We now proceed to understand mathematically the nature of the various fixed points that we found in the phase portraits. By setting dj/dt = 0 in equation (13), we see that all possible non-trivial fixed points 5 are located on (i) the lines ω = ±π/2, as in Paper II, and/or (ii) the lines ω = 0, ±π, consistent with Figures 1 and 2. Finding the j values of the fixed points requires plugging these ω values into dω/dt = 0, given by equation (12), and solving the resulting algebraic equation for j. We do this next for each of the fixed point.

Fixed points at ω = ±π/2
In §A1, we show how to calculate the j value of the fixed points at ω = ±π/2, which we call j f,π/2 , for arbitrary ǫGR and for any Γ > 0, i.e. for both types of phase portraits shown in Figures 1, 2.
The values of j f,π/2 are found as solutions to the quartic polynomial (A1), and we illustrate their behaviour in Figure 3 for several values of Γ and Θ. One can see that j f,π/2 always increases with ǫGR (see (A4)), explaining why in Figures 1,2 the fixed points at ω = ±π/2 always get pushed to lower e as ǫGR is gradually increased from zero. While the explicit expressions for j f,π/2 are too complicated to be shown here, we can gain important insights by considering two limiting cases, namely when ǫGR is much smaller/larger than a particular critical value: (In the top and bottom rows of Figure 1, ǫ π/2 takes values around 4.9 and 16.3 respectively). In the limit ǫGR ≪ ǫ π/2 , which we will call 'very weak GR' regime, the ǫGR term in (A1) is small and we find to lowest order in ǫGR/ǫ π/2 In the opposite limit ǫGR ≫ ǫ π/2 the right hand side in (A1) becomes small and we find to lowest order in ǫ π/2 /ǫGR j f,π/2 ≈ ǫGR 6(1 + 5Γ) Figure 3 shows that these asymptotic solutions match the actual j f,π/2 behaviour in the appropriate limits very well.
In §A1, we also show that for fixed points at (ω, j) = (±π/2, j f,π/2 ) to exist for given Γ > 0, the quantities Θ and ǫGR must obey the inequalities and Here Θ2 is the smallest positive real solution to equation (A5), while and we have defined a quantity that will appear repeatedly throughout this paper.
In Figure 4 we plot Θmax (equation (30)) as a function of Γ for various values of ǫGR (c.f. Figure 1 of Paper II). It is easy to check that the combinations of Γ, Θ and ǫGR that give rise to ω = ±π/2 fixed points in Figure 1 do obey the inequalities (29)-(30). Of course, in the limit ǫGR → 0 equations (30)-(32) reduce to the non-GR constraint (25). Finally we note that for sufficiently small Θ, the conditions (29), (30) reduce simply to the requirement that ǫGR < 2ǫstrong (for Θ ≪ 1).
In other words, if (33) is not satisfied then there are no fixed points even for initially orthogonal inner and outer orbits (i0 = 90 • ). This reflects what we see in Figures 1d,e 2d,e and 6k,l.
3.3.2 Fixed points at ω = 0, π Fixed points at ω = 0, π are unique to the 0 < Γ ≤ 1/5 regime (for Γ > 0). They are always saddle points, and we explore their properties mathematically in §A2. From now on, for brevity we will simply refer to them as being located at ω = 0 rather than ω = 0, ±π, because phase space locations separated in ω by multiples of π are equivalent (see equation (10)).
As we demonstrate in §A2, these saddle points are always located Note that j f,0 is independent of Θ, which can be seen in Figure 2c,h. Also, the constraint (9) implies that fixed points exist at (ω, j) = (0, j f,0 ) if and only if Obviously ǫGR must be finite for the inequality (35) to hold even for very small Θ -hence fixed points at ω = 0 do not exist for ǫGR = 0, which is why they were not found in Paper II and do not exist in Figures 1a,f or 2a,f. In addition we learn from (35) that there are never any fixed points at ω = 0 for Γ > 1/5 regardless of ǫGR, which explains the phase space structure in Figure 1. In particular this implies that for the LK problem (Γ = 1) the only possible fixed point locations are the standard ones at ω = ±π/2, regardless of the value of ǫGR. For positive Γ, fixed points at ω = 0 can be realised only for Γ ≤ 1/5 and we see from (34)-(35) that when ǫGR exceeds the threshold value 6(1 − 5Γ)Θ 3/2 (corresponding to ǫGR = 0.095 and 1.06 in the top and bottom rows of Figure 2, respectively), a fixed point appears at the limiting eccentricity e lim = √ 1 − Θ. Increasing ǫGR at fixed Γ always acts to increase j f,0 , i.e. to decrease the eccentricity e f,0 ≡ (1−j 2 f,0 ) 1/2 of this particular fixed point. As we increase ǫGR beyond the threshold value ǫGR = 6(1 − 5Γ) (which is independent of Θ and corresponds to ǫGR = 3 in Figure 2), the saddle point moves down the page, leaving only the fixed points at ω = ±π/2.
The lower limit here is exact, while the upper limit is correct to zeroth order in Θ. Within this range the qualitative behaviour resembles the Γ > 1/5 behaviour we saw in Figure 1; in particular, the maximum eccentricity of all orbits is found at ω = ±π/2. The range (36) is important because it allows for eccentricity excitation of initially near-circular binaries, which is not possible in the 0 < Γ ≤ 1/5 regime when ǫGR = 0 (see §3.4.1).

Determination of the maximum eccentricity of a given orbit
Our next goal is to calculate the maximum eccentricity emax reached by a binary given the initial conditions (ω0, e0, Θ, Γ, ǫGR). In particular, we wish to know if a binary will reach emax → 1, since this is the regime in which dissipative effects (e.g. GW emission) can become important. For Γ > 1/5, a binary's maximum eccentricity is always found at ω = π/2 regardless of whether its phase space orbit librates or circulates ( Figure 2). Plugging ω = π/2 into H * (ω, j) gives us a depressed quartic equation: We call real roots of equation (37) j(ω = π/2). In the limit ǫGR = 0, equaton (37) reduces to a quadratic for j 2 (ω = π/2) and we recover the solution (18) of Paper II. For ǫGR = 0 the real roots of (37) can still be written down analytically but they are too complicated to be worth presenting here. The minimum angular momentum jmin (corresponding to the maximum eccentricity emax ≡ 1 − j 2 min ) will then be given by the smallest physical root j(ω = π/2), i.e. the smallest root of (37) that satisfies √ Θ < j(ω = π/2) < 1. The situation is slightly more complex for 0 < Γ ≤ 1/5. In this case we must first work out whether an orbit librates or circulates (and if it circulates, to which circulating family it belongs, since it can be above or below the saddle point, as in Figures 2b,c,h). To do so we use the procedure given in §A3 to calculate j(ω = 0), which is the solution to the depressed cubic equation (A7) that results from plugging ω = 0 into H * (ω, j). If the orbit circulates 'below' a librating region and the saddle points then we have jmin = j(ω = 0). Otherwise jmin is found at ω = ±π/2 and we proceed as for Γ > 1/5 by solving equation (37).

Maximum eccentricity achieved by initially near-circular binaries
We can gain further insight and connect to the results of previous LK studies by considering the simplified case of initially near-circular binaries, e0 ≈ 0. Evaluating the integrals of motion H * and Θ with the initial condition e0 = 0 we find Note the lack of ω0 dependence in these constants. Now, for Γ > 1/5 eccentricity is always maximised at ω = π/2, so can be found by solving (37). Plugging (38) into (37) we find that jmin is the solution to the equation In the LK limit of Γ = 1, equation (39) is equivalent to e.g. equation (34) of Fabrycky & Tremaine (2007) or equation (50) of Liu et al. (2015) 6 . Note that jmin = 1 (i.e. emax = 0) is a solution to this equation. It is the correct solution in the special case of a perfectly initially circular orbit, e0 ≡ 0, which necessarily remains circular forever. This is because perfectly circular binaries feel no net torque from the external tide, which can be seen by plugging j = 1 into equation (13).
Meanwhile, an orbit that has e0 infinitesimally larger than zero can have jmin corresponding to a non-trivial solution of (39). This will be the case if and only if the ω = ±π/2 fixed points have not yet disappeared below e = 0 (panels (a)-(d) and (f)-(h) of Figure  1). Because of the constraint (29), a necessary (and for i0 → 90 • , sufficient) requirement for this is ǫGR < 6(1 + 5Γ). In that case the fixed points bound the maximum eccentricity from below, so emax > e f,π/2 ≡ (1 − j 2 f,π/2 ) 1/2 . On the other hand, if ǫGR is large enough that the fixed points have disappeared through e = 0 then we simply have emax = 0 (see panels (e), (i), (j) of Figure 1).
Next we turn to the regime 0 < Γ ≤ 1/5. By consulting Figure 2 one can see that a finite eccentricity is only achieved if ǫGR is sufficiently large that the saddle point at ω = 0 has passed 'down' the (e, ω) phase space and disappeared through e = 0, but also sufficiently small that the ω = ±π/2 fixed points still exist (as in Figure  2d,i). A necessary requirement for this (which is again sufficient in the case i0 → 90 • ) is that (36) be true. Then e is maximised at ω = ±π/2 and jmin is a non-trivial solution to equation (39). On the other hand, if (36) is not satisfied then a binary that starts at e0 ≈ 0 never increases its eccentricity 7 even for i0 = 90 • .
Overall then, we see that for near-circular binaries to reach finite emax we require fixed points to exist in the phase portrait at ω = ±π/2 but not at ω = 0, and this necessarily requires ǫGR to satisfy (36).
In Figure 5 we plot emax as a function of i0 for initially nearcircular binaries. Panels (a)-(c) are for Γ > 1/5 (c.f. Figure  1 − Θ = sin i0 (the highest possible e for an initially near-circular binary, corresponding to j = cos i0). We see that for Γ > 1/5, the effect of increasing ǫGR at a fixed i0 (and therefore a fixed Θ) is always to decrease emax. This is what we would expect by comparing the top and bottom rows of Figure 1. If we consider the most favourable orbital inclination i0 = 90 • then we can easily derive the exact solution to (39). We find that either jmin = 1 (so emax = 0), or that jmin = [(1 + 4ǫGR/ǫstrong) 1/2 − 1]/2 with ǫstrong defined in (32); in the LK limit this result reduces to equation (35) of Fabrycky & Tremaine (2007). Expanding the latter solution for ǫGR/ǫstrong ≪ 1 we find Thus we expect emax → 1 for these favourably inclined binaries when GR is negligible, but also that emax will deviate from 1 considerably when ǫGR starts approaching ǫstrong, which is what we see in Figure 5a,b,c. Obviously this means that the smaller is Γ, the smaller ǫGR needs to be to suppress the very highest eccentricities. Finally we note that there is no magenta curve -corresponding to ǫGR = 30 -in either panel (b) or panel (c). This is because for these Γ values the constraint (36) is violated for ǫGR = 30, so the only possible solution to (39) is emax = 0. Now consider the regime 0 < Γ ≤ 1/5 exhibited in panels (d)-(f). The reader will notice the diminishing number of curves in these panels. Indeed, there is not even a red curve corresponding to ǫGR = 0. This again is a consequence of the fact that for ǫGR = 0, equation (36) cannot be satisfied, so that initially circular orbits achieve no eccentricity excitation (emax = 0).
A related phenomenon is that in panel (e), the green (ǫGR = 3) curve asymptotes to the black dashed line e = e lim as i0 → 0 • . This is also as expected: since Γ = 0.1, equation (36) tells us ǫGR = 3 is precisely the lower bound on GR strength above which initially nearcircular binaries can reach non-zero eccentricities at all i0, since at this value of ǫGR the saddle point crosses e = 0 -see Figure 6i,j,k.
Note that this 0 < Γ ≤ 1/5 behaviour is completely different from that found for near-circular binaries in the Γ > 1/5 regime (and therefore to the known LK results). For Γ > 1/5, taking i0 ≈ 0 • inevitably leads to emax ≈ 0 -in other words there is no eccentricity excitation for initially coplanar (i0 = 0) orbits, regardless of ǫGR. Moreover, for Γ > 1/5 even if a binary can reach a finite maximum eccentricity for ǫGR = 0, increasing ǫGR always decreases this maximum eccentricity. On the contrary, for 0 < Γ ≤ 1/5 reaching a finite emax may be possible even for initially almost coplanar orbits, and a finite ǫGR is actually necessary to trigger the eccentricity excitation starting from a circular orbit. Despite this, comparison of the top and bottom rows of Figure 5 reinforces the idea that the 0 < Γ ≤ 1/5 regime admits far fewer high-eccentricity solutions than Γ > 1/5 as ǫGR is varied.

HIGH ECCENTRICITY BEHAVIOUR
Our next goal is to understand the impact of GR precession on the time dependence of the binary orbital elements in the important limit of very high eccentricity, e → 1. This limit is relevant in a variety of astrophysical contexts. For instance the dramatic reduction of the binary pericentre distance that occurs when e approaches unity can trigger short-range effects such as tidal dissipation (leading to hot Jupiter formation), GW emission (leading to compact object mergers), and so on. Thus we wish to understand in detail how GR precession affects not only the maximum eccentricity emax, but also the behaviour of e(t) and other orbital elements in the vicinity of emax.
In §3.4 we explained how to find emax for arbitrary Γ > 0, initial conditions (e0, i0, ω0), and value of ǫGR. Here we will examine the solutions quantitatively in the high eccentricity limit, and we will explore the time spent near highest eccentricity. To this end, we will make extensive use of equation (15) which tells us dj/dt as a function of j. It is important to note that the solutions for extrema of j at ω = 0 and ω = ±π/2 are all contained within (15). Indeed, setting the first square bracket inside the square root in (15) to zero gives the depressed quartic equation (37) whose roots correspond to extrema of j at ω = ±π/2, i.e. what we have so far called j(ω = ±π/2). Setting the other square bracket to zero gives the depressed cubic (A7)-(A9) which determines the roots at ω = 0, i.e. what we called j(ω = 0).
In this section we will focus on situations in which emax is achieved at ω = ±π/2, since this is the most common prerequisite for e → 1 ( § §3.1-3.2). The rare cases in which e approaches unity at ω = 0 are covered in Appendix B.

Phase space behaviour for
We are interested in trajectories that start with initial eccentricity e0 not close to unity, and that are capable of reaching extremely high eccentricities emax → 1, i.e. jmin → 0. For this to be possible Note that for initially circular orbits to reach a non-zero emax we require fixed points to exist in the phase portrait at ω = ±π/2 but not at ω = 0; for i 0 ≈ 90 • this corresponds to 6(1 − 5Γ) < ǫ GR < 6(1 + 5Γ) -see equation (36). Note also that for Γ < 1/5, eccentricity excitation of near-circular binaries may be possible regardless of inclination, even when i 0 = 0 • , as in panel (e). a necessary condition is that Θ ≪ 1, owing to the constraint (9). Hence it is important to understand the regime Θ ≪ 1 in detail. In Figure 6 we show phase portraits for Γ = 0.5 (top row) and Γ = 0.1 (bottom row), this time fixing Θ = 10 −3 in both cases, and adding in extra dashed contours 8 with the value H * (ω = ±π/2, e = 0.9). Note that on the vertical axis we now plot 1 − e on a logarithmic scale, with eccentricity still increasing vertically as in Figures 1, 2. This allows us to see in detail how trajectories separate from e ≈ e lim as we increase ǫGR.
In these plots, Θ is sufficiently small that to a very good approximation the requirement for fixed points at ω = ±π/2 to exist is just ǫGR < 2ǫstrong (equation (33)). This critical value is surpassed in panel (l) because there ǫGR = 10 while 2ǫstrong = 9, which is why all fixed points have disappeared. Meanwhile the criterion for a saddle point to exist at ω = 0 (equation (35)) for Γ = 0.1 and Θ = 10 −3 is approximately 10 −5 < ǫGR < 3. This is consistent with what we see in panels (f)-(l) -note in particular the transitional point ǫGR = 3 in panel (j).
Comparing the top and bottom rows of Figure 6, one observes a striking difference between behaviour in the Γ > 1/5 and 0 < Γ ≤ 1/5 dynamical regimes. For Γ = 0.5 > 1/5, an initially nearcircular binary can be driven to very high eccentricity ( 0.99) even for ǫGR = 1.0 (panel (d)). Conversely, for Γ = 0.1 < 1/5 the phase space structure simply does not allow such behaviour (panels (f)-(i)). It is therefore unsurprising that one finds fewer cluster-tide driven compact object mergers from systems such as globular clusters that have a relatively high fraction of binaries in the 0 ≤ Γ ≤ 1/5 regime (Hamilton & Rafikov 2019c).

High eccentricity behaviour for ǫGR = 0
Before embarking on a full study of high eccentricity evolution for arbitrary ǫGR, we first consider the case ǫGR = 0. In that case the non-zero roots of the polynomial on the right hand side of (15) are j±, j0, one of which will correspond to the minimum angular momentum jmin. Then we can integrate (15) with ǫGR = 0 to find t(j); the resulting expression involves an incomplete elliptical integral of the first kind (see §2.6 of Paper II for the general Γ case, and Vashkov'yak (1999); Kinoshita & Nakai (2007) in the LK case of Γ = 1). Next, assuming that j 2 ≪ 1, we can expand this elliptical integral to find 9 (see §9.2 of Paper II) j1, j2 are the two roots not corresponding to jmin, and τ is a characteristic secular timescale which is independent of e0, i0, ω0: Using the definitions of C and L one can show that τ is, up to constant factors, the same at tsec defined after equation (7). Note we have taken the origin of the time coordinate to coincide with j = jmin. Clearly tmin is the characteristic evolution timescale in the vicinity of jmin, i.e. the time it takes for j to change from jmin to √ 2jmin. Note that the solution (41) is quadratic in t for t tmin and linear when t tmin, as long as j remains ≪ 1. It provides a 8 In addition to the dashed contours already included in Figures 1 and 2. 9 Note that one can get the same result simply by expanding the right hand side of (15) for j ≪ 1.
better approximation to j(t) over a wider interval of time near the peak eccentricity than the purely quadratic approximation adopted by Randall & Xianyu (2018), in their calculation of the GW energy emitted by a binary undergoing LK oscillations ( §C3).

Modifications brought about by finite ǫGR
Before we proceed to examine the j(t) behaviour, it is important to realise that including a finite ǫGR affects the right hand side of (15), and therefore the value of jmin, in two distinct ways. First, there is the obvious explicit dependence on ǫGR that appears twice in equation (15). Second, there is also an implicit dependence on ǫGR in (15) through the values of j± and j0 (see equations (16), (17)). We will now discuss this implicit dependence, and then use the results to understand j(t) behaviour in different asymptotic ǫGR regimes.
In the limit Θ ≪ 1, and assuming that e0 is not too close to 1 and Γ is not too close to 1/5, equations (18), (19) tell us that Equation (43) implies that Σ, and hence j 2 ± , will only be significantly modified by GR if ǫGR ǫstrong, in agreement with what we saw in Figures 1, 2 & 6. In this case, a perturbative approach around the non-GR solution will fail. We therefore say that any binary with ǫGR ǫstrong exists in the regime of 'strong GR', which we explore in §4.5. Conversely, if ǫGR is in what we will call the 'weak-to-moderate GR' regime: then Σ ∼ 1, and j± will stay close to their non-GR values for Θ ≪ Σ ∼ 1, namely (see equation (16)) In other words the relative perturbations to j± induced by GR can be neglected. Note that the weak-to-moderate GR regime (44) already encompasses the very weak GR regime introduced in §3.3.1. In § §4.4.1-4.4.2 we will further delineate distinct 'weak GR' and 'moderate GR' regimes. Also, using equations (17), (19) it is easy to show that the absolute change to j0 incurred by including GR will be small (≪ 1) whenever Note that for Γ not close to 1/5 and e0 not close to unity, the condition (46) is automatically guaranteed by the weak-to-moderate GR condition (44). In that case the relative perturbation to j0 due to GR precession can be neglected (if j0 ∼ 1).

High-e behaviour in the weak-to-moderate GR limit
In the non-GR limit (ǫGR = 0), for Γ > 0 the vast majority of phase space trajectories that are capable of reaching very high eccentricities reach them at 10 ω = ±π/2. As shown in Paper II, for Γ > 0 the corresponding minimum angular momentum for these orbits is always jmin = j− ≪ 1.
We now want to see what happens to (15) for finite ǫGR ≪ ǫstrong. From the discussion in §4.3 we expect that we may neglect j 2 compared to j 2 + , j 2 0 in this regime. As a result we can write dj dt where we defined the following dimensionless numbers: with To get the second equality in (50) we used the approximation (45). Both ǫ weak and γ are manifestly positive in the weak-to-moderate GR regime given Γ > 0. Except in pathological cases, σ is also positive for the regimes we are interested in here 11 . To find the minimum j at ω = ±π/2 we require the right hand side of (47) to equal zero, which, as we mentioned earlier, means that the first square bracket inside the square root must vanish. This gives a quadratic equation for jmin, the only meaningful (positive) solution to which is Equations (47) and (51) work as long as Θ ≪ 1 and ǫGR is in the weak-to-moderate GR regime, i.e. satisfies (44) and (46). Equation (51) has been used by several authors in the LK limit of Γ = 1 -see §5.2. Importantly, it allows us to write down a solution for the maximum eccentricity reached by initially nearcircular binaries in the weak-to-moderate GR regime. Indeed, let us put Θ = cos 2 i0 and assume Θ ≪ 1 (i.e. i0 ≈ 90 • ) so that the binary is capable of reaching very high eccentricity. Then j+ ≈ 1 and from (51) we find Note that this result can also be derived by solving equation (39) in the limit j ≪ 1. Moreover, in the LK limit Γ = 1 we recover from (52) a well-known result, identical to 12 e.g. equation (8) of Miller & Hamilton (2002) and equation (52) of Liu et al. (2015).
It is now instructive to investigate separately the high-e behaviour in the asymptotic regimes of weak and moderate GR precession (still assuming eccentricity is maximised at ω = ±π/2).
11 This is true because (5Γ − 1)j 2 0 is positive in the ǫ GR = 0 limit for all the cases we care about, namely any orbit with Γ > 1/5 and librating orbits with 0 < Γ ≤ 1/5. The inclusion of GR subtracts from (5Γ − 1)j 2 0 by an amount ǫ GR /(3 1 − e 2 0 ). For j 0 ∼ 1 and e 2 0 ≪ 1, this modification will not make (5Γ − 1)j 2 0 negative as long as (46) is satisfied. 12 Note that our definition of ǫ GR differs from what Miller & Hamilton (2002) call θ PN and what Liu et al. (2015) call ε GR . Our ǫ GR is defined for any outer orbit in any axisymmetric potential, whereas their parameters are defined only in the Keplerian (LK) limit. In this limit, ǫ GR = 6θ PN = 16ε GR .  (55) for j(t) near the eccentricity peak for various values of σ (equation (49)). The horizontal dotted line shows j/j min = √ 2 and the vertical dotted line shows t = t min . Panel (b) shows β weak (σ), which is the time (in units of t min ) over which the binary's j/j min changes from 1 to √ 2, defined by setting j/j min = √ 2 in the right hand side of (55). Note that both axes are on a logarithmic scale in this panel. A dashed magenta line shows the scaling β weak ∝ σ −1/2 for σ ≫ 1.

Weak GR, ǫGR ≪ ǫ weak
In the asymptotic regime of weak GR, defined by ǫGR ≪ ǫ weak , the solution (51) becomes approximately In other words, GR causes only a slight perturbation of jmin away from the non-GR value of j− at the relative level ǫGR/ǫ weak ≪ 1.
To determine the time dependence of j(t) in the vicinity of jmin, we make use of the weak GR assumption to drop the γ term in the first square bracket in (47). The result is Integration of (54) gives an implicit solution for j(t) in the form where tmin is defined in equation (41). In Figure 7a we plot the implicit solution for j/jmin as a function of t/tmin for various values of σ. Usually in the weak GR regime we expect σ ≪ 1 (equation (49)). In that case the term σj− in the final square bracket in (54) can also be dropped compared to j. Then equation (54) takes the same functional form as its non-GR analogue; integrating, we get a solution j(t) in precisely the form (41) with 13 jmin → j− and j1, j2 → j+, j0. This is reflected in Figure 7a, in which the black line (σ = 0) is exactly the non-GR result from (41), and as expected j = √ 2jmin coincides with t = tmin in that case. However the assumption σ ≪ 1 may not always be valid even in the weak GR regime, particularly if Γ is near 1/5. Figure 7a shows that as we increase σ the behaviour of j(t) becomes more sharply peaked around jmin (when time is measured in units of tmin), although against this trend one must remember that to change σ is to change one or more of ǫGR, Γ, j0 and j−, any of which will modify tmin. We are particularly interested in the value of t weak min (σ), which is the time it takes for j to go from jmin to √ 2jmin in the weak GR regime, to compare with the solution (41). By setting j/jmin = √ 2 on the right hand side of (55) and t = t weak min on the left, we find t weak min (σ) = β weak (σ)tmin, where β weak (σ) is plotted as a function of σ in Figure 7b. As expected β weak → 1 for σ → 0, i.e. in the limit of negligible GR precession. For finite GR, typical values of β weak are ∼ 1 except for very large σ 10. For σ ≫ 1 we see that β weak falls off like ∼ σ −1/2 .
In Figures C2, C3, C5 and C6 we compare the weak GR solution for j(t), namely equation (55), to direct numerical integration of the DA equations of motion (12), (13), for binaries in different dynamical regimes. Full details are given in §C3; here we only note that the values of the key quantities Γ, ǫGR, ǫ weak , σ, etc. are shown at the top of each figure. In every example, panel (a) shows log 10 (1 − e) behaviour in the vicinity of peak eccentricity, while panel (b) shows the same thing zoomed out over a much longer time interval 14 . The weak GR solution for j(t) (equation (55)) is plotted in panels (a) and (b) with a dashed green line, while the numerical solution is shown with a solid blue line. We see that for ǫGR ≪ ǫ weak (Figures C2, C3) this weak GR solution works very well, but that substantial errors begin to set in when ǫGR approaches ǫ weak (Figures C5, C6). Finally, in each of these plots we also show with red dashed lines the 'analytic' solution (C2), which coincides with (41) provided j 4 min /Θ ≪ 1. As we have already stated, in the weak GR regime j(t) takes the form (41) provided that σ ≪ 1, so it is unsurprising that in the plot with very small σ ( Figure C2) this analytic solution (C2, red dashed lines) overlaps with the weak GR solution (55, green dashed lines).

Moderate GR, ǫ weak ≪ ǫGR ≪ ǫstrong
Perhaps more interesting is the asymptotic regime of moderate GR, defined as ǫ weak ≪ ǫGR ≪ ǫstrong. In this regime one finds from (51) that i.e. a significant perturbation of jmin away from j−, resulting in a significantly reduced maximum eccentricity emax.
To determine the time dependence of j(t) around jmin we neglect j 2 − compared to j 2 in the first square bracket in (47) and find dj dt 13 Note that j ± , j 0 depend on ǫ GR through (16)-(17) only weakly, at the relative level O(ǫ GR /ǫ weak ). 14 Note however that on the horizontal axis we plot time in units of t ′ min (equation (C3)) rather than t min , as explained in Appendix C.  (60)). In panel (b) we show β mod (κ), which is the time (in units of t min ) over which j changes from j min to √ 2j min in this regime. Dotted lines show κ = 1 and β mod = 1, while a dashed magenta line shows the scaling β mod ∝ κ −1/2 for κ ≫ 1.
Integration of (58) gives an implicit solution for j(t) in the form where Note that tmin in (59) is still defined by equation (41) but taking jmin equal to its GR-modified value, namely γj−. In Figure 8a we plot this implicit solution for various values of κ. Note that κ = 1 (red line) gives precisely the solution j/jmin in the form (41), and so unsurprisingly j/jmin = √ 2 coincides with t/tmin = 1 in that case. As we increase κ we see that the time spent near the minimum j decreases (when measured in units of tmin, which itself also depends on κ).
Analogous to §4.4.1, by setting j/jmin = √ 2 on the right hand side of (59) and t = t mod min on the left, we find that the time for j to increase from jmin to √ 2jmin in the moderate GR regime is where β mod (κ) is plotted as a function of κ in Figure 8b. Clearly when κ ∼ 1 (which is true for Γ not too close to 1/5) we have β ∼ 1 and so t GR min ∼ tmin. But for κ ≫ 1 the time spent in the high eccentricity state is somewhat reduced, with a scaling t GR min ∝ κ −1/2 tmin. However for this to be a significant effect requires rather extreme values of κ (i.e. Γ very close to 1/5).
Finally, in panels (a) and (b) of Figures C4, C7 we compare the moderate GR solution (59), shown with dashed cyan lines, to direct numerical integration of the DA equations of motion, shown in solid blue. In both cases the moderate GR solution provides an excellent fit to the numerical result despite ǫGR only being slightly larger than ǫ weak .

High-e behaviour in the strong GR limit
The final asymptotic case to consider is that of strong GR, ǫGR ≫ ǫstrong. This regime is important for understanding the later stages of evolution of shrinking compact object binaries. Indeed, GW emission should eventually bring any merging binary to a small enough semimajor axis to put it in this regime.
In this limit GR precession is the dominant effect, exceeding the secular effects of the external tide -see e.g. panels (e) and (j) of Figures 1 and 2. Thus we anticipate that at high eccentricity the lowest order solution will be one of constant eccentricity and uniform prograde precession: whereωGR(0) = (C/L)ǫGR/j 2 (0) -see equation (1). High eccentricity can therefore only be achieved if j(0) ≪ 1 to start with. This is actually a highly relevant scenario in practice because it is at very high eccentricity that GW emission, and hence the shrinkage of a and the growth of ǫGR, is concentrated. Binaries periodically torqued to very high eccentricity by cluster tides eventually become trapped in a highly eccentric orbit as they enter the strong GR regime (Hamilton & Rafikov, in prep.). Their phase space trajectories are then well described by (62). Interestingly, the minimum angular momentum jmin = j(0) predicted by the solution (62) can still be described by the expression (51) in the limit ǫGR ≫ ǫstrong. To see this we note that for ǫGR ≫ ǫstrong equation (51) gives jmin ≈ ǫGR/(j 2 + ǫstrong). Then we take the expression (45) for j 2 + , and substitute into it the value of Σ we get by taking ǫGR ≫ 1 in (18)-(19), namely Σ ≈ ǫGR(1 − e 2 0 ) −1/2 /6. Putting these pieces together we find jmin ≈ 1 − e 2 0 = j(0). Thus the solution (51) interpolates smoothly between the different asymptotic GR regimes.

Evolution of ω(t)
and Ω(t) as e → 1 So far we focused on understanding the behaviour of j(t). Once this is determined one can understand the evolution of other orbital elements as well. In particular, since the Hamiltonian (2) is conserved, we can use equations (10)-(11) to express cos 2ω entirely in terms of j(t) and conserved quantities, leading to an explicit analytical expression for ω(t). Finally one can plug this cos 2ω(t) and j(t) into the equation of motion (14) for Ω(t) and integrate the result. Together with Jz(t) = Jz(0) = const., this constitutes a complete solution to the DA, test-particle quadrupole problem with GR precession.
Unfortunately this proposed solution for ω(t) and Ω(t) is very messy for arbitrary values of j. Luckily, in the high eccentricity regime j ≪ 1, one can make substantial progress by (i) making the additional (and often well-justified) assumption given by equation (C1) and (ii) adopting the ansatz 15 (41) for j(t). Then, as we show in Appendix C, one can derive relatively simple explicit analytical solutions for ω(t) and Ω(t). These solutions work very well as long as equation (41) is a good approximation to the j(t) behaviour near the peak eccentricity, which we verified numerically ( §C3). To our knowledge, an explicit high-e solution of this form accounting for GR precession has not been derived before even for the LK problem. It can be used e.g. when exploring the short-timescale (i.e. non-DA) effects near the peak eccentricity, which are important for accurate calculation of the GW wave production.

DISCUSSION
In this paper we have studied the impact of 1PN GR precession on secular evolution of binaries perturbed by cluster tides. A single dimensionless number Γ effectively encompasses all information about the particular tidal potential and the binary's outer orbit within that potential. Meanwhile the relative strength of GR precession compared to external tides is characterised by the dimensionless number ǫGR (equation (6)).
In the main body of the paper we only discussed the systems with Γ > 0. Although the resulting dynamics are significantly complicated by bifurcations that occur at Γ = ±1/5, 0, for Γ > 1/5 our qualitative results are intuitive, falling in line with those gleaned from previous LK (Γ = 1) studies that accounted for GR precession. However, for 0 < Γ ≤ 1/5 we uncover a completely new pattern of secular evolution, which we have characterised in detail. Secular dynamics of binaries with negative Γ (possible for binaries on highly inclined outer orbits in strongly non-spherical potentials, see Paper I) and non-zero ǫGR is covered in Appendix D. As mentioned in §2.1, the Γ ≤ 0 regime splits into two further regimes, namely −1/5 < Γ ≤ 0 and Γ ≤ −1/5. We generally find that in both of these regimes the resulting phase space structures and maximum eccentricity behaviour are considerably more complex and counterintuitive than for Γ > 0.
Furthermore, we have explored the evolution of binary orbital elements in the limit of very high eccentricity ( §4). This investigation revealed a number of distinct dynamical regimes that are classified according to the value of ǫGR. In §5.1 we summarise and systematise these regimes based on their physical characteristics. In §5.2 we compare our work to the existing LK literature, and in §5.3 we discuss the limitations of our work.
Based on our findings in § §3-4 we now provide a description of the basic features of each regime. Very weak GR In this limit (below the dashed curves in Figure  9) GR precession has essentially no effect on the dynamics for Γ > 1/5. GR is too weak to affect either the locations of the fixed points at ω = ±π/2, which are given by equation (27), or the maximum eccentricity reached by the binary at the same ω in the course of its tide-driven secular evolution, given by equation (53). Thus all Γ > 1/5 results of Paper II, which were derived for ǫGR = 0, are valid. However, for 0 < Γ ≤ 1/5 an important modification arises if 6(1 − 5Γ)Θ 3/2 < ǫGR ǫ π/2 , which is that saddle points appear at ω = 0, π. These saddles do not exist for ǫGR = 0 (see equation (35) and §3.3.2), but they do change the maximum eccentricity reached by the binary -see Appendix B. Weak GR In this regime (between the dashed and dotted curves in Figure 9) GR precession starts to modify the j locations of the fixed points at ω = ±π/2, which are now given by equation (28). At the same time, GR precession does not appreciably change jmin (or equivalently emax), which stays close to its ǫGR = 0 value j− (see equation (53)). If σ ≫ 1 then GR also modifies the time spent in the high eccentricity state (equation (56)). Moderate GR In this regime (between the dotted and solid curves in Figure 9) GR precession modifies not only the locations of the fixed points but also the values of jmin (and hence of emax), now given in equation (57). GR also modifies the time spent in the high-e state (equation (61)). In other words, in this regime GR precession presents an efficient barrier suppressing the maximum eccentricity reached by the binary in the course of its secular evolution. Strong GR In this limit (above the solid curves in Figure 9) GR precession dominates the binary dynamics at all times; the quantities j± get significantly affected by GR precession (equation (43)) and all fixed points in the phase portrait disappear (equations (33), (35)). Cluster tides drive only very small eccentricity oscillations on top of uniform GR precession, so that e is roughly constant -see equation (62).
In Table 2 we summarise the main features of the asymptotic ǫGR regimes that we have found in this and previous sections.
This regime separation sheds light on the physical meaning of the characteristic ǫGR values introduced in this work. We first note that the GR precession rate (1) can be written aṡ ωGR(j) =ωGR|e=0j −2 ∼ ǫGRt −1 sec j −2 (see the definition (6)). Then, for a cluster tide-driven process occurring on a characteristic timescale t ch , GR precession becomes important wheṅ ωGR(j)t ch ∼ 1, i.e. ǫGR t ch tsec j −2 ∼ 1.
If ǫGR satisfies (64), or exceeds that value, then GR breaks the coherence of the tidal torque over the timescale ∼ t ch . This means that GR precession substantially interferes with the secular evolution. We now demonstrate how this simple physical argument leads one to the critical values ǫstrong, ǫ weak and ǫ π/2 . First, in the strong GR regime we expect GR precession to dominate binary evolution at all times, even for near-circular orbits. Setting t ch ∼ tsec and j ∼ 1 we obtain ǫGR ∼ 1, which is consistent with the definition (32) of ǫstrong up to a numerical coefficient.
Second, in the moderate GR regime, we anticipate that GR precession will present an effective barrier that stops the decrease of j if t ch is the characteristic timescale of secular evolution near the eccentricity peak. In §4 we find quite generally this timescale to be t ch ∼ tmin ∼ jmintsec -see e.g. equations (41) and (61). Plugging this into the condition (64) and evaluatingωGR at jmin we immediately find that jmin ∼ ǫGR, in agreement with equation (57). When the GR barrier first emerges at the transition between weak and moderate regimes, jmin is still well approximated by the ǫGR = 0 solution j− ∼ Θ 1/2 (see equation (45)). As a result, the ǫGR value corresponding to this transition is ∼ Θ 1/2 , in agreement with the definition (50) of ǫ weak . Third, we expect fixed points in the phase portrait at ω = ±π/2 to be substantially displaced by GR precession whenωGR(j f ) becomes comparable to the characteristic secular frequencyω of libration around a fixed point. Since we are interested in the displacement of j f by an amount ∼ j f , we take thisω from the H ⋆ = const. contour centred on the ω = ±π/2 fixed point and with vertical extent ∼ j f . Plugging j ∼ j f into the equation (12) and using the expression (24) for j f we findω ∼ Θ 1/4 t −1 sec , so that in this case t ch ∼ Θ −1/4 tsec. Substituting this into the condition (64) and again setting j ∼ j f ∼ Θ 1/4 we find ǫGR ∼ Θ 3/4 for the transition between the weak and very weak GR regimes. This agrees with the definition of ǫ π/2 in equation (26).
Note that while these considerations allow us to understand the scalings of characteristic ǫGR values with Θ, one still needs the full analysis presented in § §3-4 to obtain the numerical coefficients, which are actually quite important. Indeed, equations (26), (50), and (32) feature constant numerical factors which can substantially exceed unity, especially for the LK case of Γ = 1.

Relation to LK studies
Many authors who studied the LK mechanism and its applications have included 1PN GR precession in their calculations. The maximum eccentricity of an initially near-circular binary undergoing LK oscillations (i.e. the Γ = 1 limit of §3.4.1) was derived by Miller & Hamilton (2002); Blaes et al. (2002); Wen (2003); Fabrycky & Tremaine (2007); Liu et al. (2015). Of these, Fabrycky & Tremaine (2007) and Liu et al. (2015) also produced plots very similar to Figure 5 that show how increasing ǫGR decreases the maximum eccentricity achieved by initially near-circular binaries. Various authors have derived equations identical to, or very similar to, the quartic (37) and the weak-to-moderate maximum eccentricity solution (51) in the LK limit -see for instance equation (A7) of Blaes et al. (2002), equation (8) of Wen (2003), equation (A6) of Veras & Ford (2010), and equations (64)-(65) of Grishin et al. (2018). Of course, because these studies only work with Γ = 1, the rather non-intuitive behaviour for Γ < 1/5 revealed in §3.4.1 and Appendix D3.1 has not been unveiled before. Moreover, to our knowledge no previous study has presented a clear classification of the different ǫGR regimes (which we do in §5.1), even in LK theory.
The quantitative results in the aforementioned papers have been employed in many practical calculations. Typically one simply adds the term (1) to the singly-or doubly-averaged equations of motion along with any other short range forces or higher PN effects. In population synthesis calculations of compact object mergers (Antonini & Perets 2012;Antonini et al. 2014;Silsbee & Tremaine 2017;Liu & Lai 2018) one often puts a sensible lower limit on the semimajor axis distribution below which GR is so strong that sufficient eccentricity excitation is impossible. As explained in §2.1 of Rodriguez & Antonini (2018) there are at least two ways to decide when GR dominates. One method is to take jmin corresponding to the pericentre distance that needs to be reached according to the problem at hand, and then equateωGR(jmin) with the precession rate due to the tidal perturbations (see their equation (29)), which corresponds to equation (54) in Paper II. As discussed in §5.1, this method would set a rough upper limit of ǫGR jmin; see equation (57) for a more accurate expression. A second method is to demand that ω = ±π/2 fixed points do exist in the phase portrait (Fabrycky & Tremaine 2007) allowing for substantial eccentricity excitation to occur starting from the near-circular orbits, which is equivalent to ǫGR ǫstrong. However, this is a rather weak requirement, and does not guarantee that the majority of systems with such ǫGR would reach the required jmin -many of them will be stopped by the GR barrier at eccentricities much lower than needed. The former method of setting an upper limit on ǫGR is typically more stringent and allows more efficient selection of systems for Monte Carlo population synthesis (see Figure 9).
With regard to phase space structure, the only study we know of that resembles our §3 is that by Iwasa & Seto (2016). They considered a hierarchical triple consisting of a star on an orbit around a supermassive black hole (SMBH), with another massive black hole also orbiting the SMBH on a much larger, circular orbit and acting as the perturber of the star-SMBH 'binary'. Their §C provides a brief explanation of the phase space behaviour as a parameter they call γ, which is equivalent to our ǫGR/3, is varied. Since the LK problem has Γ = 1 > 1/5, their Figure 2 is qualitatively the same as our Figure 1.

Approximations and limitations
To derive the Hamiltonian (4) we truncated the perturbing tidal potential at the quadrupole level. This is justified if the semimajor axis of the binary is much smaller than the typical outer orbital radius. Next order corrections to the perturbing potential -so called octupole terms -are routinely accounted for in LK studies (Naoz et al. 2013;Will 2017). In Appendix E of Paper I we provide the octupole correction to (4) for arbitrary Γ. When octupoleorder effects are important, the maximum eccentricity can actually be increased by GR precession (Ford et al. 2000;Naoz et al. 2013;Antonini et al. 2014). However for the applications we have in mind, e.g. a compact object binary of a ∼ 10AU orbiting a stellar cluster at ∼ 1pc, octupole corrections are negligible.
We also employed the test particle approximation, which is valid if the outer orbit contains much more angular momentum than the inner orbit. One can relax the test particle approximation: in particular, this is often necessary for weakly-hierarchical triples. Anderson et al. (2017) made a detailed study of the 'inclination window' that allows fixed points to exist in the (quadrupole) LK phase space for different ǫGR, as one varies the ratio of inner to outer orbital angular momenta. We recover their results in the test particle limit valid for our applications. We also assumed the validity of the DA approximation, the smallness of short-timescale fluctuations ('singly-averaged effects'), etc., all of which are liable to break down at very high eccentricity. For a full discussion of these issues see Papers I and II.
Finally, several of the results derived at very high eccentricity ( §4) are rather delicate when Γ is close to ±1/5 or when the binary's phase space trajectory is close to a separatrix. These are not major caveats; for instance, in a given stellar cluster potential only a small fraction of binaries will have Γ values close enough to 1/5 to be affected (Paper I).

SUMMARY
In this paper we completed our investigation of doubly-averaged (test-particle quadrupole) cluster tide-driven binary dynamics in the presence of 1PN general relativistic pericentre precession. Throughout, we parameterised the strength of GR precession relative to tides using the dimensionless number ǫGR (equation (6)). We can summarise our results as follows: • We investigated the effect of non-zero ǫGR on phase space morphology. For values of ǫGR much less than a critical value ǫstrong, bifurcations in the dynamics happen at Γ = ±1/5, 0, so that we must consider four Γ regimes separately. We found that for Γ ≤ 1/5 a non-zero ǫGR can lead to entirely new phase space morphologies, including (previously undiscovered) fixed points located at ω = 0, ±π.
• We presented general recipes for computing the locations of fixed points in the phase portrait, for determining whether a given phase space trajectory librates or circulates, and for finding its maximum eccentricity, for arbitrary ǫGR.
• We considered how the maximum eccentricity reached by an initially circular binary is affected by GR precession. For Γ > 1/5 the intuitive picture holds that a larger ǫGR leads to a lower maximum eccentricity, but this is not always the case for Γ ≤ 1/5.
• We delineated four distinct regimes of secular evolution with GR precession depending on the value of ǫGR -'strong GR', 'moderate GR', 'weak GR', and 'very weak GR' -and provided physical justification for transitions between them.
• We also studied secular evolution with GR precession in the limit of very high eccentricity. We determined the GR-induced modifications to the minimum angular momentum jmin achieved by the binary and the time dependence of j(t) near the eccentricity peak, which can be rather non-trivial.
• We also provided an approximate analytic description for the evolution of other orbital elements -pericentre and nodal angles -near the eccentricity peak, accounting for the GR precession.
In upcoming work we will apply the results of this paper to understand the long-term evolution of compact object binaries due to GW emission, leading to their mergers and the production of LIGO-Virgo GW sources. Furthermore, these results will inform future studies on the effect of short-timescale fluctuations ('singly-averaged effects') on binaries undergoing cluster tide-driven secular evolution, as well as the population synthesis calculations of merger rates.
Next we wish to instead fix Γ and ǫGR and look for the resulting condition on Θ that allows the fixed points to exist. To begin with, we look for the critical Θ values for which j f,π/2 = 1 and j f,π/2 = √ Θ. The former is Θ1, the expression for which is given in equation (31), and the latter is Θ2 which is determined implicitly through the equation For Γ > 0 this equation can have meaningful (0 ≤ Θ2 ≤ 1) solutions only if 10Γ/(1 + 5Γ) ≤ 1, i.e. if Γ ≤ 1/5. For Γ > 1/5 the value of Θ2 has no physical significance. Next, to determine the proper constraint on Θ we begin by setting Θ = 0 in (A1) -we see that the fixed point exists and has value j f,π/2 = (ǫGR/[6(1 + 5Γ)]) 1/3 , so that j f,π/2 > √ Θ at Θ = 0. As we increase Θ, there are two possibilities. The first is that j f,π/2 increases steeply enough that it reaches unity (at Θ = Θ1) before it intersects √ Θ. In this case the constraint on Θ for fixed points to exist is Θ < Θ1. The second scenario is that j f,π/2 intersects √ Θ (at Θ2) before it reaches unity. Then for the inequality √ Θ < j f,π/2 to be satisfied one needs Θ < Θ2. Of course, since Θ2 is only physically meaningful for Γ ≤ 1/5, the second scenario can only occur in that Γ regime. This reasoning leads us to the constraint (30) for fixed points to exist at ω = ±π/2. In the limit ǫGR → 0 this constraint reduces to the non-GR constraint (25).
In summary, for Γ > 0, fixed points at ω = ±π/2 exist if the constraints on both Θ and ǫGR are satisfied simultaneously; thus they exist in the sub-volume of (Γ, Θ, ǫGR) space bounded by the inequalities (29), (30). We can use this information to understand Figure 4 in more detail. Recall that in this figure we plotted Θmax(Γ), namely the maximum value of Θ for which fixed points exist at ω = ±π/2 for a given ǫGR (equation (30)). We now seek to understand separately the behaviour for Γ > 1/5 and 0 < Γ ≤ 1/5.

A2 Fixed points at ω = 0
For Γ ≤ 1/5 we found (e.g. Figure 2) that hitherto undiscovered fixed points could arise at ω = 0. To find the eccentricity of these fixed points we plug ω = 0 into dω/dt = 0 using equation (12). The result is a cubic equation for j with no quadratic or linear terms. The solution is j = j f,0 with j f,0 given in equation (34), which is physically meaningful only for 0 < Γ ≤ 1/5 (i.e. fixed points at ω = 0 do not exist for Γ > 1/5). The determinant of the Hessian matrix of H * (ω, j) evaluated at the point (0, j f,0 ) is equal to (A6) Clearly for 0 < Γ ≤ 1/5 the determinant (A6) is negative whenever the fixed point j f,0 exists, so (ω, e) = (0, e f,0 ) is necessarily a saddle point in the phase portrait, consistent with Figures 2b,c,h.
A3 Does a given orbit librate or circulate?
Here we show how to determine whether a phase space trajectory is librating or circulating, given Γ, Θ, ǫGR and the initial phase space coordinates (ω0, e0). For Γ > 0, librating orbits cannot cross ω = 0 and so any trajectory that passes through ω = 0 must be circulating 16 . Therefore we can figure out whether an orbit librates or circulates by determining whether it crosses ω = 0.
Plugging ω = 0 into H * (ω, j) gives us a depressed cubic polynomial: where We call the real roots of this polynomial j(ω = 0). In the limit ǫGR = 0 we have a0 = 0 and we recover the expression for j(ω = 0) from equation (14) of Paper II. For ǫGR = 0, the nature of the roots of (A7) depends on the sign of the discriminant We can evaluate ∆ given (ω0, e0, Θ, Γ, ǫGR). In particular, we evaluate H * using these initial conditions, guaranteeing that the value of H * we plug in to (A8) is physically meaningful. There are then a few different cases to consider: • If ∆ < 0, equation (A7) has one real root, which may or may not be physical. If a1 < 0 then this root can be written as whereas for a1 > 0 it is given by Once j(ω = 0) has been determined, the orbit circulates if √ Θ < j(ω = 0) < 1, and librates otherwise.
• If ∆ > 0 (which necessarily requires a1 < 0) there are three distinct real roots, and they can be expressed as for k = 0, 1, 2. We also know that that the product of the three real roots of (A7) is −a0 = ǫGR/[3(5Γ − 1)] and their sum is 0. For Γ > 1/5 this implies that two roots (namely k = 1, 2) must be negative and one (k = 0) positive. Thus the orbit circulates if the k = 0 solution lies in ( √ Θ, 1), and librates otherwise. For 0 < Γ ≤ 1/5, one root (k = 2) must be negative and the other two (k = 0, 1) positive. If either or both of the two positive roots lies in ( √ Θ, 1) then the orbit circulates. If neither of them do then it librates. The case of both positive roots lying in ( √ Θ, 1) corresponds to two coexisting families of circulating orbits that share values of H * , one above e f,0 and one below, as in Figure 2b,c,h. To determine the family of circulating orbits to which the trajectory belongs we compare its initial eccentricity e0 with that of the saddle point e f,0 . If e0 > e f,0 then the orbit circulates in the family 'above' the saddle point, and vice versa.

APPENDIX B: HIGH-ECCENTRICITY BEHAVIOUR FOR ORBITS WHOSE ECCENTRICITY MAXIMA ARE FOUND AT ω = 0
When GR is switched off, the only binaries whose eccentricity is maximised at ω = 0 are those on circulating phase space trajectories in the regime 0 < Γ ≤ 1/5 (e.g. Figure 6f; see Paper II for a thorough discussion). The minimum j in this case is jmin = j0 (Paper II), which is given in equation (17). For this to correspond to very high eccentricity one needs D ≈ 1, and the orbit must sit very close to the separatrix between librating and circulating orbits, which can be hard to achieve in practice.
Nevertheless, suppose j0 ∼ Θ 1/2 ≪ 1 for ǫGR = 0; then for emax not to be changed radically when we do include GR, a necessary but insufficient condition is (46). Finding the minimum j at ω = 0 requires that we set the final square bracket in (15) to zero, which is the same as solving the depressed cubic (A7). In §A3 we explained how to determine the appropriate explicit solution to (A7) for arbitrary initial conditions. In particular, we note that the ǫGR = 0 solution jmin = j0 corresponds exactly to equation (A13) with k = 0. However, the general solutions for ǫGR = 0 are not very enlightening. We can make some analytical progress if we further assume that If (B1) is true, then the first order solution for finite ǫGR is In other words, since j 2 0 ∼ Θ ≪ 1 by construction, jmin starts to substantially deviate from j0 when the saddle points appear at ω = 0 -see equation (35) and §3.3.2. Note that the condition (B1) is very stringent and requires that the binary be deep in the very weak GR regime ( §5.1), so (B2) may not be useful in practice.
However, assumption (IV) is new. It is equivalent to the requirement that is the cosine of the binary's minimum inclination. Assumption (IV) is nearly always satisfied at high eccentricities since we normally have 17 j 2 ∼ Θ. The additional assumption (IV) allows us to take the solution for j(t) from (41) that we got using assumptions (I)-(III) and simplify the expression for tmin. The result is: 17 Indeed, equation (24), which is valid in the very weak GR regime, tells us that j 4 f ∼ Θ ≪ 1 and we know that j f then provides an upper bound on j min for Γ > 1/5.
We note that t ′ min diverges as jmin → Θ, that is as emax → e lim . This is as expected from e.g. Figure 6a, since trajectories that approach e lim become ever 'flatter' in the vicinity of ω = ±π/2, i.e. less and less sharply peaked around their eccentricity maxima, so the fraction of a secular period they spend in the vicinity of emax increases.
Next we obtain the solution for ω(t). First, using the conservation of H * (equation (2)) and assumptions (I) and (IV) we easily get an expression for cos 2 ω(j) without stipulating any particular form of j(t): Now plugging in the particular form (C2) for j(t) we find the following explicit solution for 19 ω(t): where and is a dimensionless function of time. In Figure C1a we show how χ varies as a function of jmin/Θ 1/2 ≥ 1. We see that χ > 1 always and that typical values of χ are ∼ a few. Finally we can get the solution for Ω(t) by using assumption (I) in equation (14), plugging in the solutions (C2), (C5) for j(t), ω(t), and integrating in time. The result is where we introduced jz ≡ Jz/L = j cos i. In equation (C8) the value of Ω(0) is an arbitrary constant to be prescribed. Otherwise, equations (C2)-(C3), (C5)-(C8), and the equation jz(t) = jz(0) provide a complete, explicit description of the DA dynamics in the high-e limit whenever assumptions (I)-(IV) are satisfied.
18 To see this we take the explicit expressions for j 1 = j + and j 2 = j 0 from (16)-(19) and simplify them using assumption (I). Plugging the simplified expressions into (41) and expanding the result using assumption (IV) we recover equation (C3). 19 We have included the sgn(t) factor in (C12) because for Γ > 0 the pericentre angle ω must increase towards π/2 as j decreases to j min (at t = 0), and continue to increase as j increases away from j min . We note that ω, Ω make finite 'swings' across the maximum eccentricity peak. Indeed, equation (C5) tells us that ω takes asymptotic values giving a total swing of magnitude |∆ω| = 2 cos −1 [χ −1 − ǫGR/(30Γjmin)]. Clearly the larger ǫGR, the bigger is this swing 20 , which makes sense since GR promotes fast apsidal precession.
Similarly from (C8) we find: Thus, the size of the swing in Ω across the eccentricity peak |∆Ω| = π − O(ǫGR) is reduced by GR effects.

C1 Analytic solution in the LK limit
To apply the analytic solution to the LK case of hierarchical triples, let the tertiary perturber have mass M and the outer orbit have semimajor axis ag and eccentricity eg. Then we set Γ = 1, evaluate ǫGR using equations (3), (6) and take t ′ min equal to (To derive this formula we have used the results of Paper I, Appendix B).
C2 Simplified analytic solution in the limit ǫGR = 0 One can get a simplified version of the analytic solution if one takes the non-GR limit. For ǫGR = 0 the solution for j(t) takes the same form (C2), while (C5) and (C8) simplify to In Figure C1 we show the characteristic behaviour of this non-GR solution. In panel (b) we plot j/jmin as a function of t/t ′ min : obviously j(t) is quadratic in t for t t ′ min and linear for t t ′ min . Panels (c) and (d) demonstrate how the solutions for ω(t) and Ω(t) look for various χ. We note that both angles evolve very rapidly during the interval −1 t/t ′ min 1 and rather slowly otherwise, particularly for χ ≫ 1. We see also that the behaviour of ω is quite strongly dependent on χ; it completes a swing |∆ω| = 2 cos −1 (χ −1/2 ) as t runs from −∞ to +∞. The evolution of Ω depends somewhat less strongly on χ, and its asymptotic value is independent of χ, so that that the total swing in Ω across the eccentricity peak is always |∆Ω| = π. Of course these swings in ω and Ω are not completely correct because we expect our analytic formulae to break down once e differs significantly from unity ( §C3).  Figure C1. Analytic solution to the DA equations of motion without GR at high eccentricity, for binaries that achieve maximum eccentricity at ω = ±π/2. The solution depends on the parameter χ (equation (C6)) which is plotted as a function of j min /Θ 1/2 in panel (a). Panels (b)-(d) show the evolution of j(t), ω(t) and Ω(t) respectively, where we have taken the maximum eccentricity to coincide with t = 0.

C3 Validity of the analytic solution
In this section we test the accuracy of the analytical solution (C2), (C5), (C8) and the simplified solution (C12), (C13) derived in the non-GR limit, against direct numerical integration of the DA equations of motion (12), (13), (14), in different dynamical regimes.

C3.1 Three examples with Γ = 1
First we consider some examples in the LK case of Γ = 1. Precisely, we consider a binary with component masses m1 = m2 = 1M⊙ orbiting a point mass M = 10 5 M⊙. For the outer orbit we choose a pericentre distance rp = 0.4pc and an apocentre distance ra = 0.6pc. The outer orbit is then an ellipse with semimajor axis ag = 0.5pc and eccentricity eg = 0.2. For the inner binary orbit we take the initial conditions e0 = 0.5, i0 = 89.75 • (so that Θ = 1.4 × 10 −5 ), ω0 = 0 • . When we integrate the equations of motion we will shift the time coordinate so that maximum eccentricity is achieved at t = 0; in each example we choose a value of Ω0 so that Ω(0) ≈ π/2. All that remains is to specify the initial semimajor axis a0.
In Figure C2 we Overall, in Figure C2 the analytic solution provides an excellent fit to the exact numerical integration. Errors are only noticeable once e falls below ∼ 0.9 (panels (b) and (d)). This good agreement reflects the fact that σ, j 4 min /Θ ≪ 1 and ǫGR ≪ ǫ weak , meaning that all assumptions (I)-(IV) are fulfilled. Moreover, we see that the full analytic solution (red dashed lines) and non-GR solution (blue dashed lines) overlap almost exactly in panels (c)-(f). This is unsurprising because the binary actually sits in the very weak GR regime ǫGR < ǫ π/2 , meaning GR effects are negligible ( §5.1).
In Figure C3 we use all the same system parameters as in Figure  C2 except we set a0 = 30AU. This increases ǫGR and puts the binary in the weak GR regime ǫ π/2 < ǫGR < ǫ weak . The fact that ǫGR is no longer smaller than ǫ π/2 is responsible for the disagreement between the analytic and non-GR solutions in panels (c)-(f). Nevertheless, the analytic solution still matches the numerical one very well for e 0.9, although not quite as well as in Figure C2, owing to the fact that σ is now comparable to unity (breaking assumption (III)).
Next, in Figure C4 we again run the same experiment but this time with a0 = 15AU. This puts the binary in the moderate GR regime, ǫ weak < ǫGR < ǫstrong, which violates assumption (II). Additionally we have σ ≫ 1, violating assumption (III). We see solutions (C2), (C5), (C8) largely fail to capture the high-eccentricity behaviour even over a very short timescale. At the same time, we note that the moderate GR solution (59) captures the j(t) behaviour extremely well in this case.

C3.2 Three examples with Γ = 0.235
Next we consider some examples with a different value of Γ. To achieve this we replace the Kepler potential with a Hernquist potential Φ(r) = −GM(b + r) −1 , where the total mass M = 10 5 M⊙ and the scale radius b = 1pc. (The outer orbit still has rp = 0.4pc and ra = 0.6pc, but will now fill a 2D annulus rather than forming a closed ellipse -see Paper I). As a result we find Γ = 0.235. This is close enough to Γ = 1/5 that σ and κ attain large values, putting our analytical solutions to a demanding test.
In Figure C5 we integrate exactly the same system as in Figure  C2 except for this replacement of the potential -in particular, we . Comparing analytic results at high eccentricity with exact numerical integration (see §C3 for details). In this example ǫ GR < ǫ π/2 so the binary is in the very weak GR regime. Note also that σ ≪ 1.  Figure C2 except we take a 0 = 30AU, so the binary is in the weak GR regime, ǫ π/2 < ǫ GR < ǫ weak , and the value of σ is now approaching unity.  Figure C4. As in Figure C2 except we take a 0 = 15AU, so the binary is in the moderate GR regime, ǫ weak < ǫ GR < ǫstrong. Note that σ is significantly larger than unity in this case.   Figure C6. As in Figure C5 except we take a 0 = 40AU. Again the binary is in the weak GR regime, ǫ π/2 < ǫ GR < ǫ weak , but is approaching the moderate GR regime.   Figure C7. As in Figure C5 except we take a 0 = 35AU. In this case the binary is (just) in the moderate GR regime ǫ weak < ǫ GR < ǫstrong.
again take a0 = 50AU. We see that this puts the binary in the weak GR regime ǫ π/2 < ǫGR < ǫ weak (as in Figure C3), but owing to the proximity of Γ to 1/5 the value of σ is much larger than unity (unlike in Figure C3). One consequence of this is that the analytic approximation to log 10 (1 − e) fails rather early on, with significant errors by the time e falls below 0.99 ( Figure C5b). Despite this, the analytic approximations to ω(t) and Ω(t) are still excellent (panels (c)-(f)). This is because ω and Ω are sensitive only to the eccentricity behaviour at the very peak -they change very rapidly over the interval −2.2 < t/t ′ min < 2.2 (panels (c) and (e)), but are almost constant the rest of the time. Thus, as long as j(t) is captured well near the very peak eccentricity, as it is in panel (a), the analytic solutions for ω, Ω work well despite assumption (II) being broken.
In Figure C6 we investigate the same system except with a0 reduced to 40AU. The binary is still in the weak GR regime but only just so, violating assumption (II). It also has σ ≫ 1 like it did in Figure C5, violating assumption (III). We see that the analytic fit to log 10 (1 − e) is quite poor even at the very peak (panel (a)). Interestingly though, the evolution of ω (panel (d)) is reproduced rather accurately, highlighting how sensitive ω(t) is to the value of peak eccentricity jmin (see equation (C9)), and how insensitive it is to anything else. However, the evolution of Ω(t) is not reproduced very well. The same conclusions hold for Figure C7, in which we have reduced the semimajor axis further to a0 = 35AU, putting the binary squarely in the moderate GR regime (so that both assumptions (II) and (III) are broken).

C3.3 Conclusions
While assumptions (I) and (IV) are almost always good provided we consider binaries that reach very high eccentricity (1 − e ≪ 0.1), assumptions (II) and (III) are liable to fail in some regimes.
We have seen that for log 10 (1 − e) to be accurately reproduced by the analytical solution (C2) for e 0.9, all four assumptions (I)-(IV) must be valid.
However, the analytic solution (C8) for Ω can be very accurate even for σ ≫ 1 (violating assumption (III)) provided the behaviour of j(t) in the close vicinity of jmin is reproduced reasonably well.
What is more, the solution (C5) for ω(t), and the swing ∆ω in particular, can be very accurate even if the system is in the moderate GR regime, invalidating both assumptions (III) and (IV). This is because ω is extremely sensitive to the behaviour of j around absolute peak eccentricity and largely insensitive to j otherwise.
D1 Phase space behaviour in the case −1/5 < Γ ≤ 0 Unlike for Γ > 0, the dynamical behaviour in the regime −1/5 < Γ ≤ 0 cannot be understood using only one value of Γ as an example. Thus, we consider three values. In Figures D1, D2 and D3 we plot contours of constant H * in the (ω, e) phase space for Θ = 0.1, taking Γ = −0.1, Γ = −0.15 and Γ = −0.19 respectively. The manually-added dashed contours are the same as in Figure 2. We now discuss these three figures in turn.
First we discuss Figure D1 (Γ = −0.1). From Paper II we know that when ǫGR = 0, fixed points never exist in the phase space for −1/5 < Γ ≤ 0. Thus all phase space trajectories circulate and their maximum eccentricity is found at ω = ±π/2, as in panel (a). Now we consider finite ǫGR. In panel (b), namely for ǫGR = 1.0, we see that fixed points have appeared at ω = 0, ±π, which we will refer to simply as ω = 0 from now on. These fixed points are not saddle points like they were for 0 < Γ ≤ 1/5 ( Figure 2); instead they are maxima of H * and host a region of librating orbits that is connected to e lim .
As we increase ǫGR further we see that these fixed points move down the page to lower eccentricity, and their associated librating islands become larger in area. At some threshold value of ǫGR the librating islands become disconnected from the line e = e lim , coinciding with the appearance of new saddle points at ω = ±π/2, e = e lim . As we increase ǫGR beyond this threshold the fixed point at ω = 0 continues to move down the page (panels (c) and (d)), as do the saddle points at ω = ±π/2, and a new family of high-e circulating orbits runs over the top of the librating islands, reminiscent of what we saw for 0 < Γ ≤ 1/5 in Figure 2. Partitioning the different librating islands and circulating regions in panels (c) and (d) are separatrices that cross at the saddle points. Continuing to increase ǫGR forces both kinds of fixed point to move to lower eccentricities. The saddle points move fastest and disappear first; in panel (e), the fixed point at ω = 0 remains but the saddle points at ω = ±π/2 have disappeared through e = 0. Increasing ǫGR even further still, the ω = 0 fixed point reaches e = 0 and then disappears. This leaves a phase space filled with circulating trajectories (panel (f)), which is similar to the ǫGR = 0 case shown in panel (a) except that the maximum eccentricities are now found at ω = 0 rather than ω = ±π/2, and the locations of the maxima and minima of H * are reversed (see the colourbars).
Moving on to Figure D2 (Γ = −0.15), we find a completely different picture of rather impressive dynamical diversity. In this figure we have to use twelve panels to fully illustrate the complex phase space behaviour. To begin with, panels (a) and (b) in Figure D2 have the same morphology as Figures D1a,b. However, panel (c) is very different from Figure D1c. This time, at some threshold value of ǫGR a pair of fixed points emerges from a single point at ω = π/2, e = e f,π/2 , and the same thing happens at ω = −π/2. An increase in ǫGR nudges these fixed points apart in their eccentricity values (panel (d)): one of them moves up the page and the other moves down. In each pair, the fixed point with higher e is a minimum of H * and hosts a region of librating orbits. The fixed point with lower e is a saddle point, and sits on the separatrix that surrounds the upper point's librating region. In addition we still have the usual fixed point and accompanying librating island at ω = 0. As a result, we now find two families of circulating trajectories. One runs close to e = 0 under the separatrices passing through the saddle points. The other runs above these separatrices, but below the separatrices surrounding the librating islands centred on ω = 0, ±π. This second type of circulating trajectory reaches high eccentricity by running above the upper fixed points at ω = ±π/2. Quite remarkably, these circulating trajectories also exhibit non-monotonic behaviour of ω(t), i.e.ω is > 0 at some times and < 0 at others, despite the trajectory being a circulating one.
Increasing ǫGR further, the upper fixed point (minimum) at ω = ±π/2 continues to move up the page, while the lower fixed point (saddle) moves down (panels (e)-(g)). Meanwhile the ω = 0 fixed points also move down the page, albeit much more slowly. Eventually the librating region surrounding the upper fixed point at ω = ±π/2 becomes connected to e = e lim . Simultaneously, the saddle point at ω = ±π/2 and its associated separatrices merge with the separatrices surrounding the ω = 0 librating regions (see the transition from panel (g) to panel (h)). Accompanying this transition is the change in the nature of the second family of finite eccentricity circulating orbits described above -they now run above (below) the saddle points at ω = 0, ±π (ω = ±π/2). As ǫGR continues to increase the pair of fixed points at ω = ±π/2 continue to move apart in eccentricity, until eventually the lower one disappears at e = 0 (panel (j)) followed by the upper one at e = e lim (panel (k)). In panels (k) and (i) we retain only the fixed points at ω = 0, with qualitatively the same overall phase space behaviour as in Figure D1e. The ω = 0 fixed points also disappear once ǫGR becomes sufficiently large. Figure D3 (Γ = −0.19) shows yet again a different qualitative behaviour. Like in Figures D1, D2, fixed points emerge at ω = 0 followed by additional fixed points at ω = ±π/2, e = e f,π/2 (panel(b)). However, this time the ω = ±π/2 fixed points do not come in pairs like they did in Figure D2. Instead they are minima of H * and are surrounded by a librating island that stretches to e = 0 (though e f,π/2 = 0 for any ǫGR). These fixed points move up the page as we increase ǫGR (panel(c)) until they become connected to e = e lim (panel (d)). At this stage circulating trajectories exhibit a transition similar to that in D2. Thereafter we have qualitatively the same behaviour as in Figure D2j.
As these three examples demonstrate, the qualitative dynamical behaviour in the regime −1/5 < Γ ≤ 0 is highly complex. It is also very difficult to analyse mathematically. The simplest place to start is with the fixed points at ω = 0, j = j f,0 . The formulae describing these fixed points can be carried over from §3.2 and Appendix A: the value of j f,0 is still determined by equation (34) and the fixed points exist provided equation (35) is true. The key difference for negative Γ compared to positive Γ is that the determinant of the Hessian matrix of H * evaluated at the fixed points, namely the expression (A6), is now manifestly positive rather than negative. Thus the fixed points at ω = 0, j = j f,0 are now true extrema (more precisely, maxima) of H * and host a librating island, which is reflected in Figures D1-D3.
In other words, increasing Θ decreases the eccentricity of the fixed points at ω = ±π/2 should they exist (same as Γ > 0), but increasing ǫGR increases their eccentricity (opposite to Γ > 0). The condition on ǫGR for the existence of these fixed points is the same as (29) but reversing the inequalities, i.e. replacing each '<' with '>'. The condition on Θ is the same as that given for 0 < Γ ≤ 1/5 in equation (30). Additionally, the fixed points at ω = ±π/2 are always true extrema (minima) of H * in this Γ regime since the expression (D1) is always positive. Meanwhile, the fixed points at ω = 0 follow exactly the same rules as for −1/5 < Γ ≤ 0, appearing at ω = 0, e = e lim when ǫGR reaches the critical value ǫGR = 6(1 − 5Γ)Θ 3/2 and then working their way down the page towards e = 0 as ǫGR is increased, disappearing for ǫGR > 6(1 − 5Γ) (equation (35)). The only difference is that these fixed points are maxima of H * , not saddle points, which follows from the fact that the quantity (A6) is positive for Γ < 0.
D3 Orbit families and maximum eccentricity for Γ ≤ 0 regimes Owing to the highly complex phase space morphology, working out a trajectory's orbital family analytically is often a very tedious job for negative Γ values. The same goes for finding a binary's maximum eccentricity: in practice it is best simply to take a brute-force approach by solving the cubic and quartic equations (A7), (37) numerically to get all seven possible roots, and then declaring jmin to be the real root closest to but smaller than the initial j value. We will not pursue any further technical details here.

D3.1 Maximum eccentricity for initially near-circular binaries
With this brute-force approach it is straightforward to calculate emax for a given i0, Γ and ǫGR for initially near-circular binaries when Γ ≤ 0 (c.f. §3.4.1). In Figure D5 we show emax(i0) for various ǫGR values. In each panel we use a different negative value of Γ (c.f. Figure 5).
In the top row of Figure D5 (panels (a)-(c)) we explore the regime −1/5 < Γ ≤ 0. To understand panels (a) and (b) it is worth looking back at Figures D1 and D2 and asking what we expect of the behaviour of initially near-circular orbits. We expect from Figures  D1a,b,c,d and D2a-i that for low enough ǫGR the maximum eccentricity will be zero. This immediately explains, for instance, why there is no red line (corresponding to ǫGR = 0) in panels (a) and (b) of Figure D5. However, when ǫGR takes a value such that (I) the fixed point exists at ω = 0 and (II) the librating region that this fixed point hosts is connected to e = 0, then the eccentricity of initially circular orbits is maximised at ω = 0 and is nonzero ( Figure D1e and Figure D2i,j,k).
We know that (I) is true if and only if ǫGR satisfies (35). We also know that for (II) to be true the fixed points at ω = ±π/2 must have disappeared below e = 0, which happens at the critical Θ value Θ = Θ1(Γ, ǫGR) given in equation (31) for arbitrary Γ. Putting these constraints together, taking Θ → 0 and using the fact that −1/5 < Γ ≤ 0, we find that a necessary condition for both (I) and (II) to be true is 6(1 + 5Γ) < ǫGR < 6(1 − 5Γ), which is the same as (36) if we replace '<' with '>'. For Γ = −0.05 this gives 4.5 < ǫGR < 7.5, which is why there is only a cyan line in Figure  D5a. Similarly, for Γ = −0.1 we get 3 < ǫGR < 9, hence the solo cyan line in Figure D5b.
This necessary constraint on ǫGR was derived for Θ → 0, i.e. i0 → 90 • . To find the necessary constraint on i0 for both (I) and (II) to be true, we set Θ = cos 2 i0 = Θ1(Γ, ǫGR). Then initially near-circular binaries whose Γ values produce phase portraits like in Figures D1, D2 can achieve a finite emax only if they have i0 greater than the critical value cos −1 Θ1(Γ, ǫGR). For panels (a) and (b) of Figure D5 these values are i0 = 66 • and i0 = 55 • respectively, which we show with vertical dotted lines. Plugging H * , Θ from (38) into the depressed cubic (A7), we find that the corresponding minimum j value is: which is independent of i0. In panels (a) and (b) of Figure D5, the straight horizontal cyan lines for ǫGR = 5 correspond to the solution (D3). Panels (c)-(f) of Figure D5 all share a similar morphology, so we will consider them together. In each panel, for a fixed ǫGR a finite emax arises at some critical value of i0, increases as a function of i0 until it reaches emax = e lim , and then is constant for all larger values of i0 up to 90 • . Note that on the non-constant parts of these curves we have essentially the opposite of the intuitive Γ > 1/5 result: for a fixed initial inclination, a larger ǫGR leads to a larger emax. The behaviour we see in these panels is consistent with what we expect from the Γ = −0.19 example given in Figure D3 and the Γ ≤ −1/5 example we studied in Figure D4. In those figures, the fixed points that emerge at ω = ±π/2, e = e f,π/2 host regions of librating orbits that are connected to e = 0. In each case, the maximum eccentricity of initially circular orbits is determined by eccentricity of the separatrix at the point ω = ±π/2. As we increase ǫGR, the value of e f,π/2 is increased, pushing the separatrix to higher e, and so the maximum eccentricity of initially circular orbits grows. Eventually, however, e f,π/2 is increased so much that the separatrix reaches e = e lim (dashed black line) -see the transition between Figure D4c and D4d. At the same time, the librating islands that are hosted by fixed points at ω = 0 become connected to e = 0. After that the maximum eccentricity is given by (D3) and is independent of i0 -hence the straight horizontal lines in Figure D5c-f. The main qualitative difference between panel (c) and panels (d)-(f) is that panel (d) exhibits no red line, i.e. no solution for ǫGR = 0. This is because in the regime −1/5 < Γ ≤ 0 a finite ǫGR is always required for any fixed points to exist ( Figures D1-D3).
Finally we may briefly compare Figure D5 with Figure 5. Consider what happens if we fix ǫGR and increase i0 from zero. In Figure D5 (Γ ≤ 0), the larger is ǫGR, the lower i0 is required to achieve a non-zero emax (provided ǫGR is not so large that no eccentricity excitation is possible). On the contrary, in Figure 5 (Γ > 0) the most favourable situation for eccentricity excitation is always to have ǫGR as small as possible: the larger ǫGR, the larger i0 is required to get a non-zero maximum eccentricity. Though the two regimes differ in this respect, they are similar in that for binaries with i0 ≈ 90 • a larger ǫGR always leads to a smaller maximum eccentricity (again, provided ǫGR is such that eccentricity excitation is possible). This paper has been typeset from a T E X/L A T E X file prepared by the author.