Empirical stability boundary for hierarchical triples

Abstract

The three-body problem is famously chaotic, with no closed-form analytical solutions. However, hierarchical systems of three or more bodies can be stable over indefinite timescales. A system is considered hierarchical if the bodies can be divided into separate two-body orbits with distinct time and length scales, such that one orbit is only mildly affected by the gravitation of the other bodies. Previous work has mapped the stability of such systems at varying resolutions over a limited range of parameters, and attempts have been made to derive analytic and semi-analytic stability boundary fits to explain the observed phenomena. Certain regimes are understood relatively well. However, there are large regions of the parameter space which remain unmapped, and for which the stability boundary is poorly understood. We present a comprehensive numerical study of the stability boundary of hierarchical triples over a range of initial parameters. Specifically, we consider the mass ratio of the inner binary to the outer third body ( $q_\mathrm{out}$ ), mutual inclination (i), initial mean anomaly and eccentricity of both the inner and outer binaries ( $e_\mathrm{in}$ and $e_\mathrm{out}$ respectively). We fit the dependence of the stability boundary on $q_\mathrm{ out}$ as a threshold on the ratio of the inner binary’s semi-major axis to the outer binary’s pericentre separation $a_\mathrm{in}/R_\mathrm{p, out} \leq 10^{-0.6 + 0.04q_\mathrm{out}}q_\mathrm{out}^{0.32+0.1q_\mathrm{out}}$ for coplanar prograde systems. We develop an additional factor to account for mutual inclination. The resulting fit predicts the stability of $10^4$ orbits randomly initialised close to the stability boundary with $87.7\%$ accuracy.

1. Introduction

The three-body problem had been studied for over 100 years (Poincaré & Magini 1899) and is notoriously chaotic, with no analytic solution (Poincaré 1892). For strongly chaotic systems, a statistical approach is generally preferred to yield a probabilistic outcome of the scattering in terms of the final states (Monaghan 1976a,b). Only recently, the explicit dependence on the orbital elements was found with various methods (e.g. Stone & Leigh 2019; Ginat & Perets 2021a,b; Kol 2021a,b; Manwadkar et al. 2021).

Hierarchical systems are a subset of three-body systems such that the system can be divided into an inner and an outer binary, each only weakly perturbed by the other. Hierarchical systems are deemed stable if they maintain their hierarchical structure without collision or ejection for a large number of orbits. The long-term evolution of hierarchical systems had been explored extensively over the years (Harrington 1968; Ford, Kozinsky, & Rasio 2000; Naoz et al. 2013; Toonen et al. 2020).

Hierarchical systems are the only multiple systems observed in nature, either for triple stars of comparable masses (where the majority of massive stars are actually in triples or higher multiples, Duchêne & Kraus 2013; Moe & Di Stefano 2017), or systems with extreme mass ratios where one of the masses is much smaller than the other ones. Examples of such extreme mass ratio systems include (but are not limited to) circumbinary planets (Doyle et al. 2011), moons and binary asteroids (Richardson & Walsh 2006; Naoz, Perets, & Ragozzine 2010), or binary stars/black holes around an intermediate or super-massive black hole (Fragione et al. 2019).

The boundary between stable hierarchical and unstable chaotic systems has a non-negligible width. For an inner separation $a_\mathrm{in}$ , and an outer separation $a_\mathrm{ out}$ , the stability boundary (SB) will generally depend on the separation ratio $\alpha \equiv a_\mathrm{in}/a_\mathrm{out}$ . It is expected that for $\alpha \sim 1$ the system is strongly chaotic, while it should be stable for $\alpha \ll 1$ . The exact SB depends also on the masses of the bodies, as well as other orbital elements of both binaries. Due to the large parameter space and inherent chaotic nature of the problem, traditional studies explored the SB over a limited parameter space (Eggleton & Kiseleva 1995; Holman & Wiegert 1999; Petrovich 2015; He & Petrovich 2018; Quarles et al. 2018), while new methods explore the boundary using a machine learning approach for specific choices of initial conditions (Lalande & Trani 2022; Vynatheya et al. 2022, hereafter V+22). Of particular importance are the outer mass ratio $q_\mathrm{out} \equiv m_\mathrm{in}/m_\mathrm{out}$ , and the mutual inclination between the two orbits, i. Here $m_\mathrm{in}$ is the total mass of the inner binary and $m_\mathrm{out}$ is the mass of the outer companion.

For extreme mass ratios $q_\mathrm{out}\ll1$ , the Hill radius $r_\mathrm{H} \equiv a_\mathrm{out}\!\left(q_\mathrm{out}/3\right)^{1/3}\ll a_\mathrm{out}$ defines the length scale of the SB (Hill 1878). Tight inner binary orbits with $a_\mathrm{in} \ll r_\mathrm{H}$ will be stable while wide orbits with $a_\mathrm{in} \gtrsim r_\mathrm{H}$ will be unstable. Numerical studies have repeatedly shown that the SB lies between between $0.4r_H-0.9r_H$ (Henon 1970; Hamilton & Burns 1991). The inclination dependence of the SB was examined analytically by Innanen (1979, 1980), Mylläri et al. (2018), who explained the greater stability of retrograde orbits over prograde orbits by pointing out the asymmetric role of the Coriolis acceleration in stabilising the former and destabilising the latter. Secular analysis reveals more complicated behaviours such as the Lidov-Kozai (LK) mechanism in mutually inclined orbits, in which the outer perturber exerts a torque on the inner binary, exciting oscillations in inclination and eccentricity of the inner orbit (Lidov 1962; Kozai 1962). This reduces the critical stability radius for mutually inclined orbits. Grishin et al. (2017) generated a semi-analytic fit of the SB for arbitrary inclinations which explains the orbital distribution of irregular satellites of the giant planets (Carruba et al. 2002).

Figure 1. Fraction of stable orbits out of 20 evenly spaced mean anomalies against mutual inclination for $q_\mathrm{out}=1$ (left) and $q_\mathrm{out}=10^{-5}$ (right). The left ordinate axis shows $a_\mathrm{in}$ in units of $a_\mathrm{out}$ , while the right ordinate shows $a_\mathrm{in}$ in units of the Hill radius $r_\mathrm{H}$ .

For comparable masses, $q_\mathrm{out} \lesssim 1$ , instability (and later chaotic evolution and ejection) arises due to resonant interactions between the inner orbit and harmonics of the outer orbit (Mardling 2008, 2013). Due to angular momentum exchange, the outer orbit gains eccentricity, which induces higher harmonics. The structure of overlaps of neighbouring resonances leads to chaotic interactions (Chirikov 1979) and can be used to estimate stability. Mardling & Aarseth (2001, hereafter MA01) estimated the critical separation ratio for stability in terms of the mass ratio as $\alpha \propto (1+1/q_\mathrm{out})^{-2/5} \propto (q_\mathrm{out}/(1+q_\mathrm{out}))^{2/5}$ . Although strictly valid for comparable masses, $q_\mathrm{out} \gtrsim 0.2$ , numerous studies extrapolate the MA01 stability boundary to very high or very low mass ratios. Moreover, the MA01 inclination dependence is not particularly accurate. Very recently, V+22 updated the stability boundary to include non-monotonic inclination dependence.

Each empirical fit relies on an extensive parameter study and specific assumptions in the relevant regime which limit its domain of validity. However, unifying two or more regimes is challenging due to the very distinct underlying physical processes. For example, while the comparable masses case relies on exciting eccentricity in the outer orbit, the extreme mass ratio case cannot change the eccentricity of the outer orbit at all, and the instability is due to tidally shearing the inner binary. In addition, many of the recent population studies regularly check for violation of the SB during the dynamical evolution of a triple system (Toonen et al. 2020; Hamers et al. 2021; Grishin & Perets 2022). Finally, recent observations find companions with mass as low as $0.2\,\mathrm{M}_\odot$ to massive O-stars, pushing the mass ratio to extreme values (Reggiani et al. 2022). Given the large multiplicity of massive stars, it is plausible to find triples with extreme mass ratios. A unified stability boundary for wider range of mass ratios, inclinations and other parameters is compelling.

In this paper, we find a unified empirical fit for the stability boundary for any outer mass ratio, $q_\mathrm{out}$ , and extend it to any mutual inclination. By doing so, we bridge together the $1/3$ power law for the extreme mass ratio regime and the $2/5$ regime for comparable masses. We also extend the inclination dependence fit and explore the dependence on the eccentricities and other orbital parameters. Our paper is organised as follows: Section 2 describes the numerical methods and initial conditions. Section 3 presents the results and the unified stability boundary fit, which is given in Equations (1) through (4). In Section 4, we discuss the implications of our findings, and assess the utility of our fit against other contemporary models, and Section 5 summarises our findings.

2. Methods

All integrations were computed using the REBOUND N-body code, with the IAS15 adaptive symplectic integrator, accurate to machine precision for a billion orbits (Rein & Spiegel 2014). An orbit was considered disrupted if the distance between the inner binary components $m_1$ and $m_2$ exceeded half the outer pericentre separation $r_\mathrm{p, out} / 2 \equiv a_\mathrm{out} (1 - e_\mathrm{out}) / 2$ . In orbits with $q_\mathrm{out} \le 1$ , the inner binary can be disrupted before its separation exceeds $a_\mathrm{out}/2$ . However in practice such orbits soon reach this separation, and the above criterion is sufficient to define stability, thus it is a conservative definition of stability. Reaching inner orbital separations of half the outer orbit’s pericentre also requires significant changes in the orbital energy, which coincides with other definitions of stability (e.g. V+22). We compare different stability thresholds in Section 4.2. Systems with large values of $q_\mathrm{out } \gg 1$ can be non-hierarchical and stable, where our criterion is not satisfied (Bhaskar et al. 2021, V+22).

An orbit was classified as stable if it survived for 100 outer orbital periods ( $P_\mathrm{out}$ ) without disruption, following previous work. Grishin et al. (2017) integrated for $100 P_\mathrm{out}$ to incorporate the impact of secular changes that occur on timescales of $\lesssim 10 P_\mathrm{out}$ for systems close to the stability boundary. Grishin et al. (2017) also ran one grid of initial conditions for longer times of $10^4 P_\mathrm{out}$ and obtained similar results. V+22 also used a simulation duration of $100 P_\mathrm{out}$ , arguing that 99% of the unstable orbits are classified as such within $100 P_\mathrm{out}$ (see their Figure 1). However, the actual timescale for instability may be sensitive to initial conditions. While this paper was under review, Hayashi, Trani, & Suto (2022) integrated a limited sample of initial conditions for much longer timescales to show that some orbits (especially orbits with large $e_\mathrm{out}$ and large mutual inclinations) could reach instability much later, after $\gtrsim 10^4P_\mathrm{ out}$ .

2.1 Initialisation

The system was initialised in Jacobi coordinates, described as follows. $m_1$ was initialised at the centre of the system initially with zero velocity. REBOUND initialisation tools were used to place $m_2$ in an orbit defined by $a_\mathrm{in}$ and $e_\mathrm{in}$ . Finally, $m_\mathrm{out}$ was initialised in orbit around the centre-of-mass of $m_1$ and $m_2$ , based on $a_\mathrm{out}$ , $e_\mathrm{out}$ and inclination i.

Each plot shows a high-resolution grid of systems, with one input parameter on the abscissa and a suitably scaled $a_\mathrm{in}$ on the ordinate. Eccentric orbits were initialised at apastron for greater precision. Since orbits at different mean anomalies exhibit slightly different stability, systems with 20 evenly-spaced mean anomalies were chosen for each grid coordinate, and the plots show the fraction of these systems that were stable at each coordinate. Red areas exhibited a high fraction of stable orbits, while the unstable regions are blue. For a detailed breakdown of initial conditions, see Table 1.

In order to obtain our results, we overall integrated over 5 million orbits. The last row of Table 1 lists the number of panels for each figure times the number of runs to generate each panel.

Table 1. Initial conditions for the plots presented in this paper. Variables represented by an interval are distributed uniformly unless given in exponential form (i.e. [ $10^{-2}, 10^0$ ]), in which case they are distributed logarithmically. The value of $a_\mathrm{in}$ on the ordinate of each plot is scaled as appropriate to the figure.

3. Results

Here, we explore the dependence of the SB on the various system parameters: the inclination (Section 3.1), mass ratios (Section 3.2) and eccentricity (Section 3.3). We present our empirical fits in Section 3.4. All logarithms in this paper are base 10.

3.1 Inclinations

Figure 1 shows the SB dependence on the mutual inclination. It closely follows previously reported non-monotonic fits (Hamilton & Burns 1991; Grishin et al. 2017). The boundary is essentially identical to the $q_\mathrm{out}\ll 1$ Hill approximation up to $\log\![q_\mathrm{out}]\leq-1$ , with prograde orbits stable up to around $0.5r_\mathrm{H}$ , and retrograde orbits as far as $0.8r_\mathrm{H}$ . Stability peaks at inclinations of $60^\circ$ , and reaches maximal values again near $180^\circ$ . The $\sim 90^\circ$ inclinations are the least stable, as expected from the prominent LK resonances that excite large eccentricities in this regime.

The left panel of Figure 1 shows the $q_\mathrm{out}=1$ case. Here the prograde and retrograde SB’s are similar with a small dip in the high inclination range. This is because the Hill approximation does not apply and instead the critical separation is somewhere between $\alpha=a_\mathrm{in}/a_\mathrm{out} \approx 0.3$ (for prograde) and $\alpha \approx 0.4$ – $0.55$ (for retrograde), compatible with the MA01 and V+22 stability boundaries for both prograde and retrograde cases.

The SB also peaks around inclination $50^\circ$ more smoothly in the $q_\mathrm{out}=1$ case, while the decay is more abrupt in extreme $q_\mathrm{out}$ cases and occurs close to $60^\circ$ . The increase of the critical inclination for the onset of LK oscillations for triples with mild hierarchy (close to their stability limit) has been observed by Grishin et al. (2017), was later derived analytically by Grishin, Perets, & Fragione (2018) for small $q_\mathrm{out}$ , and was recently extended to any mass ratio by Mangipudi et al. (2022). These results are consistent with the analytic expectations.

3.2 Mass ratios

Figure 2 shows the dependence of the SB of the mass ratio $q_\mathrm{out}$ . In the Hill approximation, where $\log q_\mathrm{out}<-1$ , the power law of $1/3$ dominates across all inclinations. The inclination-dependent boundary of equal-mass systems manifests in a relative increase in the gradient of the SB towards $\log q_\mathrm{out} = 0$ in the low-inclination systems (left panel), and a corresponding decrease in the gradient of the stability boundary for more inclined systems (right panel).

Figure 2. Fraction of stable orbits out of 20 evenly spaced mean anomalies against the mass ratio $q_\mathrm{out}$ , for inclinations $i=0$ (left) and $i=150^\circ$ (right). The white line indicates the Hill radius $r_\mathrm{H}$ in both panels. At extreme mass ratios, the relationship is consistent with the Hill regime, while at $q_\mathrm{out} \geq 0.1$ , high inclination orbits become relatively less stable, while low inclination orbits become more stable, although the inclination dependence is not monotonic (see Figure 1).

3.3 Inner and outer eccentricities

Previous work has demonstrated that in systems with $q_\mathrm{out} \sim 1$ , the stability boundary of hierarchical triples depends on both outer and inner eccentricity (MA01; V+22). Preliminary consideration suggests that the circular stability boundary condition may be translated into one for a general eccentric boundary by considering the ratio between semi-major axis of the inner orbit $a_\mathrm{in}$ and the Hill pericentre radius $r_\mathrm{H,p} \equiv r_\mathrm{H} (1 - e_\mathrm{out})$ , rather than $a_\mathrm{in}/r_\mathrm{H}$ . The impulse from the outer object is greatest at its pericentre. However, the eccentricity of the inner orbit is determined by the three-body interactions rather than by its initial value for $q_\mathrm{out} < 1$ (e.g., the maximum eccentricity in the Lidov-Kozai regime is set by the masses, separations and the mutual inclination Grishin et al. 2017; Mangipudi et al. 2022).

Figure 3 shows the dependence of the SB on the initial eccentricity, for $q_\mathrm{in}=q_\mathrm{out}=0.01$ . The left panel shows that the SB is largely independent of $e_\mathrm{in}$ . The right panels shows that the dependence of the SB on the outer binary’s eccentricity is indeed well approximated by using $r_\mathrm{H,p}$ in place of $r_\mathrm{H}$ in the SB threshold. Both of these approximations begin to break down when $q_\mathrm{out} \gtrsim1$ . The initial eccentricity of the inner binary is more significant when $q_\mathrm{out} \gtrsim 1$ ; in this case, the apocentre separation of the inner orbit becomes a relevant parameter $r_\mathrm{a, in} \equiv a_\mathrm{in} (1+e_\mathrm{in})$ , since $m_2$ is most weakly attracted to $m_1$ at apocentre (see (MA01; V+22) for a discussion of eccentricity in this regime).

Figure 3. Fraction of stable orbits out of 20 evenly spaced mean anomalies against $e_\mathrm{in}$ (left) and $e_\mathrm{out}$ (right), plotted for $q_\mathrm{out}=0.01$ . In the left panel, $e_\mathrm{out} = 0$ , and in the right panel $e_\mathrm{in} = 0$ .

3.4 Fits

Constructing a fit to a multivariate parameter space entails significant challenges, due to the numerous potential correlations between the various parameters. In order to tackle this problem, we first attempt to approximately describe the SB by the product of single-variable functions for each parameter, since the plots for inclination look similar at different mass ratios. Furthermore, since $e_\mathrm{in}$ is mostly irrelevant and the effects of $e_\mathrm{out}$ can be explained by considering the critical separation as a fraction of the outer periastron separation, the problem is reduced to only two variables, inclination and mass ratio:

(1) \begin{equation} \frac{a_\mathrm{in}}{R_\mathrm{p, out}} = f(q_\mathrm{out})\cdot g(i) \cdot h(q_\mathrm{out}, i). \end{equation}

A piecewise linear regression with two segments was used to describe the log–log mass ratio boundary accurately in the low-inclination limit, and a smooth interpolation between these two regimes provides the first factor, $f(q_\mathrm{out})$ :

(2) \begin{equation} f(q_\mathrm{out}) = 10^{-0.6 + 0.04q_\mathrm{out}}q_\mathrm{out}^{0.32+0.1q_\mathrm{out}}, \end{equation}

where the numerical pre-factor is chosen to optimise classification accuracy across eccentricities (see Section 4 for more details on classification accuracy).

The inclination boundary is more complicated, and requires more careful treatment. A cosine function describes the curve up to $60^\circ$ , while the high-inclination region is fit by a fourth-order polynomial. The combination of these equations gives the function g(i) (here, all inclinations are given in radians):

(3) \begin{align} g(i) = & \Bigg\{ \begin{array}{ll} -0.4 \cos{i} + 1.4 & \ \ i < \pi / 3\\ -0.1773 i^4 + 1.1211 i^3 - 1.9149 i^2 +0.5022 i + 1.6222 & \ \ i \geq \pi/3 \end{array} \end{align}

Given these two independent fits, it is necessary to account for the mutual interaction between inclination and mass ratio. The factor $h(q_\mathrm{out}, i)$ reduces stability at high inclinations for roughly equal-mass systems (Equation (4)). Thus, the empirical fit to the stability is given by Equation (1) with

(4) \begin{equation} \log h(q_\mathrm{out},i) = \frac{-iq_\mathrm{out}^{1.3}}{1500} \ . \end{equation}

4. Comparison to other stability boundary fits

Here, we compare our stability boundary with other works. We first summarise the fractions of orbits correctly classified as stable/unstable, then directly compare to the recent algebraic fit from V+22 and discuss the performance of different instability thresholds and their regions of validity.

4.1 Overall fractions of correctly classified orbits

Equation (1) was compared to other well-known hierarchical stability boundary models by predicting stability for $10^4$ orbits initialised according to Table 2. In order to focus on the regime of interest and better differentiate between the quality of the different fits, $a_\mathrm{in}/R_\mathrm{p, out}$ was restricted to be between 0.5 and 1.5 times the stability boundary as described by Equation (1). This results in less impressive classification scores than are found in other work (V+22), since we focus on the region in which it is most difficult to predict stability. As in Section 2, we classified an orbit as stable if the distance between $m_1$ and $m_2$ remained smaller than $r_\mathrm{p, out}/2$ for 100 outer orbital periods. Equation (1) predicted stability correctly for 87.7% of orbits, compared to 67.4% for the fit proposed by MA01 and 70.5% for the adjusted fit put forward in V+22. Of the misclassified orbits, 53.3% were unstable (classified stable) and the rest were stable, but classified as unstable, compared to 97.6% and 98.6% of misclassified orbits being stable (classified as unstable) for the MA01 and V+22 fits, respectively. This stark discrepancy arises from the inconsistency of the 2/5 power law with the 1/3 power law Hill regime in mass-hierarchical systems, in which both the MA01 and V+22 fits drop off too quickly with $q_\mathrm{out}$ , causing them to over-classify orbits as unstable.

4.2 Varying the instability threshold

The definition of stability used in this paper differs somewhat from that put forward by V+22, who extend instability to include any system for which $a_\mathrm{in}$ or $a_\mathrm{out}$ change by more than 10% after 100 $P_\mathrm{out}$ , in order to catch systems which have not yet been disrupted but will be in the near future. Figure 4 compares the stability boundaries put forward by this work, MA01 and V+22 against inclination, for different mass ratios and different stability criteria.

It is clear that Equation (1) most closely approximates the stability boundary when $q_\mathrm{out} \ll 1$ , due to its adherence to the Hill mass-ratio dependence in this regime. Since the MA01 and V+22 stability boundaries both approach a 2/5 power law at small $q_\mathrm{out}$ , they underestimate the maximum stable $a_\mathrm{in}$ in these systems. At $q_\mathrm{out} = 0.1$ , the stability boundary is best fit by the V+22 formula.

Table 2. Range of orbital elements used for fit evaluation. Inclination followed an isotropic distribution, mass ratios were distributed logarithmically, and all other elements were distributed uniformly. $a_{\mathrm{in, crit}}$ refers to the critical $a_\mathrm{in}$ of the stability boundary predicted by Equation (1).

Figure 4. Fraction of stable orbits out of 20 evenly spaced mean anomalies against i, plotted for $q_\mathrm{out}=0.1$ (top row) and $0.001$ (bottom row). Left panels: our threshold for instability ( $a_\mathrm{in}>a_\mathrm{out}/2$ ). Right panels: V+22’s threshold for instability (either separation varies by $>\!\!10\%$ ).

As for the different definitions of stability used, they seem to be largely consistent below $\sim\!120^{\circ}$ inclination, however in retrograde systems, the V+22 definition provides a much broader region of ambiguous stability, with a fraction of orbits being labeled stable depending on the initial mean anomaly. This is probably due to the large chaotic extent of meta-stable orbits that are unstable according to the V+22 criterion, but are stable for our criterion. This retrograde region is accurately fit by the V+22 stability boundary in regimes of moderate $q_\mathrm{out}$ , however that fit underestimates the high inclination stability boundary when using our definition.

For $10^{-6} < q_\mathrm{out} < 1$ , our definition of stability appears to provide a more clearly defined, smooth, boundary. However, for $q_\mathrm{out} > 1$ , the V+22 definition is more appropriate, as it considers changes to both $a_\mathrm{in}$ and $a_\mathrm{out}$ . In these systems, instability can mean ejection or disruption of $m_\mathrm{out}$ , which may not be detected by the definition of stability we use.

The stability boundary described by Equation (1) unifies our previous understandings of the Hill regime (Grishin et al. 2017) and the similar-mass regime (MA01). It provides an accurate prediction of stability for any system with $q_\mathrm{out} \lesssim 1$ , and captures the complexity of the inclination dependence. However, the dependence of the $a_\mathrm{in}/a_\mathrm{out}$ SB cannot be extrapolated far beyond $q_\mathrm{out} = 1$ , since $a_\mathrm{in}$ and $a_\mathrm{out}$ become comparable at the SB, and the system ceases to be hierarchical. Hence, in that regime it is necessary to adopt a fit that asymptotes toward independence from $q_\mathrm{ out}$ (V+22; MA01).

5. Conclusions

We have obtained an empirical fit to hierarchical three-body stability for systems with $q_\mathrm{out} \lesssim 1$ . This fit attempts to bridge the gap between the well understood hierarchical-in-mass Hill regime (Grishin et al. 2017) and previous understandings of the equal-mass regime (MA01; V+22). We have mapped the stability boundary extensively across mass ratio and inclination within this range, and determined a ’turn-off’ point from the Hill regime at $q_\mathrm{out} \approx 0.1$ . We have also confirmed the accuracy of our fits through an extensive test, correctly predicting stability in 87.7% of systems on or near the SB.

Our empirical fit captures the features of the non-monotonic inclination dependence of the stability boundary more accurately than other contemporary works which derive semi-analytic fits (Mardling & Aarseth 2001; Grishin et al. 2017; Vynatheya et al. 2022). Furthermore, it is the only stability boundary fit which unifies the resonant equal-mass regime and the Hill regime, dominated by tidal shearing. This is reflected in its classification performance, where it out-performs alternative algebraic stability boundary models for systems with $q_\mathrm{ out} \ll 1$ while maintaining a comparable accuracy for those with $q_\mathrm{out} \approx 1$ .

The stability boundary fit presented in this paper is limited by its treatment of $q_\mathrm{out}$ as it approaches 1. Equation (2) defines a stability boundary which increases indefinitely with $q_\mathrm{out}$ , and with an increasing gradient. This arises from a somewhat myopic analysis of the behaviour of $q_\mathrm{out}$ in co-planar prograde systems, a relationship which is not borne out in more inclined systems (see Figure 2), and cannot be sustained as $q_\mathrm{out}$ exceeds 1. Furthermore, this work deals only in passing with eccentricity as a factor impacting stability, and does not attempt to consider its interaction with inclination, nor to model its effect beyond a simple linear fit. This approximation is accurate for $q_\mathrm{out} < 0.1$ , but fails in the equal-mass regime and beyond.

Equation (1) provides a promising foundation for a generalised stability boundary for all hierarchical triples. Future work should consider the fits for systems with $q_\mathrm{out} > 1$ provided by MA01, V+22 and others, and expand our study of the parameter space to evaluate the interactions of eccentricities with other parameters, in order to construct this unified fit. Analytic future work may also reveal the exact transition between the resonant and Hill regimes, to be compared with our empirical findings.

The study of populations of triples and their evolution can be numerical (Anderson, Lai, & Storch 2017; Toonen et al. 2020; Hamers et al. 2021; Grishin & Perets 2022), or semi-analytical (Muñoz et al. 2016; Mangipudi et al. 2022). However, even numerical population studies generally require secular approximations to reduce computational costs, and therefore rely on dynamical stability checks during the evolution. Our new stability boundary fit will aid in distinguishing stable and unstable systems, which will most likely affect the overall branching ratios of the different outcomes.

Recently, machine learning tools have also been utilised to study triple stability (Lalande & Trani 2022; Vynatheya et al. 2022). Due to their complexity, large neural networks have the potential to capture non-linear multi-dimensional interactions between parameters more accurately than any algebraic fit to the stability boundary, and further study should endeavour to apply these methods to a wider range of systems for more broadly applicable classification. However, neural networks could be harder to implement in simple population-synthesis style studies, where our simple algebraic stability fit will be useful.