Extreme eccentricities of triple systems: Analytic results

Triple stars and compact objects are ubiquitously observed in nature. Their long-term evolution is complex; in particular, the von-Zeipel-Lidov-Kozai (ZLK) mechanism can potentially lead to highly eccentric encounters of the inner binary. Such encounters can lead to a plethora of interacting binary phenomena, as well as stellar and compact-object mergers. Here we find explicit analytical formulae for the maximal eccentricity, $e_{\rm max}$, of the inner binary undergoing ZLK oscillations, where both the test particle limit (parametrised by the inner-to-outer angular momentum ratio $\eta$) and the double-averaging approximation (parametrised by the period ratio, $\epsilon_{\rm SA}$) are relaxed, for circular outer orbits. We recover known results in both limiting cases (either $\eta$ or $\epsilon_{\rm SA} \to 0$) and verify the validity of our model using numerical simulations. We test our results with two accurate numerical N-body codes, $\texttt{Rebound}$ for Newtonian dynamics and $\texttt{Tsunami}$ for general-relativistic (GR) dynamics, and find excellent correspondence. We discuss the implications of our results for stellar triples and both stellar and supermassive triple black hole mergers.

Although the general three-body is notoriously chaotic (Poincaré 1892;Valtonen & Karttunen 2006), in hierarchical systems, where the inner binary is perturbed by a distant companion, perturbative techniques allow certain triple systems to be integrable. In particular, the ZLK may also provide an important channel for driving compact object binaries toward merger as gravitationalwave sources (see, e.g., Mandel & Farmer 2022;Mapelli 2021;Mandel & Broekgaarden 2021 for recent reviews). In this context, ZLK can play a key role in isolated field triples or quadruples (e.g., Silsbee & Tremaine 2017;Antonini et al. 2017;Hamers & Thompson 2019;Fragione et al. 2020;Trani et al. 2022). It may also impact the evolution of merging double compact objects in in dense environments such as globular clusters (Miller & Hamilton 2002;Antonini et al. 2016) or galactic nuclei (e.g., Hoang et al. 2018;Fragione et al. 2019).
The "standard" ZLK integrable mechanism relies on the double-averaging (DA) procedure, where the Hamiltonian is averaged over the mean anomalies of both binaries, such that the evolution occurs over long, secular, timescales (Ford et al. 2000), much longer than both orbital periods. In addition, when the interaction potential is truncated at the leading, quadrupole order, recent studies found analytic solutions in the test particle limit, when one of the inner binary members has low mass with respect to the other member Lubow 2021), and in the general case (Hamers 2021).
When either the DA or the quadrupole approximation break down, the problem becomes chaotic again. The next leading order, the octupole order, applies when the outer orbit is eccentric (i.e. the eccentric ZLK mechanism, Katz et al. 2011;) and when the inner binary masses are unequal. The timescales of octupole effects are generally longer than the quadrupole timescale.
The DA approximation can also break down and cause chaotic evolution, but it was given much less attention. Historical studies of the Lunar motion resolved the tension of additional apsidal precession in terms of additional "evection terms" (Tisserand 1894). The perturbation theory was generalised for irregular satellites on eccentric and inclined orbits (Ćuk & Burns 2004). Recently, Luo et al. (2016) found an effective potential which corrects the DA procedure and makes it more compatible with N-body simulations, while still retains the attractive features of fast integration of the secular equations. Grishin et al. (2017) used the additional corrections to generalise the Hill-stability limit for arbitrary inclination, and Grishin et al. (2018) found corrected formulae for the maximal eccentricity, e max and the critical inclination in the mildly hierarchical case for the test particle limit.
Although the Grishin et al. (2018) extension is limited to cases where the tertiary is much more massive (the Hill approximation) on a circular orbit, direct N-body integrations of more realistic binaries produce more merger rates per galaxy and more eccentric mergers than secular integrations (Fragione et al. 2019).
Motivated by the enhanced merger rates of mildlyhierarchial systems, we take a further step in the analytic description in the mildly-hierarchical three body problem. We relax the test-particle approximation and find a formula for e max for any dynamically stable triple system in an outer circular orbit. Our expression of e max reduces to known results in the literature for the relevant limits and allows analytical estimates for encounters and mergers of triple systems, with direct implications to triple stellar evolution and field triple BH mergers.
The paper is organised as follows. In Section 2 we review the preceding work on secular dynamics of triple systems. In Section 3 we motivate the need for and derive a new analytic formula for the maximal eccentricity and discuss its validity and limitations. In Section 4 we present our results, comparing the analytic expressions to three-body numerical simulations. We discuss the caveats ans astrophysical implications in Section 5. We summarise in Section 6.

OVERVIEW OF SECULAR DYNAMICS
Here we review the governing equations of the standard ZLK mechanism. We then extend the overview to include non-secular contributions that affect the DA formalism as well as relativistic corrections.

Standard von Zeipel-Lidov-Kozai mechanism
Consider an inner binary with masses m 0 and m 1 , semi-major axis a and eccentricity e. The total mass is m bin = m 1 + m 0 . The inner binary is perturbed by an outer companion of mass m out , semi-major axis a out and eccentricity e out . The normal to the outer orbital plane is denoted byn out . The inner binary's eccentricity vector is e and its specific angular momentum vector (expressed in units of the circular angular momentum of a binary with the same orbital semi-major axis) is j = √ 1 − e 2ĵ . The DA quadrupole term in the Hamiltonian (we will sometimes refer to this specific energy as the "potential") that governs the evolution of e and j is where Φ 0 = Gm out m 0 m 1 a 2 m bin a 3 out (1 − e 2 out ) 3/2 .
We omitted terms that contain the semi-major axes, corresponding to the energies of the Keplerian ellipses, which are assumed to be constant. In terms of the argument of pericentre of the inner binary ω and the mutual inclination between the orbital planes i, using j z ≡ j ·n out = √ 1 − e 2 cos i and e z ≡ e ·n out = e sin ω sin i, the quadrupole term can be written as (Naoz 2016 For large initial mutual inclinations, | cos i 0 | < 3/5, the inner eccentricity is coherently excited until a certain threshold (see Eq. (18) in section 3).
Additional extensions to the standard ZLK mechanism are i) including higher-order terms in the multipole expansion, ii) breaking down the DA hierarchical approximation, iii) going beyond the test particle approximation, iv) including additional external forces. We neglect i) for two reasons: the timescales on which the higher-order terms contribute are longer and we focus on circular outer binaries, e out = 0, which sets the octupole term to zero. Under this assumption, combining effects (ii),(iii) and (iv) allows us to find and verify a unified formula for the maximal eccentricity when these effects are taken into account together. We review each extension below.

Corrected Double Averaging
The DA approximation neglects perturbations on timescales shorter than the secular timescale Antognini 2015) τ sec ≈ 1 2π where m tot = m out + m bin . When the system is only mildly hierarchical, the accumulated errors in neglecting these perturbations may be large. It is possible to correct for these errors through the use of the singleaveraged (SA) equations of motion (Luo et al. 2016), which depend on the position of the outer body on its orbit. The parameter that quantifies the level of the hierarchy and these short-timescale perturbations is the SA strength The additional effective "corrected double averaging" potential in terms of the vector elements is which has an axisymmetric term φ circ (j z , e z , e) when the outer orbit is circular, and non-axisymmetric term φ ecc which is important when the outer orbit is significantly eccentric. Explicit expressions for both of these terms are given by Luo et al. (2016). We will discard the φ ecc term, as we are interested in axisymmetric orbits, thus the SA potential is The vector elements are expressed in terms of the Kepler elements in a reference frame where theẑ direction is along the outer angular momentum. Luo et al. (2016) showed that these corrections are consistent with N-body integration results (see e.g. their Fig. 5).
2.3. Relaxing the test particle limit In the test particle limit, the traditionally conserved quantity is the z-component of the angular momentum, j z . In the general case, a modified conserved quantity can be found from the conservation of the angular momentum L tot = L in +L out (Haim & Katz 2018;Anderson et al. 2017;Liu & Lai 2018). Consider the magnitude L out and e out are constant to quadrupole order. The ratio of the circular inner angular momentum to the outer angular momentum is defined as where µ in = m 0 m 1 /m bin and µ out = m bin m out /m tot are the reduced masses of the inner and outer binaries, respectively. Rewriting L in = jηL out , we can rewrite Eq. 8 in terms of the approximately conserved quantity K 2 : We can alternatively use In the limit L in L out , η → 0 and K 1 = K 2 reduce to the familiar conserved quantity j z . We take the initial eccentricity e 0 ∼ 0, so that K 1 ∼ cos i 0 .
We note that the procedure of Luo et al. (2016) corrects for the additional secular terms that arise due to SA . However, there are also osculating orbital elements that vary on outer orbital timescales that do not affect the secular evolution. In other words, the quantity K 2 is conserved on average, but fluctuates around this mean value on a dynamical timescale. Small fluctuations do not affect the (corrected) secular evolution, but become important when the fluctuation δe is larger than its mean value e max (cf. section 3.1).

General Relativistic corrections
When the separation between two bodies is sufficiently small, so that the gravitational radius r g = Gm bin /c 2 is not negligible relative to a(1 − e), general relativity (GR) affects the binary's evolution. Here we include the contribution of GR corrections to the interaction energy at the leading post-Newtonian (PN) order (Blaes et al. 2002;Liu et al. 2015;Liu & Lai 2018): where is the ratio of GR to ZLK precession rates. In the limit GR 1, GR precession is significant enough for the LK mechanism to be suppressed. For large precession e max → 0, and the inner binary essentially evolves in isolation. If the inner separation is small enough, the binary may merge due to gravitational-wave (GW) dissipation within a Hubble time as an isolated binary (Peters 1964;Fragione et al. 2019).
We note that similar estimates can be made for the impact of tidal interactions and rotational bulges (Liu et al. 2015).

COMPUTING THE MAXIMUM ECCENTRICITY
The total interaction term in the specific energy is obtained by adding the contributions of Eqs. 1, 6, and 12: The phase portrait under the Hamiltonian flow governed by the potential in Eq. 14 is still one dimensional, and hence the dynamics are still integrable for a circular outer orbit, e out = 0.
Similarly to Grishin et al. (2018), we evaluate Φ tot at two points on the orbit, but we use the conservation of K 2 (instead of j z ) to account for the general case of nonnegligible inner angular momentum (η > 0), similarly to the approach taken by Anderson et al. (2017) for the purely DA case ( SA = 0). This unified approach leads to an implicit equation for e max , as shown below.
For the initial point on the phase portrait, we specify i 0 and e 0 , and set ω = 0. The latter choice of ω = 0 places the initial condition in a circulating orbit as long as e 0 > 0 1 . The potential at this point is where j 0 = 1 − e 2 0 , and j z = j 0 cos i 0 is evaluated at its initial value.
The second point is evaluated at e max , which is attained at ω = ±π/2 for axisymmetric potentials (e.g. Hamilton & Rafikov 2019). The potential at this point is where j min = 1 − e 2 max and cos i m is the mutual inclination at e = e max . We assume a small initial eccentricity (e 0 1) in order to focus on how circular binaries become eccentric. Conservation of K 1 , Eq. (11), yields cos i m = (K 1 + e 2 max η/2)/j min , where K 1 ≈ cos i 0 . Equating Φ 1 = Φ 2 allows us to find an implicit expression for e max in the appropriate ZLK window (Anderson et al. 2017): where g max = cos i 0 + ηe 2 max /2. Eq. (17) allows us to extract (albeit implicitly) e max in terms of SA , η, i 0 and GR , which are all fixed and determined from initial conditions.
In the test particle limit, without corrected DA and in the purely Newtonian regime, the usual expression (Naoz 2016) for the maximal eccentricity is obtained:

Fluctuation term
We have derived an expression for e max by equating the potential, including the SA term, at the initial and maximum-eccentricity points, assuming a conserved K 2 . However, K 2 is only conserved on average, and actually fluctuates around this averaged value, where the fluctuating amplitude is (Haim & Katz 2018) This results in a fluctuating eccentricity (δe) about the mean value (e max ). The total maximal eccentricity is where δj is the fluctuation in j. Differentiating K 2 from Eq. (10), we find ∆K 2 = (cos i m + jη)δj.
The inclination i m at which e max is attained falls on the solid black lines in Fig. 1. We focus on configurations where η 0.5 (see Sec. 5.1), and we approximate cos i m as ± 3/5. The next contribution is of order ∼ η/10 O(1), so it can be neglected. The maximal fluctuation in j z is (Haim & Katz 2018;Grishin et al. 2018) Using Eqs. (19), (21) and (22) at e max , δj is δj = 1 + |K 2 |η Finally, for the fluctuating eccentricity, we square Eq. (20), neglect the (δe) 2 term (assuming that e max δe), and substitute our expression for δj from Eq. (23): The fluctuating eccentricity is a function of i 0 , η, SA , and e max , which are all determined from initial conditions. For a specific domain of initial inclination the fluctuation becomes sufficiently large such that the maximal eccentricity e max + δe → 1 and the binary merges or becomes unstable.
The treatment leading to Eq. (24) is valid for mildly hierarchical triples with small oscillations δe e max . The regime of validity of this treatment can be determined by considering considering when |K 2 | ≤ |∆K 2 |, at which time the perturbations become so large that the eccentricity is not well defined: Figure 2. Initial conditions in semi-major axis ratio and mass ratio space (m0 = m1). Dashed lines indicate fixed angular momentum ratio η (Eq. 9). Solid lines indicate fixed hierarchy parameter SA (Eq. 5). Black diamonds denote initial conditions used for the set of N-body simulations, numbered as in Table 1 Fig. 2. We set m0 = m1 = 50M and GR = 0.0001, with the inner semi-major axis a, tertiary mass mout and semi-major axis aout following from the chosen ( SA, η).
We note that in the test particle limit (η → 0), the RHS of Eq. 25 reduces to (9/8) SA as determined by Grishin et al. (2018).
We have developed an expression for e max (Eq. 17) and the fluctuation magnitude δe (Eq. 24) relative to e max . Fig. 1 shows the maximum eccentricity e max + δe of the inner binary, with initial mutual inclination i 0 to the tertiary body with inner-to-outer angular momentum ratio η (Eq. 9). The quantity 1 − (e max + δe) is plotted as a colour map on a logarithmic scale, for various choices of SA = {0.03, 0.08, 0.1, 0.15}. The vertical red lines indicate specific choices of η to compare against N-body simulations (see Fig. 2, 3). The black lines indicate the boundaries of zero and non-zero eccentricity. The solid lines are the same as in Fig. 1 of Anderson et al. (2017). The region between solid grey lines shows where the fluctuation becomes large and the eccentricity will approach unity (∆K 2 > K 2 ). In the limit of SA → 0 we expect to reproduce the top left panel of Anderson et al. (2017), corresponding to GR = 0 (see also next section). In the limit of η → 0 we get back the results of Grishin et al. (2018).
We are now ready to test our Newtonian analytic theory and relativistic extension via numerical simulations in the next section.

Numerical codes
We now seek to compare the analytic formalism we have developed in previous sections to N-body integration. We employ two different numerical integrators: rebound (Rein & Liu 2012) and tsunami (Trani 2020). Within rebound, we use ias15, a fast, adaptive, higher-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits (Rein & Spiegel 2015). tsunami is a fast and accurate few-body code designed to follow the evolution of selfgravitating systems. The integrator is based on Mikkola & Tanikawa (1999) algorithmic regularisation and can easily handle close encounters, highly hierarchical sys- Figure 3. Comparison of analytical emax (black) and emax + δe (red) prescriptions against N-body numerical simulations using rebound (Rein & Liu 2012, blue dots) and tsunami (Trani 2020, green dots). The grey region shows where ∆K2 > K2, i.e., the eccentricity approaches unity, leading the inner binary to merge or the system to become unstable. The initial conditions are shown in tems, and extreme mass ratios. tsunami also includes post-Newtonian (PN) corrections through 2.5 PN order, i.e., including both relativistic precession and lowestorder radiation reaction, allowing us to explore systems which have large GR . We set a collision radius equal to five times the Schwarzschild radius of the particles in order to avoid unphysical behavior when the PN expansion breaks down. Whenever the distance between two particles is less than the sum of their collision radii, the simulation is stopped. Because the triples we consider are close to the instability regime (see Figure 2), it is possible that some systems will break up over the course of the simulation. We stop a simulation and label it as unstable if the binding energy of either the inner or outer binary becomes negative. Otherwise, the end time of all runs is 500 times the outer orbital period.

Initial conditions
The set of initial conditions we have chosen to compare against N-body integrations are shown in Fig. 2 and Table 1. Fig. 2 shows initial conditions in semimajor axis ratio and mass ratio space for m 0 = m 1 . Dashed lines indicate fixed angular momentum ratio η (Eq. 9). Solid lines indicate fixed hierarchy parameter SA (Eq. 5). Black diamonds denote the set of N-body Figure 4. Comparison of analytical emax (black) and emax + δe (red) prescriptions against N-body numerical simulations using tsunami (Trani 2020, green dots, crosses). Dots are indicative of initial inclinations for which the system remains stable. Crosses indicate initial inclinations that lead to collision or break-up of the inner binary. The grey region shows where ∆K2 > K2, i.e., the eccentricity is formally unconstrained, leading the inner binary to merge or the system to become unstable. The initial conditions follow Table 1, with both separations scaled by a factor of 10 −3 such that GR = 0.1 in all cases. The purple dashed line is the expected limit at small GR 1, emax ≈ 1 − (8 GR/9) 2 ≈ 0.996.

simulations, numbered as in
For non-negligible values of both SA and η, the systems must be compact and close to the stability boundary. The Mardling & Aarseth (2001) stability limit is too strict for m out m 1 , so simulation 4 is still stable.
For a concrete example, we consider an equal mass inner binary of 50M BHs. The separation of the inner and outer binaries and mass of the tertiary body are then determined from the choice of η, SA and GR . All simulations start with e = 0, Ω = π/4, ω = π/2 and f = 0. The outer orbit has e out = 10 −5 , Ω out = ω out = f out = 0. Ω is the longitude of ascending node. f is the true anomaly. This ensures that the osculating j z is at its mean value. Fig. 3 shows the comparison of our analytic prescription against numerical integrations. The initial conditions are the red vertical slices in the corresponding panels of Fig. 1. The analytic curves e max (black) and e max + δe (red) are compared against N-body integrations computed using rebound (blue) and tsunami (green). The simulations closely follow the e max + δe curve for a range of values of SA and η. Furthermore, Nbody simulations match the region we have determined (Eq. 25) within which the eccentricity becomes unconstrained (grey region). The width of the grey region is proportional to the magnitude of SA , and the grey region appears to be positioned at larger initial inclinations for larger values of η.

Including GR effects
We now consider systems with GR = 0.1, shown in Fig. 4. In order to achieve such a large value of GR , the inner separation is of order 10 −3 AU, so even a circular binary with 50M components will merge in ∼ 1 year through GW emission (Peters 1964). On the other hand, the merger may be eccentric, hence it is instructive to study this case.
Comparing Figs. 3 and 4, we notice that e max curves (black) only reach values of order 1 − (8 GR /9) 2 ≈ 0.996 in the GR case as relativistic precession suppresses the ZLK mechanism, compared to eccentricity approaching unity in the Newtonian case (cf. Eq. (17) with SA = η = 0). This is consistent with the result of Liu et al. (2015).
We notice similar behaviour when examining the e max + δe (red) curve, leading to a "shrinking" of the region where ∆K 2 > K 2 . The criterion we have developed (Eq. 25) no longer follows the regions within which the eccentricity becomes unconstrained. In the case of η = 0.5, SA = 0.03, the red curve indicates that the grey region ought to disappear entirely as GR effects prevent the inner binary from flipping. Some orbits, especially the highly eccentric ones, merge through GW emission or come sufficiently close to be labeled as collisions (green crosses in Fig. 4) within the duration of the simulation. At an eccentricity of 0.99, the merger timescale shrinks by a factor of almost one million relative to a circular binary, so that a binary consisting of two 50M BHs at a separation of 0.001 AU merges in approximately 1 minute (Peters 1964;Mandel 2021). Binaries that reach more moderate maximal eccentricities will be circularised prior to the merger.
Overall, we find good correspondence between the analytic and numerical results. The only discrepancy is that occasionally maximal eccentricities appear to be slightly larger than predicted analytically. We find no difference in maximal eccentricity between 1PN only and full 2.5PN expansions, although energy dissipation through GW emission in the 2.5 PN expansion brings the system to merger.

Astrophysical implications
General picture: We derive and validate a new formula for the maximal eccentricity of the inner binary, implicitly given in Eq. (17).
Perhaps the main effect is the expansion of the parameter space in which mildly hierarchical systems can reach very high eccentricity. Grishin et al. (2018) discussed this in the test particle limit, and we find the same effect for triples of comparable masses. In the hierarchical, test-particle, Newtonian limit, when SA = GR = η = 0, Eq. (17) yields a range | cos i 0 | < (3/5)(1 − e 2 max ) of initial inclinations that can yield eccentricities e max or greater. This range in cos i 0 is only ≈ 0.07 for e max > 0.999, as illustrated by the e DA max curve in the top left panel of Fig. 3. On the other hand, relaxing the hierarchical limit, Eq. (25) predicts that arbitrarily large eccentricities can be reached for | cos i 0 | < (9/8) SA (with GR = η = 0) -i.e., a range of 0.34 in cos i 0 for SA = 0.15. This is illustrated by the thickness of the grey region in the bottom left panel of Fig. 3. Thus, high maximal eccentricities are significantly more likely for mildly hierarchical systems, assuming an isotropic (flat in cos i 0 ) distribution of initial inclination angles between the inner and outer orbits.
The extension beyond the test particle limit shifts the initial inclinations that lead to e max → 1 toward retrograde configurations. The range of cosines of initial inclination angles leading to unbound eccentricity is increased by a fraction (1 − 9 SA η/8) −1 (see Eq. 25), but it is a mild correction in most cases.
The extreme case of η 1 requires the tertiary to be much less massive than the inner binary, thus the inner binary is essentially unperturbed. This is the "inverse ZLK problem" (Naoz et al. 2017) and is not discussed here.
We tested our results both with Newtonian and PN codes, where the strength of GR is encapsulated in GR (Eq. 13). It is important to stress that this is the relative strength of GR precession to ZLK precession, and not a statement on the proximity to the gravitational radius and the breakdown of the PN expansion.
If all three masses are similar and the outer orbit is circular, Eq. (5) reduces to SA ∼ (a/a out ) 3/2 , Eq. (9) to η ∼ (a/a out ) 1/2 , and Eq. (13) to GR ∼ (a/a out ) −3 (r g /a). Stability requires that a out is at least a few times larger than a, so in this regime, non-negligible η and SA are generally obtained when a is of order 0.1a out . A significant GR then requires a 10 4 r g . Once eccentricity is driven up to large values by ZLK resonances, the inner binary periapsis in such systems becomes sufficiently small that they merge rapidly through GW emission. This can leave very limited opportunities for interactions that put systems into this regime of interest.
Stellar BH and stars: For the case of stellar mass BHs, in order for for SA , η and GR to all be non-negligible, requires extremely compact systems that merge within thousands of years. Isolated field triples could fit into a compact system for a out 100R (Vigna-Gómez et al. 2021), which corresponds to GR ∼ few × 10 −3 , so the treatment can still be Newtonian. However, such compact triples are expected to be coplanar due to interactions during stellar evolution. Stellar evolution may weaken the hierarchy of otherwise stable triples (Perets & Kratter 2012;Toonen et al. 2021). Such triples will enter the semi-secular regime before becoming dynamically unstable, which may result in an enhanced rate of mergers and collisions. Dense environments such as globular clusters can in principle construct such compact triples and lead to eccentric mergers. (Fragione et al. 2019), as potentially inferred for GW190521 (Romero-Shaw et al. 2020).
SMBHs: Supermassive triple BHs are more promising in having larger values of GR for longer, since the GW-driven inspiral duration scales with the mass when the separation in units of the gravitational radius, a/r g , is fixed. This leaves more time for dynamical or gas interactions to insert systems into the regime of interest. Bonetti et al. (2016) show that ∼ 10 8 M BHs may stall at radii of around a pc and triple-induced mergers are plausible. The fiducial system in their Fig. 6 has GR ≈ 0.7, η ≈ 0.2, SA ≈ 0.02. In a follow-up Monte Carlo study, Bonetti et al. (2018) found that between 20 − 30% of the initially sampled triple BH systems are merging, with preference for equal mass inner binary (see their Fig. 6.). Without PN evolution, the merger rate is around 40 − 60%, roughly twice than the rate with PN terms.

Limitations and caveats
Outer eccentricity: We limited our study to outer circular binaries, e out = 0. For non-zero eccentricity, the corrected averaging term has a non-axisymmetric contribution, hence K 2 is no longer conserved and the dynamics are chaotic. Moreover, for unequal inner masses, the octupole term may also induce chaotic evolution at longer timescales. In the DA test-particle limit, e max is then unconstrained for a wider range of initial incli-nations, depending on the octupole strength (e.g Naoz et al. 2011;Katz et al. 2011;Muñoz et al. 2016). We expect the maximal eccentricity to be also increasing compared to outer circular rbits on average. A systematic study is deferred for future work.
Tides, stellar evolution and other complications: A similar treatment can account for tides when the inner binary contains a non-degenerate star. Equilibrium tides (Hut 1981) cause extra precession and can also quench ZLK oscillations, similarly to GR. The strength of these tides depends on the radii and apsidal motion constants of the stars, and e max had been calculated analytically in the DA regime (Liu et al. 2015). It is possible to extend the study to main sequence or giant stars which have a convective envelope. Similarly to GW dissipation, tidal dissipation could also be important and reduce the separation. We do not consider it here.

Future work
The interplay of eccentricity, spin and orbital apdisal precession affects the last stages of compact BH evolution (Phukon et al. 2019). It is possible that the librating argument of pericentre around π/2 in the ZLK regime could serve as a physically motivated prior and improve eccentric GW templates for the analysis of mergers from dynamical formation channels.
Analytical results can play a key role in other largescale numerical modelling, such as population synthesis of triple systems. Our results could serve as a prescription and save the computational cost of numerically integrating orbital dynamics.

SUMMARY
We have derived an analytical expression for the maximal eccentricity of hierarchical triple systems with an outer circular orbit. The expression is valid for any angular momentum ratio and any orbital period ratio (provided that the system is dynamically stable), and can thus be applied to any triple system which satisfies the assumptions of a circular outer orbit and dynamical stability. The main results are Eqs. (17), (24) and Fig. 1. This contributes to the analytic understanding of the long term stability and evolution of triples of field stars or BHs of any mass, and can be used as a prescription in population synthesis codes.
After verifying and recovering previous results, we focus on mild hierarchy ( SA 0.05, Eq. 5) and comparable masses (η 0.1, Eq. 9). These systems are inherently close to being dynamically unstable (Fig. 2). Nevertheless, we are able to obtain accurate results with good agreement with numerical simulations (Fig. 3, 4) in both Newtonian and PN limits. Similarly to the test particle case, mild hierarchy (moderately large SA ) is mainly responsible for the increase in e max and the range of initial inclinations over which it can be achieved, while the angular momentum ratio η mostly shifts the initial inclinations contributing to maximum eccentricities towards retrograde orbits.
We expect that GR corrections to triple dynamics are small for stellar mass BHs, and Newtonian treatments should be sufficient. On the other hand, SMBH triples do require PN corrections. Other stars and compact objects may also experience enhanced rates of collisions or mergers via the ZLK mechanism in the mildly hierarchical, non-test-particle limit, though additional aspects of tides and stellar evolution need to be taken into account.