# Empirical Stability Boundary for Hierarchical Triples

Max Tory  School of Physics and Astronomy, Monash University, VIC 3800,
Australia  Evgeni Grishin  School of Physics and Astronomy, Monash University,
VIC 3800, Australia  [ Ilya Mandel  School of Physics and Astronomy, Monash
University, VIC 3800, Australia

(dd Mmm YYYY; dd Mmm YYYY; dd Mmm YYYY)

###### 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 un-
mapped, 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_{\rm out}$), mutual
inclination ($i$), initial mean anomaly and eccentricity of both the inner and
outer binaries ($e_{\rm in}$ and $e_{\rm out}$ respectively). We fit the
dependence of the stability boundary on $q_{\rm out}$ as a threshold on the
ratio of the inner binary’s semi-major axis to the outer binary’s pericentre
separation $a_{\rm in}/R_{\rm p,out}\leq 10^{-0.6+0.04q_{\rm out}}q_{\rm
out}^{0.32+0.1q_{\rm 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.

###### doi:

10.1017/pasa.2022.xx

###### keywords:

Hierarchical systems - Stability - Gravitation - Celestial Mechanics

OzGrav: Australian Research Council Centre of Excellence for Gravitational
Wave Discovery, Clayton, VIC 3800, Australia E.
<EMAIL_ADDRESS>OzGrav: Australian Research Council Centre
of Excellence for Gravitational Wave Discovery, Clayton, VIC 3800, Australia
30 August 2022

##  1 Introduction

The three-body problem had been studied for over one hundred 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, 2021b, a; Kol, 2021b, a; 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 et al., 2000; Naoz et al.,
2013; Toonen et al., 2020a).

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 et al., 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_{\rm in}$, and an outer
separation $a_{\rm out}$, the stability boundary (hereafter SB) will generally
depend on the separation ratio $\alpha\equiv a_{\rm in}/a_{\rm 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_{\rm out}\equiv m_{\rm in}/m_{\rm out}$, and the mutual
inclination between the two orbits, $i$. Here $m_{\rm in}$ is the total mass
of the inner binary and $m_{\rm out}$ is the mass of the outer companion.

For extreme mass ratios $q_{\rm out}\ll 1$, the Hill radius $r_{\rm H}\equiv
a_{\rm out}\left(q_{\rm out}/3\right)^{1/3}\ll a_{\rm out}$ defines the length
scale of the SB (Hill, 1878). Tight inner binary orbits with $a_{\rm in}\ll
r_{\rm H}$ will be stable while wide orbits with $a_{\rm in}\gtrsim r_{\rm 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).

For comparable masses, $q_{\rm 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_{\rm
out})^{-2/5}\propto(q_{\rm out}/(1+q_{\rm out}))^{2/5}$. Although strictly
valid for comparable masses, $q_{\rm 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., 2020b; Hamers et
al., 2021; Grishin & Perets, 2022). Finally, recent observations find
companions with mass as low as $0.2M_{\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_{\rm 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 7. 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_{\rm p,out}/2\equiv a_{\rm out}(1-e_{\rm
out})/2$. In orbits with $q_{\rm out}\leq 1$, the inner binary can be
disrupted before its separation exceeds $a_{\rm 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 sec. 4.2. Systems with large values of $q_{\rm 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_{\rm out}$) without disruption, following previous work. Grishin et al.
(2017) integrated for $100P_{\rm out}$ to incorporate the impact of secular
changes that occur on timescales of $\lesssim 10P_{\rm 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_{\rm out}$ and obtained similar
results. (V+22) also used a simulation duration of $100P_{\rm out}$, arguing
that $99\%$ of the unstable orbits are classified as such within $100P_{\rm
out}$ (see their Fig. 1). However, the actual timescale for instability may be
sensitive to initial conditions. While this paper was under review, Hayashi et
al. (2022) integrated a limited sample of initial conditions for much longer
timescales to show that some orbits (especially orbits with large $e_{\rm
out}$ and large mutual inclinations) could reach instability much later, after
$\gtrsim 10^{4}P_{\rm out}$.

###  2.1 Initialisation

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_{\rm in}$ on the ordinate of each plot is scaled as appropriate to the figure. Variable | Fig. 1 | Fig. 2 | Fig. 3 | Fig. 4
---|---|---|---|---
$i$ | [0,180] | 0, 150 | 0 | [0,180]
$e_{\rm in}$ | 0 | 0 | [0, 0.85] | 0
$e_{\rm out}$ | 0 | 0 | [0, 0.85] | 0
$q_{\rm in}$ | $10^{-2}$ | $10^{-2}$ | $10^{-2}$ | [$10^{-2}$, $10^{0}$]
$q_{\rm out}$ | 1, $10^{-5}$ | [$10^{-4}$, $10^{0}$] | $10^{-2}$ |  $10^{-1}$, $10^{-3}$
$\\#$ of runs/$10^{5}$ | $2\times 4.68$ | $2\times 12$ | $2\times 1.36$ | $4\times 4.68$

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

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_{\rm in}$ and $e_{\rm in}$. Finally, $m_{\rm out}$ was
initialised in orbit around the centre-of-mass of $m_{1}$ and $m_{2}$, based
on $a_{\rm out}$, $e_{\rm 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_{\rm 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.

##  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 2: Fraction of stable orbits out of 20 evenly spaced mean anomalies
against the mass ratio $q_{\rm out}$, for inclinations $i=0$ (left) and
$i=150^{\circ}$ (right). The white line indicates the Hill radius $r_{\rm H}$
in both panels. At extreme mass ratios, the relationship is consistent with
the Hill regime, while at $q_{\rm 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).

Fig. 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_{\rm out}\ll 1$
Hill approximation up to $\log[q_{\rm out}]\leq-1$, with prograde orbits
stable up to around $0.5r_{\rm H}$, and retrograde orbits as far as $0.8r_{\rm
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 Fig. 1 shows the $q_{\rm 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_{\rm in}/a_{\rm 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_{\rm
out}=1$ case, while the decay is more abrupt in extreme $q_{\rm 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 et al. (2018) for small $q_{\rm 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_{\rm out}$. In
the Hill approximation, where $\log q_{\rm 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_{\rm 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).

###  3.3 Inner and Outer Eccentricities

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

Previous work has demonstrated that in systems with $q_{\rm out}\sim 1$, the
stability boundary of hierarchical triples depends on both outer and inner
eccentricity (; ). 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_{\rm in}$ and the Hill pericentre radius $r_{\rm H,p}\equiv
r_{\rm H}(1-e_{\rm out})$, rather than $a_{\rm in}/r_{\rm 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_{\rm 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).

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

###  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_{\rm in}$ is mostly irrelevant and the effects
of $e_{\rm 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:

$\frac{a_{\rm in}}{R_{\rm p,out}}=f(q_{\rm out})\cdot g(i)\cdot h(q_{\rm
out},i).$ (1)

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_{\rm
out})$:

$f(q_{\rm out})=10^{-0.6+0.04q_{\rm out}}q_{\rm out}^{0.32+0.1q_{\rm out}},$
(2)

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):

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

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

$\log h(q_{\rm out},i)=\frac{-iq_{\rm out}^{1.3}}{1500}\ .$ (7)

##  4 Comparison to other stability boundary fits

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

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_{\rm in}/R_{\rm 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_{\rm 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_{\rm 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_{\rm in}$ or $a_{\rm out}$ change by more than 10% after 100 $P_{\rm 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_{\rm 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_{\rm out}$, they underestimate the
maximum stable $a_{\rm in}$ in these systems. At $q_{\rm out}=0.1$, the
stability boundary is best fit by the V+22 formula.

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_{\rm out}$, however that fit
underestimates the high inclination stability boundary when using our
definition.

For $10^{-6}<q_{\rm out}<1$, our definition of stability appears to provide a
more clearly defined, smooth, boundary. However, for $q_{\rm out}>1$, the V+22
definition is more appropriate, as it considers changes to both $a_{\rm in}$
and $a_{\rm out}$. In these systems, instability can mean ejection or
disruption of $m_{\rm 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_{\rm out}\lesssim 1$, and captures the complexity of the inclination
dependence. However, the dependence of the $a_{\rm in}/a_{\rm out}$ SB cannot
be extrapolated far beyond $q_{\rm out}=1$, since $a_{\rm in}$ and $a_{\rm
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_{\rm out}$ (; ).

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_{\rm in,crit}$ refers to the critical $a_{\rm in}$ of the stability boundary predicted by Equation 1. $a_{\rm in}/a_{\rm in,crit}$ | $\log q_{\rm in}$ | $\log q_{\rm out}$ | $i$ | $e_{\rm in}$ | $e_{\rm out}$
---|---|---|---|---|---
[0.5,1.5] | [-2,0] | [-6,0] | [0,$\pi$] | [0,1] | [0,1]

##  5 Conclusions

We have obtained an empirical fit to hierarchical three-body stability for
systems with $q_{\rm 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 (; ). 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_{\rm 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_{\rm out}\ll 1$ while maintaining a comparable
accuracy for those with $q_{\rm out}\approx 1$.

The stability boundary fit presented in this paper is limited by its treatment
of $q_{\rm out}$ as it approaches 1. Equation 2 defines a stability boundary
which increases indefinitely with $q_{\rm out}$, and with an increasing
gradient. This arises from a somewhat myopic analysis of the behaviour of
$q_{\rm out}$ in co-planar prograde systems, a relationship which is not borne
out in more inclined systems (see Fig. 2), and cannot be sustained as $q_{\rm
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_{\rm 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_{\rm 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 et al., 2017; Toonen et al., 2020b; 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.

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

## References

  * Anderson et al. (2017) Anderson, K. R., Lai, D., & Storch, N. I. 2017, MNRAS, 467, 3066
  * Bhaskar et al. (2021) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2021, AJ, 161, 48
  * Carruba et al. (2002) Carruba, V., Burns, J. A., Nicholson, P. D., & Gladman, B. J. 2002, Icarus, 158, 434
  * Chirikov (1979) Chirikov, B. V. 1979, Phys. Rep., 52, 263
  * Doyle et al. (2011) Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602
  * Duchêne & Kraus (2013) Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269
  * Eggleton & Kiseleva (1995) Eggleton, P., & Kiseleva, L. 1995, ApJ, 455, 640
  * Ford et al. (2000) Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385
  * Fragione et al. (2019) Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47
  * Ginat & Perets (2021a) Ginat, Y. B., & Perets, H. B. 2021a, Physical Review X, 11, 031020
  * Ginat & Perets (2021b) —. 2021b, arXiv e-prints, arXiv:2108.01085
  * Grishin & Perets (2022) Grishin, E., & Perets, H. B. 2022, MNRAS, 512, 4993
  * Grishin et al. (2018) Grishin, E., Perets, H. B., & Fragione, G. 2018, MNRAS, 481, 4907
  * Grishin et al. (2017) Grishin, E., Perets, H. B., Zenati, Y., & Michaely, E. 2017, MNRAS, 466, 276
  * Hamers et al. (2021) Hamers, A. S., Rantala, A., Neunteufel, P., Preece, H., & Vynatheya, P. 2021, MNRAS, 502, 4479
  * Hamilton & Burns (1991) Hamilton, D. P., & Burns, J. A. 1991, Icarus, 92, 118
  * Harrington (1968) Harrington, R. S. 1968, AJ, 73, 190
  * Hayashi et al. (2022) Hayashi, T., Trani, A. A., & Suto, Y. 2022, arXiv e-prints, arXiv:2209.08487
  * He & Petrovich (2018) He, M. Y., & Petrovich, C. 2018, MNRAS, 474, 20
  * Henon (1970) Henon, M. 1970, A&A, 9, 24
  * Hill (1878) Hill, G. W. 1878, American Journal of Mathematics, 1, 5
  * Holman & Wiegert (1999) Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621
  * Innanen (1979) Innanen, K. A. 1979, AJ, 84, 960
  * Innanen (1980) —. 1980, AJ, 85, 81
  * Kol (2021a) Kol, B. 2021a, Celestial Mechanics and Dynamical Astronomy, 133, 17
  * Kol (2021b) —. 2021b, arXiv e-prints, arXiv:2107.12372
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Lalande & Trani (2022) Lalande, F., & Trani, A. A. 2022, arXiv e-prints, arXiv:2206.12402
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
  * Mangipudi et al. (2022) Mangipudi, A., Grishin, E., Trani, A. A., & Mandel, I. 2022, ApJ, 934, 44
  * Manwadkar et al. (2021) Manwadkar, V., Kol, B., Trani, A. A., & Leigh, N. W. C. 2021, MNRAS, 506, 692
  * Mardling (2008) Mardling, R. A. 2008, in The Cambridge N-Body Lectures, ed. S. J. Aarseth, C. A. Tout, & R. A. Mardling, Vol. 760, 59
  * Mardling (2013) —. 2013, MNRAS, 435, 2187
  * Mardling & Aarseth (2001) Mardling, R. A., & Aarseth, S. J. 2001, MNRAS, 321, 398
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15
  * Monaghan (1976a) Monaghan, J. J. 1976a, MNRAS, 176, 63
  * Monaghan (1976b) —. 1976b, MNRAS, 177, 583
  * Muñoz et al. (2016) Muñoz, D. J., Lai, D., & Liu, B. 2016, MNRAS, 460, 1086
  * Mylläri et al. (2018) Mylläri, A., Valtonen, M., Pasechnik, A., & Mikkola, S. 2018, MNRAS, 476, 830
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155
  * Naoz et al. (2010) Naoz, S., Perets, H. B., & Ragozzine, D. 2010, ApJ, 719, 1775
  * Petrovich (2015) Petrovich, C. 2015, ApJ, 808, 120
  * Poincaré (1892) Poincaré, H. 1892, Les méthodes nouvelles de la mécanique céleste, doi:10.3931/e-rara-421
  * Poincaré & Magini (1899) Poincaré, H., & Magini, R. 1899, Il Nuovo Cimento, 10, 128
  * Quarles et al. (2018) Quarles, B., Satyal, S., Kostov, V., Kaib, N., & Haghighipour, N. 2018, ApJ, 856, 150
  * Reggiani et al. (2022) Reggiani, M., Rainot, A., Sana, H., et al. 2022, A&A, 660, A122
  * Rein & Spiegel (2014) Rein, H., & Spiegel, D. S. 2014, MNRAS, 446, 1424
  * Richardson & Walsh (2006) Richardson, D. C., & Walsh, K. J. 2006, Annual Review of Earth and Planetary Sciences, 34, 47
  * Stone & Leigh (2019) Stone, N. C., & Leigh, N. W. C. 2019, Nature, 576, 406
  * Toonen et al. (2020a) Toonen, S., Portegies Zwart, S., Hamers, A. S., & Bandopadhyay, D. 2020a, A&A, 640, A16
  * Toonen et al. (2020b) —. 2020b, A&A, 640, A16
  * Vynatheya et al. (2022) Vynatheya, P., Hamers, A. S., Mardling, R. A., & Bellinger, E. P. 2022, arXiv e-prints, arXiv:2207.03151