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

ABSTRACT

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 generalize 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.

1 INTRODUCTION

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 2019, 2020). These discoveries certainly warrant an astrophysical explanation, which is complicated by the fact that the time-scale 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 |$m_1=m_2=30\, {\rm M}_{\odot }$| will only merge within 10¹⁰ yr if its initial semimajor axis a is ≲0.2 au. 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¹ 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 (Kozai 1962; Lidov 1962), or simply the mean field potential of the star cluster or galaxy in which the binary resides (Heisler & Tremaine 1986; Brasser, Duncan & Levison 2006; Bub & Petrovich 2019; Hamilton & Rafikov 2019a,b).

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 primarily 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 those arising from rotational or tidal bulges; see e.g. Liu, Muñoz & Lai 2015; Muñoz, Lai & Liu 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, Murray & Mikkola 2014; Rodriguez et al. 2015; Liu et al. 2015; Hamers et al. 2018; Liu & Lai 2018; Samsing, Hamers & Tyles 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 (Blaes, Lee & Socrates 2002; Miller & Hamilton 2002; Wen 2003; Veras & Ford 2010; Liu et al. 2015; Anderson, Lai & Storch 2017; Grishin, Perets & Fragione 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; Bub & Petrovich 2019; Hamilton & Rafikov 2019c) where we expect similar considerations to apply.

The main purpose of this paper is to explore systematically the effect of GR precession on the underlying phase space dynamics and eccentricity evolution of a tidally perturbed binary. We will focus exclusively on the ‘doubly averaged’ (DA)² 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 (Naoz 2016). 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 Section 5.3. The DA approximation will be relaxed in an upcoming paper (Hamilton & Rafikov, in preparation). This 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 Section 2, we write down the DA perturbing Hamiltonian and establish the notation that we will use for the rest of the paper. In particular, we introduce the key parameter ϵ_GR that measures the strength of GR precession relative to tidal torques. In Section 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 Section 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 time-scale of high eccentricity episodes. In Section 5, we discuss our results in the light of the existing literature, and comment on the limitations of our study. We summarize in Section 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.

2 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.

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.

2.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 mid-plane of a thin disc (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ζ²)⁻¹, where ζ ≡ b/R [this special case was considered by Brasser et al. (2006); compare their disturbing function (A.5) with our equation (4)].

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.

3 PHASE SPACE BEHAVIOUR

To gain a qualitative understanding of the dynamics driven by the Hamiltonian equations (12)–(13), one can fix the values of Γ, Θ, and ϵ_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’.

In Section 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²)cos ²i is conserved, a necessary condition for this is that Θ ≪ 1. Fixed points are crucial to this investigation because they may force an initially low-e binary to very high maximum eccentricity e_max. 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 e_max. 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’ e_max 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 (Section 3.1) and for 0 < Γ ≤ 1/5 (Section 3.2). We describe some properties of the fixed points in Section 3.3, and then show how to calculate the maximum eccentricity of a given binary in Section 3.4. Details of the mathematical results that we quote throughout Sections 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

Fig. 1 shows phase portraits for the case Γ = 0.5 > 1/5. In the top (bottom) row, we set |$\Theta = 0.1\, (0.5)$|⁠. From left to right, we vary ϵ_GR taking ϵ_GR = 0, 2, 5, 10, and 30. In each panel, a black horizontal dashed line shows the limiting possible eccentricity |$e_\mathrm{lim} = \sqrt{1-\Theta }$| (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 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; 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 Fig. 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 Section 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 Section 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 Fig. 2, we show phase portraits similar to Fig. 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 Fig. 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 (Figs 2a and f) and Γ > 1/5 (Figs 1a and 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 host 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 that we determine in Section 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). 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_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.

Physically, these new features might have been anticipated. First of all, in Paper II we saw that for 0 < Γ ≤ 1/5, the cluster-driven ω evolution of circulating trajectories is always retrograde (contrary to the Γ > 1/5 case in which it is prograde). Since GR always promotes prograde precession, its effect in this Γ regime is initially to slow down the overall precession rate |$\dot{\omega }$|⁠, allowing for a more coherent torque compared to the case of ϵ_GR = 0. This leads to the appearance of new librating solutions. This is first true for binaries that previously lie just below the separatrix in the (ω, e) phase space, as these circulate the slowest (recall that the secular period diverges on the separatrix itself), giving |$\dot{\omega } \sim 0$| for a relatively small value of ϵ_GR. On the other hand, at the highest binary eccentricities (near e = e_lim) GR may dominate the dynamics, causing the binary’s pericentre angle ω to precess rapidly, leading to the appearance of new high-e circulating solutions.

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 and π. 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 Section 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 Figs 2(c) and (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 Fig. 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 (Figs 2e and j). The amplitude of eccentricity oscillations decreases correspondingly until cluster tides are completely negligible and only GR apsidal precession remains.

3.3 Fixed points

We now proceed to understand mathematically the nature of the various fixed points that we found in the phase portraits in Sections 3.1–3.2. By setting dj/dt = 0 in equation (13), we see that all possible non-trivial fixed points⁶ are located on (i) the lines ω = ±π/2, as in Paper II, and/or (ii) the lines ω = 0, ±π, consistent with Figs 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 points.

3.3.1 Fixed points at ω = ±π/2

In Section 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 Figs 1 and 2. The values of j_f,π/2 are found as solutions to the quartic polynomial (A1), and we illustrate their behaviour in Fig. 3 for several values of Γ and Θ. One can see that j_f,π/2 always increases with ϵ_GR (see A4), explaining why in Figs 1 and 2 the fixed points at ω = ±π/2 always get pushed to lower e as ϵ_GR is gradually increased from zero.

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 Section 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).

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 Fig. 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 realized 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 Fig. 2, respectively), a fixed point appears at the limiting eccentricity |$e_\mathrm{lim} = \sqrt{1-\Theta }$|⁠. Increasing ϵ_GR at fixed Γ always acts to increase j_f,0, i.e. to decrease the eccentricity |$e_\mathrm{f,0} \equiv (1-j_\mathrm{f,0}^2)^{1/2}$| of this particular fixed point. As we increase ϵ_GR to the threshold value ϵ_GR = 6(1 − 5Γ) (which is independent of Θ and corresponds to ϵ_GR = 3 in Fig. 2), the saddle point vanishes through e = 0, leaving only the fixed points at ω = ±π/2.

3.4 Determination of the maximum eccentricity of a given orbit

Our next goal is to calculate the maximum eccentricity e_max reached by a binary given the initial conditions (ω₀, e₀, Θ, Γ, ϵ_GR). In particular, we wish to know if a binary will reach e_max → 1, since this is the regime in which dissipative effects (e.g. GW emission) can become important.

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 Figs 2b, c, and h). To do so, we use the procedure given in Section 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’ the librating regions and the saddle point, then we have j_min = j(ω = 0). Otherwise, j_min is found at ω = ±π/2 and we proceed as for Γ > 1/5 by solving equation (37).

3.4.1 Maximum eccentricity achieved by initially near-circular binaries

Meanwhile, an orbit that has e₀ infinitesimally larger than zero can have j_min 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 Fig. 1). Because of the constraint (29), a necessary (and for i₀ → 90°, sufficient) requirement for this is ϵ_GR < 6(1 + 5Γ). In that case, the fixed points bound the maximum eccentricity from below, so |$e_\mathrm{max}\gt e_{\mathrm{f},\pi /2} \equiv (1-j^2_{\mathrm{f},\pi /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 e_max = 0 (see panels e, i, and j of Fig. 1).

Next, we turn to the regime 0 < Γ ≤ 1/5. By consulting Fig. 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 Figs 2d and i). A necessary requirement for this (which is again sufficient in the case i₀ → 90°) is that (36) be true. Then, e is maximized at ω = ±π/2 and j_min is a non-trivial solution to equation (39). On the other hand, if (36) is not satisfied then a binary that starts at e₀ ≈ 0 never increases its eccentricity⁸ even for i₀ = 90°.

Overall then, we see that for near-circular binaries to reach finite e_max we require fixed points to exist in the phase portrait at ω = ±π/2 but not at ω = 0, and this necessarily requires ϵ_GR to satisfy (36).

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 (e_max = 0).

A related phenomenon is that in panel (e), the green (ϵ_GR = 3) curve asymptotes to the black dashed line e = e_lim as i₀ → 0°. This is also as expected: Since Γ = 0.1, equation (36) tells us that ϵ_GR = 3 is precisely the lower bound on GR strength above which initially near-circular binaries can reach non-zero eccentricities at all i₀, since at this value of ϵ_GR the saddle point crosses e = 0 – see Figs 6(i)–(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 i₀ ≈ 0° inevitably leads to e_max ≈ 0 – in other words, there is no eccentricity excitation for initially coplanar (i₀ = 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 e_max 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 Fig. 5 reinforces the idea that the 0 < Γ ≤ 1/5 regime admits far fewer high-eccentricity solutions than Γ > 1/5 as ϵ_GR is varied.

4 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 example, 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 e_max but also the behaviour of e(t) and other orbital elements in the vicinity of e_max.

In Section 3.4, we explained how to find e_max for arbitrary Γ > 0, initial conditions (e₀, i₀, ω₀), and value of ϵ_GR. Here, we will examine the solutions quantitatively in the high eccentricity limit, and 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)–(A8) that determines the roots at ω = 0, i.e. what we called j(ω = 0).

In this section, we will focus on situations in which e_max is achieved at ω = ±π/2, since this is the most common prerequisite for e → 1 (Sections 3.1–3.2). The rare cases in which e approaches unity at ω = 0 are covered in Appendix  B.

4.1 Phase space behaviour for |$\Theta \ll 1, \, \, \, \Gamma \gt 0$|

We are interested in binaries that start with initial eccentricity e₀ not close to unity, and that are capable of reaching extremely high eccentricities e_max → 1, i.e. j_min → 0. For this to be possible, a necessary condition is that Θ ≪ 1, owing to the constraint (9). Hence, it is important to understand the regime Θ ≪ 1 in detail.

In Fig. 6, we show phase portraits for Γ = 0.5 (top row) and Γ = 0.1 (bottom row), this time fixing Θ = 10⁻³ in both cases, and adding in extra dashed contours⁹ 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 Figs 1 and 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), since in that case ϵ_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⁻³ is approximately 10⁻⁵ < ϵ_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 Fig. 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 near-circular 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). More precisely, for 0 < Γ ≤ 1/5, the eccentricity of the saddle point (34) acts as a hard boundary on the maximum eccentricity of low-e orbits, and most of them do not get close even to that value. Even when ϵ_GR is increased so that a new family of circulating orbits appears, and the librating region is significantly enlarged, the system admits very few solutions that start at low e and achieve high e. 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).

4.2 High eccentricity behaviour for ϵ_GR = 0

Note that the solution (42) is quadratic in t for t ≲ t_min and linear when t ≳ t_min, as long as j remains ≪1. It provides a 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 (Section C3).

4.3 Modifications brought about by finite ϵ_GR

Before we proceed to examine the j(t) behaviour, it is important to realize that including a finite ϵ_GR affects the right-hand side of (15), and therefore the value of j_min, in two distinct ways. First, there is the obvious explicit dependence on ϵ_GR that appears twice in equation (15). Secondly, there is also an implicit dependence on ϵ_GR in (15) through the values of j_± and j₀ (see equations 16 and 17). We will now discuss this implicit dependence, and then use the results to understand j(t) behaviour in different asymptotic ϵ_GR regimes.

4.4 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¹¹ ω = ±π/2. As shown in Paper II, for Γ > 0 the corresponding minimum angular momentum for these orbits is always j_min = j₋ ≪ 1.

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 maximized at ω = ±π/2).

4.4.1 Weak GR, ϵ_GR ≪ ϵ_weak

In the case σ ≪ 1, the term σj₋ in the final square bracket in (55) can also be dropped compared to j. Then, equation (55) takes the same functional form as its non-GR analogue; integrating, we get a solution j(t) in precisely the form (42) with¹⁵j_min → j₋ and j₁, j₂ → j₊, j₀. This is reflected in Fig. 7(a), in which the black line (σ = 0) is exactly the non-GR result from (42), and as expected |$j=\sqrt{2}j_\mathrm{min}$| coincides with t = t_min in that case.

In Figs C2, C3, C5, and C6, we compare the weak GR solution for j(t), namely equation (56), to direct numerical integration of the DA equations of motion (12) and (13), for binaries in different dynamical regimes. Full details are given in Section 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₁₀(1 − e) behaviour in the vicinity of peak eccentricity, while panel (b) shows the same thing zoomed out over a much longer time interval.¹⁶ The weak GR solution for j(t) [equation (56)] 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 (Figs C2 and C3) this weak GR solution works very well, but that substantial errors begin to set in when ϵ_GR approaches ϵ_weak (Figs C5 and C6). Finally, in each of these plots we also show with red dashed lines an ‘analytical’ solution, equation (C2), which coincides with (42) provided |$j_\mathrm{min}^4/\Theta \ll 1$|⁠. As we have already stated, in the weak GR regime j(t) takes the form (42) provided that σ ≪ 1, so it is unsurprising that in the plot with very small σ (Fig. C2) this analytical solution (equation C2, denoted with red dashed lines) overlaps with the weak GR solution (equation 56, shown with green dashed lines).

4.4.2 Moderate GR, ϵ_weak ≪ ϵ_GR ≪ ϵ_strong

Finally, in panels (a) and (b) of Figs C4 and C7 we compare the moderate GR solution (62), 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.

4.5 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 eventually brings any merging binary to a small enough semimajor axis to put it in this regime.

Interestingly, the minimum angular momentum j_min = j(0) predicted by the solution (67) can still be described by the expression (52) in the limit ϵ_GR ≫ ϵ_strong. To see this, we note that for ϵ_GR ≫ ϵ_strong equation (52) gives |$j_{\min }\approx \epsilon _\mathrm{GR}/(j_+^2\epsilon _\mathrm{strong})$|⁠. Then, we take the expression (46) for |$j_+^2$|⁠, and substitute into it the value of Σ we get by taking ϵ_GR ≫ 1 in (18)–(19), namely |$\Sigma \approx \epsilon _\mathrm{GR}(1-e_0^2)^{-1/2}/6$|⁠. Putting these pieces together, we find |$j_\mathrm{min}\approx \sqrt{1-e_0^2}=j(0)$|⁠. Thus, the solution (52) interpolates smoothly between the different asymptotic GR regimes.

4.6 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 J_z(t) = J_z(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¹⁷ (42) 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 (42) is a good approximation to the j(t) behaviour near the peak eccentricity, as we verify numerically in Section 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 for instance in order to explore the short-time-scale (i.e. non-DA) effects near peak eccentricity, which are important for accurate calculation of the LK-driven merger rate (Grishin et al. 2018).

5 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 characterized 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 uncovered a completely new pattern of secular evolution, which we characterized 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 Section 2.1, the Γ ≤ 0 regime splits into two further regimes, namely −1/5 < Γ ≤ 0 and Γ ≤ −1/5. We found that in both of these regimes the resulting phase space structures and maximum eccentricity behaviour are considerably more complex and counter-intuitive than those for Γ > 0.

Furthermore, we have explored the evolution of binary orbital elements in the limit of very high eccentricity (Section 4). This investigation revealed a number of distinct dynamical regimes that are classified according to the value of ϵ_GR. In Section 5.1, we summarize and systematize these regimes based on their physical characteristics. In Section 5.2, we compare our study to the existing LK literature, and in Section 5.3 we discuss its limitations.

5.1 Summary of ϵ_GR regimes and their physical interpretation

These characteristic values of ϵ_GR allow us naturally to delineate four important regimes of secular dynamics:

• very weak GR: ϵ_GR ≲ ϵ_π/2,

• weak GR: ϵ_π/2 ≲ ϵ_GR ≲ ϵ_weak,

• moderate GR: ϵ_weak ≲ ϵ_GR ≲ ϵ_strong,

• strong GR: ϵ_GR ≳ ϵ_strong.

Based on our findings in Sections 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 Fig. 9), GR precession has essentially no effect on the dynamics for Γ > 1/5. More precisely, 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 (54). 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 Section 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 Fig. 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 j_min (or equivalently e_max), which stays close to its ϵ_GR = 0 value j₋ (see equation 54). If σ ≫ 1, then GR also modifies the time spent in the high-eccentricity state (equation 59).

Moderate GR:  In this regime (between the dotted and solid curves in Fig. 9), GR precession modifies not only the locations of the fixed points but also the values of j_min (and hence of e_max), now given in equation (60). GR also modifies the time spent in the high-e state (equation 66). 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 Fig. 9), GR precession dominates the binary dynamics at all times; the quantities j_± are significantly affected by GR precession (equation 44) and all fixed points in the phase portrait disappear (equations 33 and 35). Cluster tides drive only very small eccentricity oscillations on top of uniform GR precession, so that e is roughly constant – see equation (67).

In Table 2, we summarize the main features of the asymptotic ϵ_GR regimes that we have found in this and previous sections.

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 ∼ t_sec and j ∼ 1, we obtain ϵ_GR ∼ 1, which is consistent with the definition (32) of ϵ_strong up to a numerical coefficient.

Secondly, 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 time-scale of secular evolution near the eccentricity peak. In Section 4, we find quite generally this time-scale to be t_ch ∼ t_min ∼ j_mint_sec – see e.g. equations (42) and (66). Plugging this into the condition (69) and evaluating |$\dot{\omega }_\mathrm{GR}$| at j_min, we immediately find that j_min ∼ ϵ_GR, in agreement with equation (60). When the GR barrier first emerges at the transition between weak and moderate regimes, j_min is still well approximated by the ϵ_GR = 0 solution j₋ ∼ Θ^1/2 (see equation 46). As a result, the ϵ_GR value corresponding to this transition is ∼Θ^1/2, in agreement with the definition (51) of ϵ_weak.

Thirdly, we expect fixed points in the phase portrait at ω = ±π/2 to be substantially displaced by GR precession when |$\dot{\omega }_\mathrm{GR}(j_\mathrm{f})$| becomes comparable to the characteristic secular frequency |$\dot{\omega }$| of libration around a fixed point. Since we are interested in the displacement of j_f by an amount ∼j_f, we take this |$\dot{\omega }$| from the H^⋆ = const. contour centred on the ω = ±π/2 fixed point and with vertical extent ∼j_f. Plugging j ∼ j_f into equation (12) and using expression (24) for j_f, we find |$\dot{\omega }\sim \Theta ^{1/4}t_{\rm sec}^{-1}$|⁠, so that in this case t_ch ∼ Θ^−1/4t_sec. Substituting this into the condition (69) 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 Sections 3–4 to obtain the numerical coefficients, which are actually quite important. Indeed, equations (26), (51), and (32) feature constant numerical factors that can substantially exceed unity, especially for the LK case of Γ = 1.

5.2 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 Section 3.4.1) was derived by Miller & Hamilton (2002), Blaes et al. (2002), Wen (2003), Fabrycky & Tremaine (2007), and Liu et al. (2015). Of these, Fabrycky & Tremaine (2007) and Liu et al. (2015) also produced plots very similar to Fig. 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 (52) 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 Section 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 Section 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 section 2.1 of Rodriguez & Antonini (2018), there are at least two ways to decide when GR dominates. One method is to take j_min corresponding to the pericentre distance that needs to be reached according to the problem at hand, and then equate |$\dot{\omega }_\mathrm{GR}(j_\mathrm{min})$| with the precession rate due to the tidal perturbations (see their equation 29), which corresponds to equation (54) in Paper II. As discussed in Section 5.1, this method would set a rough upper limit of ϵ_GR ≲ j_min; see equation (60) 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 not a very stringent requirement, and does not guarantee that the majority of systems with such ϵ_GR would reach the required j_min – 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 Fig. 9).

With regard to phase space structure, the only study we know of that resembles our Section 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 section 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 fig. 2 is qualitatively the same as our Fig. 1.

5.3 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 octupole-order effects are important, the maximum eccentricity can actually be increased by GR precession (Ford, Kozinsky & Rasio 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 ∼ 10 au orbiting a stellar cluster at ∼1 pc, 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-time-scale 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 (Section 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).

6 SUMMARY

In this paper, we completed our investigation of DA (test-particle quadrupole) cluster tide-driven binary dynamics in the presence of 1PN general-relativistic pericentre precession. Throughout, we parametrized the strength of GR precession relative to tides using the dimensionless number ϵ_GR (equation 6). We can summarize 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 and ±π.

• 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 j_min 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 analytical 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-time-scale fluctuations (‘singly averaged effects’) on binaries undergoing cluster tide-driven secular evolution, as well as the population synthesis calculations of merger rates.

ACKNOWLEDGEMENTS

We thank the anonymous referee for several insightful comments on the manuscript. CH is funded by a Science and Technology Facilities Council (STFC) studentship. RRR acknowledges financial support through the STFC grant ST/T00049X/1, NASA grant 15-XRP15-2-0139, and John N. Bahcall Fellowship.

DATA AVAILABILITY

No new data were generated or analysed in support of this research.

Footnotes

REFERENCES

APPENDIX A: MATHEMATICAL DETAILS OF PHASE SPACE BEHAVIOUR FOR Γ > 0

In this appendix, we provide some mathematical details for the results quoted in Section 3.

A1 Fixed points at ω = ±π/2

The criteria for these ω = ±π/2 fixed points to exist can be found by demanding that the condition (9) is obeyed, i.e. that |$\sqrt{\Theta } \lt j_\mathrm{f,\pi /2} \lt 1$|⁠. Let us begin by fixing Γ and Θ; then, owing to the monotonic behaviour of j_f,π/2(ϵ_GR) (equation A4), we simply look for the ϵ_GR values that correspond to |$j_{\mathrm{f},\pi /2}=\sqrt{\Theta }$| and j_f,π/2 = 1. Doing so, we arrive straightforwardly at the condition (29) on ϵ_GR.

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) and (30). We can use this information to understand Fig. 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.

For Γ > 1/5, the lines in Fig. 4 correspond to Θ_max = Θ₁ (equation 31). Then ∂Θ₁/∂Γ = Γ⁻²(ϵ_GR − 6)/60, so that Θ₁ increases (decreases) monotonically with Γ for ϵ_GR > 6 (ϵ_GR < 6). For the special value ϵ_GR = 6, we have Θ₁ = 1/2 = const., hence the straight horizontal brown line in Fig. 4.

On the other hand, for 0 < Γ ≤ 1/5 we have Θ_max = min[Θ₁, Θ₂]. By equating Θ₁ = Θ₂ in equations (31) and (A5), it is straightforward to show that Θ₂ becomes smaller than Θ₁ when Γ is reduced below a critical value Γ_crit = (6 − ϵ_GR)/30, and that this happens at Θ₁ = Θ₂ = 1. This is reflected in Fig. 4 – as we decrease Γ starting from 1/5, the red (ϵ_GR = 0), yellow (ϵ_GR = 3), and green (ϵ_GR = 5) lines transition from solid (Θ₁) to dotted (Θ₂) at the points (1/5, 1), (1/15, 1), and (1/30, 1), respectively. For ϵ_GR > 6 (blue, pink, and black lines in Fig. 4), we have Γ_crit < 0, so this transition never occurs for positive Γ. Finally, for the special value ϵ_GR = 6 we have from equation (A5) that Θ₂ = (1 − 5Γ)⁻¹ > 1, which is obviously greater than Θ₁ = 1/2, so Θ_max = Θ₁. Hence, the brown horizontal line in Fig. 4 extends all the way to Γ → 0.

A2 Fixed points at ω = 0

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 Figs 2(b), (c), and (h).

A3 Does a given orbit librate or circulate?

We can evaluate Δ given ω₀, e₀, Θ, Γ, and ϵ_GR. 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 |$j_0^2 \gt 0$|⁠, then this root can be written as

  $$\begin{eqnarray*} j(\omega =0)&=& -2\frac{\vert q \vert }{q}\sqrt{\frac{j_0^2}{3}} \nonumber \\ & \times &\cosh \left(\frac{1}{3}\mathrm{arccosh}\left[\frac{3 \vert q\vert }{2j_0^2}\sqrt{\frac{3 }{j_0^2}}\right]\right), \end{eqnarray*}$$

  (A10)

  whereas for |$j_0^2 \lt 0$| it is given by

  $$\begin{eqnarray*} j(\omega =0)= -2\sqrt{\frac{-j_0^2}{3}}\sinh \left(\frac{1}{3}\mathrm{arcsinh}\left[\frac{-3q}{2j_0^2}\sqrt{\frac{-3}{j_0^2}}\right]\right). \end{eqnarray*}$$

  (A11)

  Once j(ω = 0) has been determined, the orbit circulates if |$\sqrt{\Theta } \lt j(\omega =0) \lt 1$|⁠, and librates otherwise.
• If Δ > 0 (which necessarily requires |$j_0^2 \gt 0$|⁠), there are three distinct real roots, and they can be expressed as

  $$\begin{eqnarray*} j(\omega =0) = 2\sqrt{\frac{j_0^2}{3}}\cos \left(\frac{1}{3} \cos ^{-1} \left[\frac{-3q}{2j_0^2} \sqrt{\frac{3}{j_0^2}}\right] - \frac{2\pi k}{3} \right), \end{eqnarray*}$$

  (A12)

  for k = 0, 1, and 2. From the theory of polynomial equations, we also know that that the product of the three real roots of (A7) is −q = ϵ_GR/[3(5Γ − 1)] and their sum is 0.

  For Γ > 1/5, this implies that two roots (namely k = 1 and 2) must be negative and one (k = 0) positive. Thus, the orbit circulates if the k = 0 solution lies in |$(\sqrt{\Theta },1)$|⁠, and librates otherwise.

  For 0 < Γ ≤ 1/5, one root (k = 2) must be negative and the other two (k = 0 and 1) positive. If either or both of the two positive roots lie in |$(\sqrt{\Theta },1)$|⁠, then the orbit circulates. If neither of them do then it librates. The case of both positive roots lying in |$(\sqrt{\Theta },1)$| corresponds to two coexisting families of circulating orbits that share values of H*, one above e_f,0 and one below, as in Figs 2(b), (c), and (h). To determine the family of circulating orbits to which the trajectory belongs, we compare its initial eccentricity e₀ with that of the saddle point e_f,0. If e₀ > 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 maximized at ω = 0 are those on circulating phase space trajectories in the regime 0 < Γ ≤ 1/5 (e.g. Fig. 6f; see Paper II for a thorough discussion). The minimum j in this case is j_min = j₀ (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.

APPENDIX C: ANALYTICAL SOLUTION FOR ORBITAL ELEMENTS AT HIGH ECCENTRICITY

In this appendix, we present an analytical solution to the DA equations of motion for all orbital elements in the limit of high eccentricity, assuming Γ > 0. To do this, we will make the following four assumptions:

• j = j_min ≪ 1 is realized at ω = ±π/2,

• Weak or very weak GR, i.e. ϵ_GR ≪ ϵ_weak,

• σ ≪ 1 (equation 50),

• j⁴/Θ ≪ 1.

Assumptions (I)–(III) are familiar from Section 4.4.1 [and of course if we set ϵ_GR = 0 then (II) and (III) are satisfied automatically]. Taken together, assumptions (I)–(III) imply, in particular, that j(t) takes the form (42) with j_min = j₋, j₁ = j₊, and j₂ = j₀, which can be seen by expanding the weak GR equation (55) for σ ≪ 1.

In equation (C8), the value of Ω(0) is an arbitrary constant to be prescribed. Otherwise, equations (C2)–(C3), (C5)–(C8), and the equation j_z(t) = j_z(0) provide a complete, explicit description of the DA dynamics in the high-e limit whenever assumptions (I)–(IV) are satisfied.

C1 Analytical solution in the LK limit

C2 Simplified analytical solution in the limit ϵ_GR = 0

C3 Validity of the analytical 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 |$m_1=m_2=1\, {\rm M}_{\odot }$| orbiting a point mass |$\mathcal {M} = 10^5\, {\rm M}_{\odot }$|⁠. For the outer orbit, we choose a pericentre distance r_p = 0.4 pc and an apocentre distance r_a = 0.6 pc. The outer orbit is then an ellipse with semimajor axis a_g = 0.5 pc and eccentricity e_g = 0.2. For the inner binary orbit, we take the initial conditions e₀ = 0.5, i₀ = 89.75° (so that Θ = 1.4 × 10⁻⁵), and ω₀ = 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 Ω₀ so that Ω(0) ≈ π/2. All that remains is to specify the initial semimajor axis a₀.

In Fig. C2, we take a₀ = 50 au. Then, in each panel we plot the result of the direct numerical integration with a black line and we show the analytical solution (C2), (C5), (C8) with a dashed red line. Panel (a) shows the evolution of log₁₀(1 − e) as a function of time t over a short time interval |$-2.2 \le t/t_\mathrm{min}^{\prime } \le 2.2$| centred on the eccentricity peak. Panel (b) shows the same solution zoomed out over a much longer time interval |$-200 \le t/t_\mathrm{min}^{\prime } \le 200$|⁠. Analogously, panels (c) and (d) show the numerical and analytical solutions for the apsidal angle ω(t) over these same time intervals, while panels (e) and (f) show the evolution of the nodal angle Ω(t). Finally, in panels (c)–(f) we plot blue dashed lines that correspond to the simple non-GR form of the analytical solution, namely equations (C12) and (C13), though to evaluate it we still use the GR-modified value of j_min. (There are also green dashed lines in panels (a) and (b) – see Section 4.4.1.) At the top of the figure, we show the values of various key quantities that allow us to check the validity of the assumptions (I)–(IV).

Overall, in Fig. C2 the analytical 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 |$\sigma , j_\mathrm{min}^4/\Theta \ll 1$| and ϵ_GR ≪ ϵ_weak, meaning that all assumptions (I)–(IV) are fulfilled. Moreover, we see that the full analytical 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 (Section 5.1).

In Fig. C3, we use all the same system parameters as in Fig. C2 except we set a₀ = 30 au. 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 analytical and non-GR solutions in panels (c)–(f). Nevertheless, the analytical solution still matches the numerical one very well for e ≳ 0.9, although not quite as well as in Fig. C2, owing to the fact that σ is now comparable to unity [breaking assumption (III)].

Next, in Fig. C4 we again run the same experiment but this time with a₀ = 15 au. 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), and (C8) largely fail to capture the high eccentricity behaviour even over a very short time-scale. At the same time, we note that the moderate GR solution (62) 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 |$\Phi (r) = -G\mathcal {M}(b+r)^{-1}$|⁠, where the total mass |$\mathcal {M}=10^5\, {\rm M}_{\odot }$| and the scale radius b = 1 pc. (The outer orbit still has r_p = 0.4 pc and r_a = 0.6 pc, but will now fill a 2D annulus rather than forming a closed ellipse – see Paper I.) As a result, we find Γ = 0.235. Also, in this case both σ and κ attain large values, putting our analytical solutions to a demanding test.

In Fig. C5, we integrate exactly the same system as in Fig. C2 except for this replacement of the potential – in particular, we again take a₀ = 50 au. We see that this puts the binary in the weak GR regime ϵ_π/2 < ϵ_GR < ϵ_weak (as in Fig. C3), but that σ is much larger than unity (unlike in Fig. C3). One consequence of this is that the analytical approximation to log₁₀(1 − e) fails rather early on, with significant errors by the time e falls below 0.99 (Fig. C5b). Despite this, the analytical 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 \lt t/t_\mathrm{min}^{\prime } \lt 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 analytical solutions for ω, Ω work well despite assumption (II) being broken.

In Fig. C6, we investigate the same system except with a₀ reduced to 40 au. The binary is still in the weak GR regime but only just so, violating assumption (II). It also has σ ≫ 1 like it did in Fig. C5, violating assumption (III). We see that the analytical fit to log₁₀(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 j_min (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 Fig. C7, in which we have reduced the semimajor axis further to a₀ = 35 au, 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₁₀(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 analytical solution (C8) for Ω can be very accurate even for σ ≫ 1 [violating assumption (III)] provided the behaviour of j(t) in the close vicinity of j_min 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.

Lastly, if all assumptions (I)–(IV) are valid and we additionally have ϵ_GR ≲ ϵ_π/2, then one can employ the simpler non-GR form of the solution for ω, Ω (equations C12 and C13).

APPENDIX D: PHASE SPACE BEHAVIOUR AND MAXIMUM ECCENTRICITY IN Γ ≤ 0 REGIMES

In this appendix, we discuss the dynamical behaviour that arises in negative Γ regimes. This behaviour can be significantly more complicated than for positive Γ. In what follows, we offer an overview of the phase space dynamics for −1/5 < Γ ≤ 0 (in Section D1) and Γ ≤ −1/5 (in Section D2). Lastly, we consider the eccentricity maxima of binaries with negative Γ (Section D3), focusing mainly on initially near-circular orbits.

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 Figs 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 Fig. 2. We now discuss these three figures in turn.

First, we discuss Fig. 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 (Fig. 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 Fig. 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 points 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 Fig. D2 (Γ = −0.15), we find a completely different picture of rather impressive dynamical diversity. In this figure, we have to use 12 panels to fully illustrate the complex phase space behaviour. To begin with, panels (a) and (b) in Fig. D2 have the same morphology as Figs D1(a) and (b). However, panel (c) is very different from Fig. D1(c). 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. |$\dot{\omega }$| 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 and ±π (ω = ±π/2). As ϵ_GR continues to increase, the pair of fixed points at ω = ±π/2 continues 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 Fig. D1(e). The ω = 0 fixed points also disappear once ϵ_GR becomes sufficiently large.

Fig. D3 (Γ = −0.19) shows yet again a different qualitative behaviour. Like in Figs D1 and 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 Fig. 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 Fig. D2. Thereafter, we have qualitatively the same behaviour as in Fig. D2(j).

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 Section 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 Figs D1–D3.

Understanding the fixed points at ω = ±π/2, j = j_f,π/2 is much harder. Just like for Γ > 0, to find j_f,π/2 we must solve the depressed quartic (A1). From this equation, we can once again derive a necessary but insufficient condition for fixed points to exist at ω = ±π/2; however, since 10Γ/(1 + 5Γ) < 0 in this Γ regime, rather than the upper bound ϵ_GR < 6(1 + 5Γ) that we found for Γ > 0 (Section 3.1) we instead get a lower bound, ϵ_GR > 6(1 + 5Γ)Θ^3/2. Unfortunately, it is not easy to write down analogues of the sufficient conditions (29)–(30).²⁴ Indeed, as we saw in Fig. D2, for −1/5 < Γ ≤ 0 fixed points can arise in pairs at ω = π/2 (with another, separate pair at ω = −π/2), corresponding to there being two physical solutions to the quartic (A1).

D2 Phase space behaviour in the case Γ ≤ −1/5

We now turn to the final Γ regime, Γ ≤ −1/5, which luckily is not as complicated as 0 < Γ ≤ 1/5. We need to only illustrate it with a single example, namely Fig. D4, which is for Θ = 0.1 and Γ = −0.5.

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) and (37) numerically to get all seven possible roots, and then declaring j_min 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 e_max for a given i₀, Γ, and ϵ_GR for initially near-circular binaries when Γ ≤ 0 (cf. Section 3.4.1). In Fig. D5, we show e_max(i₀) for various ϵ_GR values. In each panel, we use a different negative value of Γ (cf. Fig. 5).

In the top row of Fig. D5 (panels a–c), we explore the regime −1/5 < Γ ≤ 0. To understand panels (a) and (b), it is worth looking back at Figs D1 and D2 and asking what we expect of the behaviour of initially near-circular orbits. We expect from Figs D1(a)–(d) and D2(a)–(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 Fig. 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 maximized at ω = 0 and is non-zero (Figs D1e and D2j–l).

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. By examining equations (A1) and (A4) in the limit of j_f,π/2 ≈ 1 and small Θ, we find that (ii) becomes true when ϵ_GR is increased beyond the threshold value 6(1 + 5Γ). Putting these constraints together and using the fact that −1/5 < Γ ≤ 0, we find that in the limit Θ → 0 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 Fig. D5(a). Similarly, for Γ = −0.1 we get 3 < ϵ_GR < 9, hence the solo cyan line in Fig. D5(b).

Panels (c)–(f) of Fig. D5 all share a similar morphology, so we will consider them together. In each panel, for a fixed ϵ_GR a finite e_max arises at some critical value of i₀, increases as a function of i₀ until it reaches e_max = e_lim, and then is constant for all larger values of i₀ 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 e_max. The behaviour we see in these panels is consistent with what we expect from the Γ = −0.19 example given in Fig. D3 and the Γ ≤ −1/5 example we studied in Fig. 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 the 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 Figs D4(c) and (d). 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 i₀ – hence the straight horizontal lines in Figs D5(c)–(f). The main qualitative difference between panels (c) and (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 (Figs D1–D3).

Finally, we may briefly compare Fig. D5 with Fig. 5. Consider what happens if we fix ϵ_GR and increase i₀ from zero. In Fig. D5 (Γ ≤ 0), the larger is ϵ_GR, the lower i₀ is required to achieve a non-zero e_max (provided ϵ_GR is not so large that no eccentricity excitation is possible). On the contrary, in Fig. 5 (Γ > 0) the most favourable situation for eccentricity excitation is always to have ϵ_GR as small as possible: The larger ϵ_GR, the larger i₀ 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 i₀ ≈ 90° a larger ϵ_GR always leads to a smaller maximum eccentricity (again provided ϵ_GR is such that eccentricity excitation is possible).

Author notes