1. Introduction
The three-body problem had been studied for over 100 years (Poincaré & Magini Reference Poincaré and Magini1899) and is notoriously chaotic, with no analytic solution (Poincaré Reference 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 Reference Monaghan1976a,b). Only recently, the explicit dependence on the orbital elements was found with various methods (e.g. Stone & Leigh Reference Stone and Leigh2019; Ginat & Perets Reference Ginat and Perets2021a,b; Kol Reference Kol2021a,b; Manwadkar et al. Reference Manwadkar, Kol, Trani and Leigh2021).
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 Reference Harrington1968; Ford, Kozinsky, & Rasio Reference Ford, Kozinsky and Rasio2000; Naoz et al. Reference Naoz, Farr, Lithwick, Rasio and Teyssandier2013; Toonen et al. Reference Toonen, Portegies Zwart, Hamers and Bandopadhyay2020).
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 Reference Duchêne and Kraus2013; Moe & Di Stefano Reference Moe and Di Stefano2017), 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. Reference Doyle2011), moons and binary asteroids (Richardson & Walsh Reference Richardson and Walsh2006; Naoz, Perets, & Ragozzine Reference Naoz, Perets and Ragozzine2010), or binary stars/black holes around an intermediate or super-massive black hole (Fragione et al. Reference Fragione, Grishin, Leigh, Perets and Perna2019).
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 Reference Eggleton and Kiseleva1995; Holman & Wiegert Reference Holman and Wiegert1999; Petrovich Reference Petrovich2015; He & Petrovich Reference He and Petrovich2018; Quarles et al. Reference Quarles, Satyal, Kostov, Kaib and Haghighipour2018), while new methods explore the boundary using a machine learning approach for specific choices of initial conditions (Lalande & Trani Reference Lalande and Trani2022; Vynatheya et al. Reference Vynatheya, Hamers, Mardling and Bellinger2022, 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 Reference Hill1878). 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 Reference Henon1970; Hamilton & Burns Reference Hamilton and Burns1991). The inclination dependence of the SB was examined analytically by Innanen (Reference Innanen1979, Reference Innanen1980), Mylläri et al. (Reference Mylläri, Valtonen, Pasechnik and Mikkola2018), 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 Reference Lidov1962; Kozai Reference Kozai1962). This reduces the critical stability radius for mutually inclined orbits. Grishin et al. (Reference Grishin, Perets, Zenati and Michaely2017) 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. Reference Carruba, Burns, Nicholson and Gladman2002).
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 Reference Mardling2008, Reference Mardling2013). 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 Reference Chirikov1979) and can be used to estimate stability. Mardling & Aarseth (Reference Mardling and Aarseth2001, 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. Reference Toonen, Portegies Zwart, Hamers and Bandopadhyay2020; Hamers et al. Reference Hamers, Rantala, Neunteufel, Preece and Vynatheya2021; Grishin & Perets Reference Grishin and Perets2022). 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. Reference Reggiani2022). 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 Reference Rein and Spiegel2014). 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. Reference Bhaskar, Li, Hadden, Payne and Holman2021, 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. (Reference Grishin, Perets, Zenati and Michaely2017) 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. (Reference Grishin, Perets, Zenati and Michaely2017) 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 (Reference Hayashi, Trani and Suto2022) 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.
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 Reference Hamilton and Burns1991; Grishin et al. Reference Grishin, Perets, Zenati and Michaely2017). 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. (Reference Grishin, Perets, Zenati and Michaely2017), was later derived analytically by Grishin, Perets, & Fragione (Reference Grishin, Perets and Fragione2018) for small
$q_\mathrm{out}$
, and was recently extended to any mass ratio by Mangipudi et al. (Reference Mangipudi, Grishin, Trani and Mandel2022). 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).
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. Reference Grishin, Perets, Zenati and Michaely2017; Mangipudi et al. Reference Mangipudi, Grishin, Trani and Mandel2022).
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).
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:
\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})$
:
\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):
\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
\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.
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).
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. Reference Grishin, Perets, Zenati and Michaely2017) 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. Reference Grishin, Perets, Zenati and Michaely2017) 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 Reference Mardling and Aarseth2001; Grishin et al. Reference Grishin, Perets, Zenati and Michaely2017; Vynatheya et al. Reference Vynatheya, Hamers, Mardling and Bellinger2022). 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 Reference Anderson, Lai and Storch2017; Toonen et al. Reference Toonen, Portegies Zwart, Hamers and Bandopadhyay2020; Hamers et al. Reference Hamers, Rantala, Neunteufel, Preece and Vynatheya2021; Grishin & Perets Reference Grishin and Perets2022), or semi-analytical (Muñoz et al. Reference Muñoz, Lai and Liu2016; Mangipudi et al. Reference Mangipudi, Grishin, Trani and Mandel2022). 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 Reference Lalande and Trani2022; Vynatheya et al. Reference Vynatheya, Hamers, Mardling and Bellinger2022). 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.
Acknowledgements
We thank Pavan Vynatheya for useful discussions, Rosemary Mardling for helpful discussions and comments on the manuscript, and the anonymous referee for useful feedback. E. G. and I. M. acknowledge support from the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE17010004. I. M. is a recipient of the Australian Research Council Future Fellowship FT190100574. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. The participation of I. M. at the Aspen Center for Physics was partially supported by the Simons Foundation. This work made use of the OzSTAR high performance computer at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS).
Open access
