#  A Stability Timescale for Non-Hierarchical Three-Body Systems

Eric Zhang Department of Physics and Astronomy, University of California, Los
Angeles, CA 90095, USA  Mani L. Bhaumik Institute for Theoretical Physics,
Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, USA
Department of Physics and Astronomy, University of California, Riverside, CA,
92507, USA  Smadar Naoz Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA  Mani L. Bhaumik Institute for
Theoretical Physics, Department of Physics and Astronomy, UCLA, Los Angeles,
CA 90095, USA  Clifford M. Will Department of Physics, University of Florida,
Gainesville, FL 32611, USA  Institut d’Astrophysique, Sorbonne Université,
75014 Paris, France

###### Abstract

The gravitational three-body problem is a fundamental problem in physics and
has significant applications to astronomy. Three-body configurations are often
considered stable as long the system is hierarchical; that is, the two orbital
distances are well-separated. However, instability, which is often associated
with significant energy exchange between orbits, takes time to develop.
Assuming two massive objects in a circular orbit and a test particle in an
eccentric orbit, we develop an analytical formula estimating the time it takes
for the test particle’s orbital energy to change by an order of itself. We
show its consistency with results from N-body simulations. For eccentric
orbits in particular, the instability is primarily driven not by close
encounters of the test particle with one of the other bodies, but by the
fundamental susceptibility of eccentric orbits to exchange energy at their
periapsis. Motivated by recent suggestions that the galactic center may host
an intermediate-mass black hole (IMBH) as a companion to the massive black
hole Sgr A*, we use our timescale to explore the parameter space that could
harbor an IMBH for the lifetime of the S-cluster of stars surrounding Sgr A*.
Furthermore, we show that the orbit of an S-star can be stable for long
timescales in the presence of other orbital crossing stars, thus suggesting
that the S-cluster may be stable for the lifetimes of its member stars.

Exoplanets (498), Black holes (162), Three-body problem (1695), Supermassive
black holes (1663), Intermediate-mass black holes (816), Galactic center
(565), Galaxies (573), High energy astrophysics (739)

##  1 Introduction

The stability of triple-body systems is a pervasive problem in astrophysics.
The three-body problem describes the dynamics of systems ranging from
planetary systems to the orbits of stars and compact objects. A “stable” bound
three-body system is loosely defined as one in which the energy of each orbit
stays roughly the same, and “instability” is associated with systems in which
the two orbits exchange energy. As a dramatic example, a bound orbit in a
three-body system becomes unbound if its energy changes from negative to
positive, and such systems are said to be unstable. Another example for
instability can be considered an exchange process between the two orbits,
which also dramatically changes the energy of each orbit.

The general approach in the literature often relies on a criterion by which
the system could be considered stable “in the long run”. A common approach to
developing such criteria is to designate hierarchical systems, or those in
which the two orbital distances are well-separated, as long-term stable. The
question of stability is then equivalent to determining a critical distance
between the two orbits, where the system switches from hierarchical to non-
hierarchical.

The most straightforward such criterion involves the Hills mechanism (Hills,
1988). For a system composed of a binary and a tertiary, the Hills critical
distance is the separation between the tertiary and the inner binary at which
the tertiary’s gravitational potential exceeds the primary’s gravitational
potential. In this case, the secondary body in the inner binary may hop
between the primary and the tertiary, yielding a new configuration of a
tertiary-secondary and a primary. In the co-planar case, the functional
expression of this place represents the first Lagrange point, which indicates
the position where the gradient of the potential in the rotating frame is zero
(e.g., Murray & Dermott, 2000a; Binney & Tremaine, 2008). A stability
criterion may be obtained by requiring that the tertiary is never closer to
either of the other bodies than the Hills critical distance. The criterion has
been improved by Hamilton & Burns (1991) and Grishin et al. (2016).

The stability of three massive objects has a generalized, hierarchy-based,
stability condition often used in the literature (e.g., Mardling & Aarseth,
2001), and similar stable-unstable boundaries were derived by Eggleton &
Kiseleva (1995); Petrovich (2015); Tory et al. (2022); Vynatheya et al.
(2022); Hayashi et al. (2022). Considering hierarchical systems, where one
mass orbits on a tight configuration about the primary and a tertiary is on a
wider orbit, a condition is often used to estimate the long-term stability
against high-eccentricity excitations due to secular dynamics (e.g., Ivanov et
al., 2005; Lithwick & Naoz, 2011; Katz & Dong, 2012; Naoz & Silk, 2014;
Antonini et al., 2014; Bode & Wegg, 2014) and non-secular perturbations to
secular dynamics (e.g., Antognini et al., 2014; Luo et al., 2016; Grishin et
al., 2018; Bhaskar et al., 2020).

However, instability is a time-sensitive concept, and not every non-
hierarchical system is instantaneously unstable. Mylläri et al. (2018) noted
the dependence of stability on time, albeit without an explicit timescale.
Additionally, Mushkin & Katz (2020) developed a stability timescale for the
outer orbit in hierarchical systems, based on formulae for secular energy
exchange (e.g., Roy & Haddow, 2003). Further studies on the time-dependence of
stability were done by Hayashi et al. (2022, 2023), using N-body simulations
of mildly hierarchical triples.

In this paper, we develop an analytical stability timescale, expanding the
parameter space to non-hierarchical systems. We consider three-body systems
consisting of a primary and companion of masses $m_{p}$ and $m_{c}$,
respectively, on a circular orbit, and a test particle $m_{t}$
$(m_{p}>m_{c}\gg m_{t})$, on a highly eccentric orbit either about the primary
body or about the massive binary as a whole, such that the system is not
necessarily hierarchical. We denote these two cases the “external companion”
and “internal companion” cases, respectively. Interestingly, the stability of
such systems is primarily sensitive to a single parameter, $\alpha$, defined
as the ratio of the periapsis distance of the test particle to the companion’s
semimajor axis:

$\alpha=\frac{r_{\rm peri}}{a_{c}}=\frac{a_{t}(1-e_{t})}{a_{c}}\ ,$ (1)

where $a_{t}$ and $e_{t}$ are the semimajor axis and eccentricity of the test
particle’s orbit. The system is most unstable for $\alpha=1$, when the
companion orbits at the same distance as the test particle’s periapsis.

We derive our stability timescale based on how long it takes for the test
particle’s orbital energy to change. We also show that the changes in energy
are driven by the fundamental susceptibility of the eccentric test particle to
energy changes at its periapsis, and not by close encounters with the
companion or the primary. Notably, we do not use the secular approximation, in
which the phases of each orbit are averaged over timescales much longer than
each orbital period, or any perturbations to the secular approximation. In
that approximation, often used to describe the evolution of hierarchical
systems, the long-term change in the energy of each orbit is zero. Thus, in
order to explore the instability associated with orbital energy exchange,
which is a feature of non-hierarchical systems, we must describe systems by
other means.

Figure 1:  An illustration of the possible configurations. Left Panel: A
system where the companion orbits external to the test particle’s periapsis.
The central body is the primary, so that the test particle’s orbital elements
are defined via the position vector ${\bf r}={\bf r}_{tp}$ and the velocity
vector ${\bf v}={\bf v}_{tp}$. Right Panel: An internal companion. The
effective central “body” is at the center of mass of the primary and the
companion, and is fixed at the origin. The test particle’s orbit is calculated
using ${\bf r}={\bf r}_{t}$ and ${\bf v}={\bf v}_{t}$.

In particular we are motivated by the configuration at the center of our
galaxy, which includes many crossing orbits. As an example, we focus on the
hypothetical existence of an intermediate mass black hole (IMBH) that may
orbit around the supermassive black hole (SMBH). The existence of such an IMBH
has been investigated based on a combination of theoretical and observational
studies (e.g., Hansen & Milosavljević, 2003; Maillard et al., 2004; Gürkan &
Rasio, 2005; Gualandris & Merritt, 2009; Chen & Liu, 2013; Generozov &
Madigan, 2020; Fragione et al., 2020; Zheng et al., 2020; Naoz et al., 2020;
GRAVITY Collaboration et al., 2020; Rose et al., 2022). The IMBH may have a
crossing orbit with the well studied star, S0-2. This star is located close to
the SMBH Sgr A*. It has an orbital period of $16$ years and an eccentricity of
about $0.88$. The recent closest approach of this star has been used to test
and confirm the prediction of general relativity (GR) for the relativistic
redshift (e.g., GRAVITY Collaboration et al., 2018; Do et al., 2019) and the
advance of the periapsis (GRAVITY Collaboration et al., 2020). Thus, in
principle, S0-2 may be used to constrain the possibility of the existence of
such an IMBH. In addition, many of the S-star cluster orbits Sgr A* at similar
distances, with many potential orbital crossings (e.g., Ghez et al., 2005a;
Gillessen et al., 2009; Yelda et al., 2014). One may wonder how unstable these
orbits truly are, and if the S-star cluster is in fact a stable configuration.

The paper is organized as follows: we develop in Section 2 an analytic
timescale $T_{\rm stab}$ for the time it takes for the energy of the test
particle to change significantly. In Section 3 we then adopt, as a proof of
concept, the possible existence of an intermediate-mass black hole (IMBH) at
the center of our galaxy, which we treat numerically, and compare the results
with the analytic timescale. In Section 4 we discuss applications of our
timescale to various astrophysical settings. Finally, in Section 5 we
summarize our results.

##  2  Analytic Stability Timescale

In this section, we present an analytic derivation of the stability timescale.
For pedagogical purposes, we derive the timescale using two approaches: via
the test particle’s energy evolution, and via its semimajor axis evolution.
The two quantities are defined by the equation,

$E=\frac{1}{2}{\bf v}^{2}-\frac{GM}{r}=-\frac{GM}{2a_{t}}\ ,$ (2)

where $E$ is the (specific) energy111This definition of the energy includes
only the binding energy between the test particle and the central body, and
ignores the interaction energy with other bodies. This definition of the
orbital energy allows the second equality, involving the osculating semimajor
axis $a_{t}$, to hold. of the test particle’s orbit, $a_{t}$ is the test
particle’s semimajor axis, $M$ is the mass of the central body or system, $r$
is the distance between the test particle and the central body, and ${\bf v}$
is the velocity vector of the test particle relative to the central body. By
this relation, the energy and semimajor axis contain the same information, so
that the a change in one of the quantities by a given factor necessarily
accompanies a change in the other by the same factor.

Throughout the derivation, and the remainder of the paper, we adopt the
following notation: the masses of the primary, companion, and test particle
will be denoted by $m_{p}$, $m_{c}$, and $m_{t}$ respectively. The semimajor
axis, eccentricity, inclination, longitude of ascending node, argument of
periapsis, and true anomaly of both orbits will be denoted by $a$, $e$,
$\iota$, $\Omega$, $\omega$, and $f$, with a subscript $c$ indicating the
companion’s orbital elements, and the subscript $t$ indicating the test
particle’s orbital elements. The companion’s and test particle’s orbital
periods are denoted $P_{c}$ and $P_{t}$, respectively.

We fix the origin at the center of mass of all three bodies. Note that, as
$m_{t}\rightarrow 0$, this is equivalent to the center of mass of the primary
and companion. Position vectors of each of the bodies will be denoted by ${\bf
r}$, and velocity vectors will be denoted by ${\bf v}$, with subscripts $p$,
$c$, and $t$ for the primary, companion, and test particle. The relative
position and vectors between bodies are denoted by ${\bf r}_{ij}={\bf
r}_{i}-{\bf r}_{j}$, where $i$ and $j$ are the subscripts of the bodies. The
relative velocity vectors are similarly denoted by ${\bf v}_{ij}$. We will
write $r$ and $v$ for the magnitudes of these vectors.

As stated in the introduction, we assume a mass hierarchy between the three
bodies, $m_{p}>m_{c}\gg m_{t}$, with $m_{t}\rightarrow 0$ as it is a test
particle. We assume a highly eccentric222$e_{t}\gtrsim 0.5$; see Appendix A
for a discussion of the low-eccentricity regime. test particle orbit and a
circular companion orbit ($e_{c}$ = 0).

Because we are working in the non-hierarchical regime, where the orbital
distances are not well-separated, defining the orbits of the companion and the
test particle can be ambiguous. In particular, the mass of the central body,
$M$, defined in Eq. (2), depends on the configuration of the system. We define
the companion’s orbit to be external (internal) to the test particle’s orbit
if its semimajor axis333Since the companion’s orbit is circular, $a_{c}$ is
the same as the orbital distance, which is constant., $a_{c}$, is greater
(less) than the test particle’s periapsis distance, $r_{\rm
peri}=a_{t}(1-e_{t})$. If the companion is external, then the test particle’s
orbit is calculated with respect to the primary, and the central mass is
$M=m_{p}$. If the companion is internal, then the test particle’s orbit is
calculated with respect to the center of mass of the primary and the
companion, i.e., the origin, and the central mass is $M=m_{p}+m_{c}$. Fig. 1
shows the two possible configurations.

###  2.1  Energy Approach

Our derivation is based on the realization that when the test particle has a
highly eccentric orbit, its energy is especially susceptible to gravitational
perturbations near its periapsis. Accordingly, to evaluate the effect of these
encounters, we approximate the energy evolution as a series of discrete,
small, jumps, where the jumps occur at each periapsis passage. The energy may
either increase or decrease at each jump, depending upon the location of the
companion body at that time, and so we treat its evolution as a random walk.
We justify these assumptions analytically later in this section, and
numerically in Section 3 (specifically, see Fig. 2).

Indeed, it is possible for the test particle to undergo a large change in its
energy while it is not at its periapsis, such as in the case of a close
encounter between the test particle and the companion, e.g., if the test
particle enters within a few Hill radii of the companion. However, for the
purposes of our analytic derivation, we will ignore the effect of such events,
which are arguably much rarer than periapsis passages. In principle, this
assumption could cause our analytic timescale to overestimate the stability of
the system, but our numerical results in Section 3 indicate that this is, in
general, not the case.

Considering only small jumps in the energy, the energy evolves in a diffusive
manner, as previously assumed. Though the expected overall energy change after
$N$ such jumps averages to zero, the mean-square energy change is nonzero.
This behavior leads to a spread in the distribution of possible energies after
$N$ periapsis passages, so that for sufficiently large $N$, it is probable
that the energy has changed significantly, and thus the system has become
unstable. In particular, we say that the energy has changed significantly when
the standard deviation $\sigma$ of the distribution of energy changes $\Delta
E$ after $N$ periapsis passages is comparable to the initial energy itself –
in other words, when

$\sigma\left(\Delta E\right)=\sqrt{N}\sigma\left(\delta E\right)\sim E\ ,$ (3)

where $\sigma(\delta E)$ denotes the standard deviation of the distribution of
energy jumps $\delta E$ over a single periapsis passage.

The change in energy over any short time interval is given by

$\delta E=\frac{\partial E}{\partial{\bf v}}\cdot{\delta{\bf v}}={\bf
v}\cdot{\bf f}_{\rm pert}\,\delta t\ ,$ (4)

where ${\bf f}_{\rm pert}$ is the perturbing force (per unit test mass) due to
the companion’s presence. From this equation, we see that the energy change
over any short time interval depends only on two quantities; the test
particle’s speed (or magnitude of its velocity), and the component of the
perturbing force directed along the velocity vector, $f_{\rm
pert}^{\parallel}$. For highly eccentric orbits, the speed will be much
greater near periapsis than at other times in the orbit. On the other hand,
the force component $f_{\rm pert}^{\parallel}$ depends strongly on the
relative position of the test particle and companion, which varies quasi-
randomly so that $f_{\rm pert}^{\parallel}$ does not peak reliably. In
particular, the parallel component of the force takes on essentially random
values each time the test particle approaches its periapsis. Thus, our initial
assumptions are justified; jumps primarily drive the change in energy at
periapsis, and such jumps will behave like a random walk. The standard
deviation of the distribution of energy changes is estimated by the scale of
these changes, i.e.,

$\sigma\left(\delta E\right)\sim v_{\rm peri}\left|{\bf f}_{\rm
pert}\right|\delta t\ .$ (5)

We note that this depends only on the vector magnitudes, since the relative
orientation of the vectors controls the sign and the relative strength of each
energy jump, but not their overall scale.

Then Eq. (3) can be solved for $N$ to obtain the number of test particle
orbits before a significant energy change occurs. Noting that the stability
timescale $T_{\rm stab}$ is of order $NP_{t}$, we may find $T_{\rm stab}$ via
the expression,

$\frac{T_{\rm stab}}{P_{t}}\sim\frac{E^{2}}{v_{\rm peri}^{2}\left|{\bf f}_{\rm
pert}\right|^{2}\delta t^{2}}\ ,$ (6)

from which we may obtain $T_{\rm stab}$ in terms of known orbital parameters
by estimating each of the terms $E$, $v_{\rm peri}$, $\left|{\bf f}_{\rm
pert}\right|$, and $\delta t$ for various regimes.

The energy, $E$, is given by Eq. (2), and $v_{\rm peri}$ can be estimated by
the periapsis velocity of the test particle in a standard Keplerian orbit,
given by

$v_{\rm peri}\approx\sqrt{\frac{GM}{a_{t}}\frac{1+e_{t}}{1-e_{t}}}\ ,$ (7)

where, as in Eq. (2), $M=m_{p}$ for external companions, and $M=m_{p}+m_{c}$
for internal companions.

Regime | Companion Orbit | Companion Period | $M$ | $\left|{\bf f}_{\rm pert}\right|$ | $\delta t$
---|---|---|---|---|---
$r_{\rm peri}<a_{c}$, $T_{\rm peri}\ll P_{c}$ | External | Long | $m_{p}$ | $Gm_{c}/a_{c}^{2}$ | $T_{\rm peri}$
$r_{\rm peri}<a_{c}$, $T_{\rm peri}\gtrsim P_{c}$ | External | Short | $m_{p}$ | $Gm_{c}/a_{c}^{2}$ | $P_{c}/2\pi$
$r_{\rm peri}>a_{c}$, $T_{\rm peri}\gtrsim P_{c}$ | Internal | Short | $m_{p}+m_{c}$ | $Gm_{p}m_{c}a_{c}^{2}/Mr^{4}_{\rm peri}$ | $P_{c}/2\pi$
Table 1: A summary of estimates of $M$, $\left|{\bf f}_{\rm pert}\right|$, and
$\delta t$ for each regime. $\left|{\bf f}_{\rm pert}\right|$ is given by Eqs.
(9) and (11), and $\delta t$ is given by Eq. (17).

Depending on whether the companion is external or internal, there are two
possible expressions for the perturbing force, ${\bf f}_{\rm pert}$. When the
companion is external, the perturbing force is given by

${\bf f}_{\rm pert}=-\frac{Gm_{c}{\bf r}_{tc}}{r_{tc}^{3}}-\frac{Gm_{c}{\bf
r}_{cp}}{r_{cp}^{3}}\ ,$ (8)

where the first term is the gravitational force on the test particle from the
companion, and the second term is due to the non-inertial motion of the
primary about the primary-companion mutual center of mass. Under the
approximation that the majority of the encounters are weak, and noting that
$a_{c}>r_{\rm peri}$ for external companions, the distance between the two
bodies, $r_{tc}$, will typically be on the order of $a_{c}$. In particular,
$r_{tc}\sim r_{cp}\sim a_{c}$, so that the magnitude of the total perturbing
force can be estimated as

$\left|{\bf f}_{\rm pert}\right|\sim\frac{Gm_{c}}{a_{c}^{2}}\ .$ (9)

When the companion is internal, we calculate the test particle’s orbit with
respect to the center of mass of the primary and companion. In this case,
${\bf f}_{\rm pert}$ is

${\bf f}_{\rm pert}=-\frac{Gm_{p}{\bf r}_{tp}}{r_{tp}^{3}}-\frac{Gm_{c}{\bf
r}_{tc}}{r_{tc}^{3}}+\frac{G(m_{p}+m_{c}){\bf r}_{t}}{r_{t}^{3}}\ ,$ (10)

where the first two terms are the total gravitational force on the test
particle, and the third term subtracts the non-perturbing component from the
total force from the primary and companion, respectively. The magnitude of the
perturbing force is best estimated by the leading-order multipole term of the
force, which is the quadrupole force. Thus, we can approximate the magnitude
of Eq. (10) as444 A more detailed derivation of this estimate can be found in
Appendix B.2.

$\left|{\bf f}_{\rm pert}\right|\sim\frac{Gm_{p}m_{c}a_{c}^{2}}{Mr^{4}_{\rm
peri}}\ .$ (11)

Roughly speaking, one can consider two possible cases for the interaction
interval $\delta t$. The first, is the aforementioned periapsis interval,
under the approximation that the most significant perturbations (though still
weak) occur over some interval of time at which the test particle is near its
periapsis, $T_{\rm peri}$. Since the test particle’s susceptibility to energy
exchange is determined by its velocity (as per Eq. (4)), we estimate the
length of this interval by the time integral of the cross-track component of
the relative test particle velocity, $v_{\phi}$, over one period of its orbit,
divided by its maximum value, i.e.,

$T_{\rm peri}=\int_{0}^{P_{t}}\frac{v_{\phi}}{v_{\phi,\rm
max}}\,dt=(1-e_{t})P_{t}\ .$ (12)

If the companion does not move significantly during the interval of periapsis
passage, then the perturbing force can be treated as constant over $T_{\rm
peri}$. In this case, the energy change over the whole periapsis passage is
well-described by an impulse approximation, and the interaction interval
$\delta t$ is the same as $T_{\rm peri}$.

Label |  $m_{c}$ [$M_{\odot}$] |  $a_{c}$ [au] |  $r_{\rm peri}$ [au] | $e_{t}$ | No. Runs
---|---|---|---|---|---
DiffComp | $50-10^{6}$ | $10-10^{4}$ | $118.32$ | $0.884$ | $8308$
DiffTestP | $5000$ | $100$ | $16-4000$ | $0.2-0.997$ | $18440$
Table 2: The relevant numerical initial conditions. The results are shown in
Figs. 3 and 5. Note that in these Figures the square points represent an
average over $10$ realizations in terms of the initial orbital angles,
$\omega_{c,t},\Omega_{c,t},i_{c,t}$ and $f_{c,t}$.

On the other hand, if the companion moves during the test particle’s periapsis
passage, we can no longer use the impulse approximation. In particular, if the
companion’s period is sufficiently short, it may complete a significant
fraction of its orbit while the test particle is at periapsis, so that the
perturbing force ${\bf f}_{\rm pert}$ cannot be regarded as constant over the
interval $T_{\rm peri}$. In this case, the impulse approximation is, strictly
speaking, invalid. Nonetheless, an interaction timescale $\delta t$ can still
be estimated. In this case, we evaluate the change in energy per periapsis
passage by dividing the interval of periapsis passage into shorter intervals,
$dt$, for which the perturbing force is constant. Then, by integrating over
each of these intervals, the total energy change, $\delta E$, is

$\delta E={\bf v}\cdot\int_{t_{i}}^{t_{f}}{\bf f}_{\rm pert}\,dt=v_{\rm
peri}\int_{t_{i}}^{t_{f}}f_{\rm pert}^{\parallel}\,dt\ ,$ (13)

where $[t_{i},t_{f}]$ is the interval of periapsis passage, and $f_{\rm
pert}^{\parallel}$ is the component of the perturbing force parallel to the
test particle’s velocity during the interaction. Thanks to the companion’s
circular orbit, this component of the force varies roughly sinusoidally with
the companion’s mean anomaly, so that it may be approximated as

$f_{\rm pert}^{\parallel}\sim\left|{\bf f}_{\rm pert}\right|\cos n_{c}t\ ,$
(14)

where $n_{c}$ is the mean motion of the companion. Then Eq. (13) can be
approximated as

$\displaystyle\delta E$ $\displaystyle\,\sim\,$ $\displaystyle v_{\rm
peri}\int_{t_{i}}^{t_{f}}\left|{\bf f}_{\rm pert}\right|\cos n_{c}t\,dt$
$\displaystyle\,=\,$ $\displaystyle v_{\rm peri}n_{c}^{-1}\left|{\bf f}_{\rm
pert}\right|\sin n_{c}t\Big{|}_{t_{i}}^{t_{f}}\sim v_{\rm
peri}n_{c}^{-1}\left|{\bf f}_{\rm pert}\right|\ ,$

where the sine term represents the dependence on the orientation and thus
causes the spread in the distribution of possible $\delta E$. From Eq. (2.1),
it follows that the interaction timescale when the companion’s period is short
is

$\delta t\sim n_{c}^{-1}=\frac{P_{c}}{2\pi}\ .$ (16)

In general, the interaction timescale is given by the shorter of the two
cases,

$\delta t\sim\min\left(T_{\rm peri},\,\frac{P_{c}}{2\pi}\right)\ ,$ (17)

so that the former (latter) value is used if the companion’s movement is (is
not) negligible over the periapsis interval.

There are three regimes to consider when evaluating Eq. (6): (1) external
companion with a long period, (2) external companion with a short period, and
(3) internal companion with a short period. The three possible cases are
summarized in Table 1. Note that Kepler’s law $P_{c}\propto a_{c}^{3/2}$
forbids the case of an internal companion with a long period.

Evaluating Eq. (6) for the two external companion cases and the one internal
companion case then gives

$\displaystyle\frac{T_{\rm stab}}{P_{t}}$ $\displaystyle\sim$
$\displaystyle\frac{q^{-2}}{4(1+e_{t})}$ (18) $\displaystyle\times$
$\displaystyle\left\\{\begin{array}[]{lr}\alpha^{-4}(1-e_{t})^{3}/(2\pi)^{2}&\quad
T_{\rm peri}<P_{c}/2\pi,~{}\alpha<1\\\ \alpha^{-1}(1-e_{t})^{2}(1+q)&\quad
T_{\rm peri}>P_{c}/2\pi,~{}\alpha<1\\\
\alpha^{7}(1-e_{t})^{2}(1+q)^{4}&\quad\alpha\geq 1\end{array}\right.\ ,$ (22)

where $q$ is the mass ratio $m_{c}/m_{p}$ and $\alpha=a_{t}(1-e_{t})/a_{c}$
(see Eq. (1)).

###  2.2 Semimajor Axis Approach

The energy of the test particle orbit is proportional to the inverse of its
osculating semimajor axis; that is, $E\propto a_{t}^{-1}$. Thus a change in
$E$ by an order of itself necessarily means that $a_{t}$ has changed by an
order of itself as well. Similar to Eq. (4), a significant change in the test
particle’s semimajor axis occurs when

$\sigma(\Delta a_{t})=\sqrt{N}\sigma(\delta a_{t})\sim a_{t}\ .$ (23)

The change in semimajor axis at each periapsis passage, $\delta a_{t}$, can be
estimated by

$\delta a_{t}\approx\frac{da_{t}}{dt}\Big{|}_{\rm peri}\delta t\ ,$ (24)

where $da_{t}/dt$, given by the Lagrange Planetary Equation for semimajor axis
evolution (e.g., Murray & Dermott (2000b); Poisson & Will (2014); where we
adopt similar notation to Will (2021)), is:

$\displaystyle\frac{da_{t}}{dt}=$ $\displaystyle
2\sqrt{\frac{a_{t}^{3}}{GM(1-e_{t}^{2})}}\times$
$\displaystyle\left\\{e_{t}\sin f_{t}\,{\bf f}_{\rm
pert}\cdot{\bf\hat{i}}\,+(1+e_{t}\cos f_{t})\,{\bf f}_{\rm
pert}\cdot{\bf\hat{j}}\,\right\\}\ ,$

where $\bf{\hat{i}}$ and $\bf{\hat{j}}$ are the radial and cross-track unit
vectors of the test particle with respect to the central body.

Given the interaction timescale mentioned above (Eqs. (12) and (16)), $\delta
a_{t}$ is, at periapsis,

$\delta a_{t}\approx 2\sqrt{\frac{a_{t}^{3}(1+e_{t})}{GM(1-e_{t})}}\,{\bf
f}_{\rm pert}\cdot{\bf\hat{j}}\,\delta t\ .$ (26)

Information about the relative orientation of the two orbits is contained by
the ${\bf f}_{\rm pert}\cdot{\bf\hat{j}}$ term. The standard deviation of the
distributions of possible $\delta a_{t}$ is then on the order of

$\sigma(\delta a_{t})\sim
2\sqrt{\frac{a_{t}^{3}(1+e_{t})}{GM(1-e_{t})}}\left|{\bf f}_{\rm
pert}\right|\delta t\ ,$ (27)

which, when plugged into Eq. (23), for the various cases of $\left|{\bf
f}_{\rm pert}\right|$ and $\delta t$, yields the same timescale as Eq. (18).

Figure 2:  Sample time evolution of the test particle’s semimajor axis. The
semimajor axis evolution contains the same information as the energy
evolution. Left Panel; External companion. An external companion IMBH of 9000
$\textup{M}_{\odot}$ orbiting Sgr A* at 380 au, and a test particle star
orbiting Sgr A* with an initial semimajor axis of 1020 au and eccentricity of
0.884. The semimajor axis can increase and decrease with no clear bias, thus
justifying the approximation of the energy evolution as a random process. We
also show a “zoomed in” time evolution of the same system. Substantive changes
to the semimajor axis of the star only occur during the “peaks” in the true
anomaly evolution, corresponding to the periapsis passage of the star. Right
Panel; Internal companion.  An internal companion IMBH of 27000
$\textup{M}_{\odot}$ orbiting Sgr A* at 64 au, and a test particle star
orbiting the Sgr A*-IMBH binary with an initial semimajor axis of 1020 au and
eccentricity of 0.884. The system’s behavior is similar to the external case.

##  3 Comparison with Numerical Results

###  3.1 Description of the Initial Conditions and the Numerical Approach

As a proof of concept, we consider a system consisting of the galactic center
black hole Sgr A* as the primary mass $m_{p}$, a star in orbit around Sgr A*
as the test particle $m_{t}$, and a hypothetical IMBH companion to Sgr A* as
the companion $m_{c}$. Using HNBody (Rauch & Hamilton, 2002), we directly
integrate over 26,000 variations of this system assuming Newtonian
gravity555We do not include any post-Newtonian (PN) precession for the
purposes of this comparison. The 1st pN precession can stabilize the system
against secular perturbations (e.g., Naoz et al., 2013b) even in the cases
where the orbits are more compact than the hierarchical limit (e.g., Wei et
al., 2021; Faridani et al., 2021). However, the purpose of this numerical
analysis is to test our analytical timescale.. From the results, we may
calculate $T_{\rm stab}$ numerically by finding the first time at which the
test particle’s semimajor axis (energy) changes by a factor of 2 from its
initial value. Comparing this value to the analytic prediction then provides a
test of the analytic timescale. We note that not every destabilization, as
defined above, necessarily leads to an ejection of the particle from the
system. Ejection occurs when the energy changes from negative to positive
(i.e., $\Delta E>+E$), and is thus a subset of destabilization.

Figure 3:  Stability timescale as a function of companion parameters. Left:
$T_{\rm stab}$ as a function of the companion’s semimajor axis $a_{c}$ and its
mass $m_{c}$, as determined by direct N-body integration. The upper (lower)
panels describe an external (internal) companion. A colored dot represents a
change in S0-2’s energy by a factor of more than 2 within the maximum
integration time of $10^{7}$ yr, where the color determines $T_{\rm stab}$. A
colored square represents 10 systems with the same $m_{c}$ and $a_{c}$, such
that all 10 systems have a change in the test particle energy of more than a
factor of 2 within the integration time. The color of the square is the
geometric mean of $T_{\rm stab}$ over all 10 systems. A gray square represents
no such change in the test particle energy within the integration time.
Overlaid are lines of constant $T_{\rm stab}$, as predicted by Eq. (18). Some
systems outside the regime where $T_{\rm stab}<10^{7}$ yr still exhibit
significant change in their orbit, if they are within the region where the
test particle can enter the Hill sphere of the companion, i.e., below the
solid black line. The dashed black line marks where the Kozai quadrupole
timescale is ten times greater than the test particle period. Secular
approximations are reasonably applicable for systems above this line, causing
them to exhibit angular momentum oscillations which can suppress the random
walk nature of the energy evolution, leading to more stable systems than
predicted by the analytic timescale. Right: The normalized change in test
particle angular momentum at $t=T_{\rm stab}$. As predicted by the quadrupole
timescale, systems above this line exhibit significant changes in angular
momentum.

We run a grid of systems, varying the IMBH mass and separation. In all runs,
we use $m_{p}=4\times 10^{6}~{}\textup{M}_{\odot}$ and
$m_{t}=10~{}\textup{M}_{\odot}$, and set the companion’s eccentricity to
$e_{c}=10^{-3}\approx 0$.666 When $e_{c}=0$ exactly, the argument of periapsis
$\omega_{c}$ is ill-defined. Eccentricities of exactly zero are unlikely in
practice, so we avoid this pathological case by setting $e_{c}=10^{-3}$.  We
present two sets of numerical runs, as summarized in Table 2. In the runs
labeled DiffComp, we set the test particle’s initial orbital parameters to be
the same in all runs, namely those of the star S0-2, i.e., $a_{t}=1020~{}{\rm
au}\,,~{}e_{t}=0.884$, and $m_{t}=10$ M⊙, and vary the companion’s mass and
semimajor axis systematically. The test particle’s argument of periapsis
$\omega_{t}$, longitude of ascending node, $\Omega_{t}$, and inclination
$\iota_{t}$, are chosen such that they will have the observed values of S0-2
on the sky, and then projected into the invariable plane. The companion for
this set of runs is varied on a grid of $m_{c}\in[50,10^{6}]$ M⊙, and
$a_{c}\in[10,10^{4}]$ au.

In the second set of runs (labeled DiffTestP), we set the companion’s mass and
semimajor axis to be $m_{c}=5000$ M⊙ and $a_{c}=100$ au, respectively, and its
eccentricity to be $e_{c}=10^{-3}$. In this case, we arbitrarily choose the
inclination of the companion’s orbit to be close to zero777 When $\iota_{c}=0$
exactly, the longitude of ascending node $\Omega_{c}$ is ill-defined.
($0.001$ rad). This choice eliminates dependencies on the companion’s argument
of periapsis and longitude of ascending nodes. In this set of runs, the test
particle’s initial periapsis distance is then systematically varied from
$a_{t}(1-e_{t})\in[16,4000]$, and its eccentricity from $e_{t}\in[0.2,0.997]$.

When the orbital configurations are varied, in runs DiffComp or DiffTestP, we
chose the initial longitudes of ascending nodes, arguments of periapsis, and
mean anomalies from a uniform distribution between $0$ and $2\pi$, and the
mutual inclination from an isotropic distribution (i.e., uniform in cosine).

In Section 3.2 we discuss the numerical time evolution of the test particle
orbit, and how it justifies the assumptions used in deriving the analytic
timescale. In Sections 3.3 and 3.4 we separately analyze each set of numerical
runs and compare them to the analytic timescale predictions.

###  3.2 Time Evolution of the Test Particle Orbit

The numerical results justify our approximation of the test particle’s energy
evolution as a series of discrete, random “jumps” occurring each time the test
particle reaches its periapsis. Fig. 2 depicts two representative cases of the
test particle’s time evolution in the presence of an external and internal
companion (left and right panels, respectively). In particular, the inset
panels show a zoom-in on the part of the evolution that exhibits the energy
jumps described in Section 2. These jumps occur at the test particle’s
periapsis, as indicated by the peaks in its true anomaly evolution. Moreover,
Fig. 2 demonstrates that over time the test particle’s semimajor axis (energy)
can increase and decrease in a diffusive manner so that its evolution is well
described by a random process. We provide additional examples of the time
evolution of the test particle’s orbit in the Appendix (see Fig. 8).

###  3.3 Dependence of the Timescale on the Companion Parameters

In the DiffComp set of runs, we integrated a total of $8308$ configurations,
initially starting with the same test particle orbital configuration and
varying systematically the orbit of the companion. All integrations are
stopped at $10^{7}$ years.

As stated above, we record the time at which the test particle’s energy has
changed from its initial value by a factor of $2$. The results are shown in
the left hand side of Fig. 3, where the color coding shows the time at which
the test particle’s energy changed by a factor of $2$, ranging from $1$ to
$10^{7}$ years. Systems whose test particle energy never changes by a factor
of $2$ within $10^{7}$ yr are depicted as a gray square. We depict this value
as a function of the companion’s semimajor axis $a_{c}$ and mass $m_{c}$
(filled circles in Fig. 3). For each point on the grid, we initially run one
system with a value of $m_{c}$ and $a_{c}$. The top panel shows systems with
an external companion ($a_{c}>r_{\rm peri}$), and the bottom panel shows
systems with an internal companion ($a_{c}<r_{\rm peri}$).

On the right hand side of Fig. 3, we plot the fractional change in angular
momentum (i.e., $\Delta L/L$) at $t=T_{\rm stab}$. The symbols have the same
meaning as the left hand side panel.

Note that in the part of the parameter space where the system is very non-
hierarchical ($a_{c}\sim r_{\rm peri})$ and where $m_{c}\ll m_{p}$, the
dynamical behaviour is rather chaotic (e.g., Stone & Leigh, 2019; Ginat &
Perets, 2021a, b; Kol, 2021; Manwadkar et al., 2021). Therefore, in this part
of the parameter space we ran an additional 10 configurations with the same
values of $m_{c}$ and $a_{c}$, and numerically calculated $T_{\rm stab}$ for
these values by taking the geometric mean of the times at which the energy
changes by a factor of $2$. These grid points are marked as boxes in Fig. 3.
Over-plotted are the contours of the analytic timescales from Eq. (18).

For the majority of the parameter space, in particular where the system is
“least hierarchical” ($\alpha\sim 1$), we find good agreement between the
analytic and numeric results. This suggests both that our analytic result is a
good estimate of the stability timescale, and also that the instability is
sufficiently explained by energy jumps at periapsis. Close encounters of the
test particle with the other two bodies, while possible in principle, are
evidently not frequent enough to drive the instability of the system in this
regime.

There are two regions of the parameter space for which there is disagreement
between the analytic and numerical results: systems above the black dashed
line in Fig. 3, which are more stable than predicted, and some systems to the
left of the blue $T_{\rm stab}=10^{6}$ yr line, which are less stable than
expected.

In the former region, the semimajor axes of the two orbits are not well-
separated, but the secular approximation is nonetheless somewhat applicable.
Thus, the semimajor axis evolution of these systems cannot be treated as
random walks, as in Fig. 2, but oscillations that are roughly consistent with
secular behavior (i.e., quasi-secular).

Figure 4:  Time evolution of a sufficiently hierarchical system.  A system
consisting of an external IMBH with mass $9000$ $\textup{M}_{\odot}$ orbiting
Sgr A* at $a_{c}=980$ au, and a test particle star with $a_{t}=1020$ au,
$e_{t}=0.884$. Although the system is not strictly hierarchical, the system
exhibits behavior predicted by the secular approximation. The test particle’s
semimajor axis does not evolve like a random walk, but instead oscillates.
Large eccentric Kozai-Lidov oscillations in the test particle’s angular
momentum occur as well. Because the system does not behave like a random walk,
the analytic timescale fails.

As an example, we show the evolution of such a system in Fig. 4. In this
system, the energy (i.e., the semimajor axis) evolves not as a random walk,
but as an oscillation, never straying far from its initial value. 888 These
energy oscillations are similar to those noted by Antognini et al. (2014); Luo
et al. (2016); Grishin et al. (2016, 2018); Bhaskar et al. (2020), which take
place roughly on an orbital timescale and may eventually cause non-secular
evolution.  As expected, thanks to the Eccentric Kozai-Lidov mechanism (EKL
Kozai, 1962; Lidov, 1962; Naoz et al., 2013a), the system also exhibits large
oscillations in the test particle’s angular momentum, shown in the upper
region of the top right panel of Fig. 3. These angular momentum oscillations
take place roughly on the timescale of

$T_{\rm
quad}\sim\frac{16}{30\pi}\frac{m_{p}+m_{c}+m_{t}}{m_{c}}\frac{P_{c}^{2}}{P_{t}}(1-e_{c}^{2})^{3/2}\
,$ (28)

(e.g., Antognini, 2015). Systems for which this timescale is much longer than
the test particle’s period, i.e., by a factor of 10 ($T_{\rm quad}>10P_{t}$),
are typically governed by quasi-secular rather than random walk behavior.
Thus, their stability is better described by secular criteria (e.g., Ivanov et
al., 2005; Lithwick & Naoz, 2011; Katz & Dong, 2012; Naoz & Silk, 2014;
Antonini et al., 2014; Bode & Wegg, 2014; Antognini et al., 2014; Luo et al.,
2016; Grishin et al., 2018; Bhaskar et al., 2020). On the upper right panel of
Fig. 3, these systems are shown to have changes in angular momentum $\gtrsim$
$50\%$ and lie above the black dashed line.

For sufficiently large $\alpha$ (bottom panel, internal companion), the system
is also hierarchical, so that secular approximations can likewise be applied,
and the energy does not behave like a random walk. Here we expect that the
secular inverse Eccentric Kozai-Lidov mechanism (iEKL Naoz et al., 2017, 2020)
or the hierarchical stability timescale of Mushkin & Katz (2020) provides a
better stability criterion.

In the latter region, some systems undergo changes in their orbital energy on
shorter timescales than predicted by our analytical criterion. In such cases,
close encounters between the test particle and another body are important to
the system’s evolution. Such close encounters may only occur if the star’s
apoapsis is within a few Hill radii of the IMBH’s orbit. This condition can be
expressed as

$a_{c}<a_{t}(1+e_{t})+kR_{\rm Hill}\ ,$ (29)

where $R_{\rm Hill}$ is the Hill radius of the companion, defined as

$R_{\rm Hill}=a_{c}\left(\frac{m_{c}}{m_{p}}\right)^{1/3}\ ,$ (30)

and $k$ is a factor of order unity (similarly to planet-planet scattering,
e.g., Chatterjee et al., 2008a). In Fig. 3 this condition is shown as the
solid black curve (where we adopt $k=3$). Indeed, systems above this curve
very rarely exhibit a large change in energy. While such systems in which
close encounters are important do exist, they are relatively few and far
between in Fig. 3, whereas the general trend in Fig. 3 obeys the timescale as
predicted by energy jumps at periapsis. This provides further evidence that
energy jumps at periapsis, not close encounters, are the primary drivers of
instability.

Figure 5:  Stability timescale as a function of test particle parameters. Note
that unlike in Fig. 3, we vary the test-particle’s orbital parameters, not the
companion’s parameters, and thus the positions of the external and internal
cases are flipped. The meaning of the symbols is the same as in Fig. 3. The
solid black line shows the Hill sphere boundary; very few systems below this
line exhibit a significant energy change within the maximum integration time
of $10^{4}$ test particle periods.

###  3.4 Dependence of the Timescale on the Test-Particle Parameters

In this set of runs, DiffTestP, we integrated a total of $18440$
configurations, this time starting with the same companion parameters (see
Table 2) and systematically varying the parameters of the test-particle. We
again record the first time at which the test particle’s energy changes by a
factor of $2$, and depict this time as a function of the initial test particle
periapsis distance $r_{\rm peri}$ and its eccentricity $e_{t}$. We integrated
each system for $10^{4}P_{t}$, where $P_{t}$ is the initial test particle
period. Similarly, as for the DiffComp run, we run additional $10$
realizations for systems that have $\alpha\sim 1$. The results are shown in
Fig. 5. The color code here represents the time for changing the energy by an
order of itself, normalized by the test particle period, i.e., $T_{\rm
stab}/P_{t}$. In Fig. 5 the internal case is shown in the top panel and the
external case is shown in the bottom panel.

We again find good agreement between the analytic and numerical results. As in
DiffComp, we find that the analytic and numerical timescales agree for a
majority of the parameter space. Similarly to the DiffComp set of runs, the
agreement with the analytic timescale is best in the regime $\alpha\sim 1$,
while the analytic timescale begins to fail in the regimes $\alpha\gg 1$ or
$\alpha\ll 1$, i.e. when the system is hierarchical. Additionally, we see
systems to the right of the blue $T_{\rm stab}=10^{3}\,P_{t}$ line which
become unstable much faster than the analytic timescale predicts, which we
attribute to close encounters between the companion and the test particle.

##  4 Applications

###  4.1 The S-cluster in the presence of an IMBH

Recent gravitational-wave observations by the LIGO-Virgo collaboration have
now confirmed the existence of IMBHs (e.g., GW190521; The LIGO Scientific
Collaboration et al., 2020a, b). Specifically, our galactic center may harbor
IMBHs as the result of a possible minor merger with a low-mass or dwarf
galaxy, or even with a globular cluster. Such a scenario was considered by
Rashkov & Madau (2013), who suggested that if IMBHs serve as the seeds of
SMBHs in the center of galaxies, hierarchical galaxy evolution would yield
many IMBHs in our galaxy. Additionally, a combination of theoretical and
observational arguments suggest that an IMBH is expected to exist in the
central parsec of our galaxy (e.g., Hansen & Milosavljević, 2003; Maillard et
al., 2004; Gürkan & Rasio, 2005; Gualandris & Merritt, 2009; Chen & Liu, 2013;
Generozov & Madigan, 2020; Fragione et al., 2020; Zheng et al., 2020; Naoz et
al., 2020). Of particular interest are the proposed efficient formation of
IMBHs as a result of black-hole mergers (e.g., Fragione et al., 2021), or via
black-holes collisions with main sequence stars (e.g., Rose et al., 2022).

To constrain the parameter space of mass and semimajor axis of a hypothetical
IMBH , Naoz et al. (2020) used the long baseline of observations of the star
S0-2, located close to the SMBH Sgr A*. S0-2 has been observed for more than
two decades, and its orbit is sufficiently regular that, if there is a
companion to Sgr A*, it is either quite close to the main black hole, or well
outside the orbit of S0-2.

The stability analysis done here (see Fig. 3) suggests that the parameter
space for which the orbit of S0-2 will be stable over its lifetime (6 Myr, Lu
et al., 2013) is constrained to the left side of Fig. 3. Another way to
visualise the parameter space is depicted in Fig. 6. In the top panel of Fig.
6, we consider S0-2 parameters for a wide range of a companion semimajor axis
and mass ratio (with respect to the SMBH’s mass), depicted as the color-coded
lines. The y-axis represents the stability timescale normalized to the period
of S0-2. As shown in the Figure, the stable regime for a given mass ratio $q$
is wherever the corresponding colored line lies above the grey line denoting
the age of S0-2. Thus, for example, for $q=0.1$ ($m_{c}=4\times
10^{5}~{}M_{\odot}$), the allowed regions are $\alpha\lesssim 10^{-2}$
(companion IMBH well outside the periapsis of S0-2), and $\alpha\gtrsim 10$
(companion well inside the periapsis of S0-2). For $q=10^{-4}$
($m_{c}=400~{}M_{\odot}$) (very low-mass companion), the orbit of S0-2 is
stable for essentially any value of $\alpha$.

Figure 6:  Stability Parameter. Top: Number of stable orbits of S0-2 as a
function of $\alpha$, for companion-primary mass ratios between $10^{-4}$ and
$10^{-1}$, assuming an eccentricity of 0.884 for S0-2. The gray line shows the
age of S0-2. Regions where the gray line lies above a chosen colored line
(corresponding to a given companion-primary mass ratio) are excluded values
for $\alpha$, as in these cases S0-2’s orbit would be too unstable to survive
to its current age. Bottom: Number of stable orbits of a test particle with
the same semimajor axis as S0-2 as a function of $\alpha$ for test particle
eccentricities ranging from $0$ to $0.999$, assuming a companion-primary mass
ratio of 0.01 ($m_{c}=4\times 10^{4}~{}M_{\odot}$). The star in the bottom
panel represents the approximate position of Kepler-419 (see text for
discussion).

###  4.2 The stability of the S-star cluster

The S-star cluster is a collection of stars on nearly isotropic orbits, within
the innermost arcsecond of Sgr A* (e.g., Lu et al., 2009; Sabha et al., 2012;
Gillessen et al., 2012). These stars are eccentric and orbit-crossing, and are
estimated to be relatively young (e.g., Ghez et al., 2004, 2005b; Lu et al.,
2009, 2013; Gillessen et al., 2009, 2012; Yelda et al., 2014; Gillessen et
al., 2017; Do et al., 2019). Thus, in the context of this investigation the
question of the stability of the S-cluster translates to the timescale of the
stability. We simplify the question by considering perturbations to S0-2’s
orbit by other stars, one at a time.

Examining the top panel of Fig. 6 we see that the low mass ratio curve
($q=10^{-4}$) implies that even an object as massive as $400$ M⊙ will be able
to lurk in the presence of S0-2 without destabilizing its orbit during its
lifetime. Thus, other stars, naturally will not destabilize S0-2’s orbit as
well. We thus suggest that during the lifetime of the S-cluster, star-star
interactions do not change the stellar orbits’ energies by order of
themselves999Note that other processes such as collisions and tidal capture
also take place on longer timescales than the S-cluster lifetime (e.g., Rose
et al., 2020, 2022), thus the conclusion here can be extended to other
physical processes beyond two-body interactions..

###  4.3 Stellar and Planetary Systems

Since the stability timescale developed here depends only on products of
ratios and other unitless parameters, it is not limited to galactic nuclei and
may be applied in a wide range of systems. Consider, for example, a planetary
system. Interestingly, the forces between the S-stars and the SMBH are similar
to those in a planetary system with a Sun-like star in the presence of a
companion. In particular, the mass ratio between Sgr A* and the stars in the
S-cluster is similar to that of a star and its planetary system. Further, the
distance scales at the galactic center (SMBH and S-Star, and SMBH-companion)
are roughly the square of those for a star-planet, and a star-companion. Thus,
our stability timescale is also relevant in the case of planetary systems.

Significantly, observations suggest that most of the massive stars reside in
binaries or higher multiples (e.g., Raghavan et al., 2006; Raghavan &
Steprāns, 2010; Sana et al., 2012; Moe & Di Stefano, 2017; Moe & Kratter,
2021). In the planetary field one also often considers orbital crossing as a
sign of instability (e.g., Chatterjee et al., 2008b; Nagasawa & Ida, 2011;
Denham et al., 2019; Wei et al., 2021; Faridani et al., 2021). However, we
suggest that the stability question needs to be cast as the timescale to
instability in exoplanetary systems as well.

As an example, consider the bottom panel of Fig. 6. This Figure depicts a
system with mass ratio of $q=0.01$ between the primary and the companion, and
varies the test-particle eccentricity. Consider a system of a star and a brown
dwarf companion on a circular orbit. It is thus easy to see that a circular
planet (the test-particle), or even a planet with a non-negligible
eccentricity on a crossing orbit, will not have its energy changed by a factor
of itself for a timescale between $2000$ and $10^{9}$ orbital periods of the
planet, depending on the value of $\alpha$. This implies that some planetary
systems may even be observed near instability. Of course, a smaller companion
results in an even longer stability timescale.

A potentially relevant example may reside in the observed Kepler-419 system.
This system has two massive Jupiters orbiting a $1.4\,M_{\odot}$ star, with
$m_{t}=2.5\,M_{J}$, $e_{t}=0.83$, $a_{t}=0.37$ au, and $m_{c}=7.3\,M_{J}$,
$e_{c}=0.18$, $a_{c}=1.7$ au, respectively. (e.g. Dawson et al., 2012, 2014).
The system’s stability has been questioned recently by Denham et al. (2019)
and Jackson et al. (2019). For this system, $P_{t}=0.19$ yr, $P_{c}=1.84$ yr,
$T_{\rm peri}=0.11(P_{c}/2\pi)$, and our hierarchical parameter is
$\alpha=0.037$. Thus, we may ask on what timescale the system’s energy will
change by an order of itself. Using the first line of Eq. (18), we obtain
$T_{\rm stab}\sim 7\times 10^{4}$ yr. This is shown as the red star on the
bottom panel of Fig. 6. Therefore, at face value, this suggests that a major
change in the inner planet’s energy could occur within fewer than $10^{5}$
years.

We note that the inner planet cannot really be considered as a test particle
in this system. However, our numerical investigation was done for a non-
negligible mass inner star and is roughly consistent with the analytical
criterion (for comparable masses, our criterion underestimates the stability
timescale). Thus, our stability timescale for Kepler-419 may not be far from
reality.

##  5 Conclusion

The long-term stability of three-body systems is a fundamental problem in
astrophysics with many applications, from clusters of stars in galactic nuclei
to planetary systems. Furthermore, the instability of triple systems is
associated with energy exchange between the two orbits, and is a time-
sensitive concept. In this paper, we develop an analytic timescale at which
instability sets in for non-hierarchical triple systems. We find that
numerical N-body experiments corroborate this analytic result. In summary:

  1. 1.

The timescale for which the system is stable, $T_{\rm stab}$, is given by Eq.
(18), assuming a circular companion and highly eccentric test particle orbit
in a non-hierarchical configuration. This analytic result is validated
numerically in Figs. 3 and 5.

  2. 2.

For systems with such non-hierarchical configurations, the instability is
primarily driven by the high susceptibility to energy change the test particle
experiences near its periapsis, rather than by close encounters between the
companion and the test particle. The evolution of the test particle’s orbital
energy behaves like a random walk, as in Fig. 2.

  3. 3.

In the hierarchical limit, the random walk assumption on the energy evolution
does not hold, and so hierarchical systems are stable for much longer than
predicted by Eq. (18). In particular, their stability is predicted by secular
evolution. Systems do not need to be strictly hierarchical for the random walk
assumption to fail; they only need to be sufficiently well-described by the
secular approximation, such as in Fig. 4. Large oscillations in the angular
momentum caused by the Eccentric Kozai-Lidov effect are an indicator of
secular-like behavior.

EZ thanks Evgeni Grishin for helpful comments and discussion. SN acknowledges
partial support from NASA ATP 80NSSC20K0505 and thanks Howard and Astrid
Preston for their generous support. CMW is grateful for the hospitality of the
Institut d’Astrophysique de Paris, where parts of this work were carried out,
and acknowledges partial support from the National Science Foundation, Grant
Nos. PHY 19-09247 and 22-07681.

##  Appendix A The Stability Timescale at Lower Eccentricities

Figure 7:  Stability Timescale for Low Eccentricities. $T_{\rm stab}$ as a
function of $m_{c}$ and $a_{c}$ as determined by direct N-body integration,
for a test particle star of semimajor axis 1020 au, and with eccentricities
0.2 and 0.5. The analytic timescale overestimates the numerical timescale.

Our analytic calculation relies on the assumption that significant changes to
the test particle’s energy occur near its periapsis. For this assumption to
hold, it is necessary to assume a highly eccentric test particle orbit,
because of the dependence of $dE/dt$ on the test particle’s velocity, which is
sharply peaked about the periapsis only for eccentric orbits. When these
assumptions are violated, the analytic timescale does not strictly hold.
Nevertheless, we show numerically that the analytic timescale still has some
applicability at lower eccentricities.

To study the behavior of systems with low test particle eccentricity, we
integrate two sets of 1554 systems, with fixed test particle orbital
configuration and systematically varied companion orbital parameters. In one
set of systems, the test particle has an initial eccentricity of 0.2, and the
initial companion semimajor axis varies from 400 to 4000 au. In the other set,
the test particle has an initial eccentricity of 0.5, and the initial
companion semimajor axis varies from 200 to 4000 au. In both sets, the initial
test particle semimajor axis is 1020 au, the companion mass is varied from
$500$ to $10^{6}$ $\textup{M}_{\odot}$, and the initial orbital angles of both
the test particle and companion are randomly drawn uniformly in
$\Omega,\omega,f$ from $[0,2\pi]$ and uniformly in $\cos i$ from $[-1,1]$. All
systems are integrated for $10^{6}$ yr.

Our numerical results show that, for lower eccentricities, the analytic
timescale overestimates the numerically calculated timescale in regions of the
parameter space where orbital crossing is possible. Thus, the energy jumps at
periapsis, as described in the main body of the paper, cannot fully explain
the instability in the low eccentricity regime. For a test particle
eccentricity $e_{t}=0.5$, this discrepancy is not major. For $e_{t}$ as low as
$0.2$, the disagreement becomes more severe, with the analytic result
overestimating the numeric timescale by roughly a full order of magnitude.

##  Appendix B Supplemental Equations and Plots

We present supplemental equations and plots to our analytic derivation in
Section 2. This section will assume familiarity with the logical framework of
Section 2, but will provide additional rigor and justification for several
arguments made in the section. In particular, we will derive Eq. (4), the
energy jump per periapsis passage, and Eqs. (8) and (10), the perturbing force
vectors, directly from the Newtonian equations of motion, and provide
additional detail for our estimate of $\left|{\bf f}_{\rm pert}\right|$ for
internal companions in Eq. (11). We also provide an additional, more detailed,
description of the time evolution of the systems studied in the paper.

###  B.1  The Perturbing Force

The time derivative of $E$, as defined by Eq. (2), is, by direct
differentiation,

$\frac{dE}{dt}={\bf v}\cdot\left[\frac{d^{2}{\bf r}}{dt^{2}}-\frac{GM{\bf
r}}{r^{3}}\right]\ ,$ (B1)

where ${\bf r}$ and ${\bf v}$ are the position and velocity vectors of the
test particle relative to the central body, and $r$ is the magnitude of ${\bf
r}$.

The first term in brackets is the relative acceleration of the test particle,
and the second term is equivalent to the acceleration of the test particle due
to the central body. Thus the entire quantity in brackets, which is the
difference of the two, is the component of the test particle acceleration that
is not due to the central body. We call this term the perturbing force (per
test mass), i.e.,

${\bf f}_{\rm pert}=\frac{d^{2}{\bf r}}{dt^{2}}-\frac{GM{\bf r}}{r^{3}}\ .$
(B2)

When the perturbing force is defined this way, Eq. (B1) is equivalent to Eq.
(4).

Eqs. (8) and (10) then follow immediately from the Newtonian equations of
motion for ${\bf r}$. For an external companion, ${\bf r}={\bf r}_{tp}$, and
the corresponding equation of motion is

$\frac{d^{2}{\bf r}_{tp}}{dt^{2}}=-\frac{Gm_{c}{\bf
r}_{tc}}{r_{tc}^{3}}-\frac{Gm_{p}{\bf r}_{tp}}{r_{tp}^{3}}-\frac{Gm_{c}{\bf
r}_{cp}}{r_{cp}^{3}}\ .$ (B3)

For an internal companion, ${\bf r}={\bf r}_{t}$, and

$\frac{d^{2}{\bf r}_{t}}{dt^{2}}=-\frac{Gm_{p}{\bf
r}_{tp}}{r_{tp}^{3}}-\frac{Gm_{c}{\bf r}_{tc}}{r_{tc}^{3}}\ .$ (B4)

Figure 8:  Additional example time evolutions of the test particle’s orbit.
Each column shows the evolution of a system similar to those in the DiffComp
set of runs. For all systems $m_{c}=27000\,\textup{M}_{\odot}$, while $a_{c}$
varies between each system. Top Row: The evolution of the test particle’s
semimajor axis $a_{t}$ (red) and true anomaly $f_{t}$ (green). Middle Row: The
magnitude of the test particle’s relative velocity, $v$. Bottom Row: The ratio
between $f^{\parallel}$, the component of ${\bf f}_{\rm pert}$ parallel to the
relative test particle velocity $\bf v$, and $|{\bf f}_{\rm pert}|$, given by
Eqs. (9) and (11).

###  B.2  Estimating $\left|{\bf f}_{\rm pert}\right|$ for Internal Companions

When the companion is internal, the gravitational potential experienced by the
test particle can be expanded in a multipole series, by

$\Phi=-\frac{Gm_{p}}{\left|{\bf r}-{\bf
r}_{p}\right|}-\frac{Gm_{c}}{\left|{\bf r}-{\bf
r}_{c}\right|}=-\sum_{l=0}^{\infty}\frac{Ga_{c}^{l}\eta_{l}}{r^{\,l+1}}P_{l}(\cos\theta)\
,$ (B5)

where $P_{l}$ are the Legendre polynomials, $\theta$ is the angle between the
companion and the test particle’s position vectors, and

$\eta_{l}=\frac{m_{p}(-m_{c})^{l}+m_{c}m_{p}^{l}}{(m_{p}+m_{c})^{l}}\ .$ (B6)

The $l=0$ term is the potential due to the central body, so the perturbing
force is the negative spatial gradient of all remaining terms, which is

${\bf f}_{\rm
pert}=\sum_{l=2}^{\infty}{\bm{\nabla}}\left[\frac{Ga_{c}^{l}\eta_{l}}{r^{\,l+1}}P_{l}(\cos\theta)\right]=\sum_{l=2}^{\infty}\frac{Ga_{c}^{l}\eta_{l}}{r^{\,l+2}}{\bm{\xi}}_{l}\
,$ (B7)

where ${\bm{\xi}}_{l}$ is a vector, with magnitude of order unity, specifying
the direction of the force, given by

${\bm{\xi}}_{l}=P^{\prime}_{l}(\hat{\bf i}\cdot\hat{\bf I})\left[\hat{\bf
I}-(\hat{\bf i}\cdot\hat{\bf I})\hat{\bf i}\right]-(l+1)P_{l}(\hat{\bf
i}\cdot\hat{\bf I})\hat{\bf i}\ ,$ (B8)

where $\hat{\bf I}$ is the unit vector of the companion’s position, and
$\hat{\bf i}$ is the unit vector of the test particle’s position.

The leading term of the remaining series is the $l=2$ or quadrupole term,
since $\eta_{1}=0$. Then $\left|{\bf f}_{\rm pert}\right|$ is well estimated
by the leading order coefficient, $Ga_{c}^{2}\eta_{2}/r^{4}$, evaluated at
$r=r_{\rm peri}$, yielding Eq. (11). This estimate is only good to an order of
magnitude; the higher order terms, while decreasing, are not negligible since
the ratio in which the series is expanded in, $a_{c}/r$, is less than unity
but not small. Thus, the higher order terms contribute nontrivially to the
total force, but not at a larger scale than the leading order.

###  B.3  The Time Evolution

Fig. 8 shows additional examples of the time evolution of the test particle’s
orbit, for various regimes.

As was shown in Section 3.2, in all regimes, significant changes in the test
particle’s semimajor axis occur when its velocity peaks; that is, in a short
interval near its periapsis. However, as discussed in Section 3.3, such jumps
need not behave like a random walk. The two right columns in Fig. 8 show
systems in the quasi-secular regime. In the first system, changes in the
semimajor axis are dominated by periapsis jumps, but these jumps behave in an
oscillatory manner. In the second system, the jumps are small as the companion
is distant from the test particle’s periapsis, so that the semimajor axis
evolution is dominated by quasi-secular effects.

The bottom panels of each column indicate that the ratio between
$f^{\parallel}$, the component of ${\bf f}_{\rm pert}$ parallel to the
relative test particle velocity $\bf v$, and our estimate $\left|{\bf f}_{\rm
pert}\right|$ of the perturbing force, given by Eqs. (9) and (11), is of order
unity. This indicates that $\left|{\bf f}_{\rm pert}\right|$ is indeed a good
order-of-magnitude estimate of the perturbing force. Note that for external
companions, the estimate is valid at all points in the test particle’s orbit,
but for internal companions (represented by the leftmost panel in Fig. 8), the
estimate is only valid near periapsis, since it was derived with that
assumption.

## References

  * Antognini et al. (2014) Antognini, J. M., Shappee, B. J., Thompson, T. A., & Amaro-Seoane, P. 2014, MNRAS, 439, 1079, doi: 10.1093/mnras/stu039
  * Antognini (2015) Antognini, J. M. O. 2015, MNRAS, 452, 3610, doi: 10.1093/mnras/stv1552
  * Antonini et al. (2014) Antonini, F., Murray, N., & Mikkola, S. 2014, ApJ, 781, 45, doi: 10.1088/0004-637X/781/1/45
  * Bhaskar et al. (2020) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2020, The Astronomical Journal, 161, 48, doi: 10.3847/1538-3881/abcbfc
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bode & Wegg (2014) Bode, J. N., & Wegg, C. 2014, MNRAS, 438, 573, doi: 10.1093/mnras/stt2227
  * Chatterjee et al. (2008a) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008a, ApJ, 686, 580, doi: 10.1086/590227
  * Chatterjee et al. (2008b) —. 2008b, ApJ, 686, 580, doi: 10.1086/590227
  * Chen & Liu (2013) Chen, X., & Liu, F. K. 2013, ApJ, 762, 95, doi: 10.1088/0004-637X/762/2/95
  * Dawson et al. (2012) Dawson, R. I., Johnson, J. A., Morton, T. D., et al. 2012, ApJ, 761, 163, doi: 10.1088/0004-637X/761/2/163
  * Dawson et al. (2014) Dawson, R. I., Johnson, J. A., Fabrycky, D. C., et al. 2014, ApJ, 791, 89, doi: 10.1088/0004-637X/791/2/89
  * Denham et al. (2019) Denham, P., Naoz, S., Hoang, B.-M., Stephan, A. P., & Farr, W. M. 2019, MNRAS, 482, 4146, doi: 10.1093/mnras/sty2830
  * Do et al. (2019) Do, T., Hees, A., Ghez, A., et al. 2019, Science, 365, 664, doi: 10.1126/science.aav8137
  * Eggleton & Kiseleva (1995) Eggleton, P., & Kiseleva, L. 1995, ApJ, 455, 640, doi: 10.1086/176611
  * Faridani et al. (2021) Faridani, T., Naoz, S., Wei, L., & Farr, W. M. 2021, arXiv e-prints, arXiv:2107.07529.  https://arxiv.org/abs/2107.07529
  * Fragione et al. (2021) Fragione, G., Kocsis, B., Rasio, F. A., & Silk, J. 2021, arXiv e-prints, arXiv:2107.04639.  https://arxiv.org/abs/2107.04639
  * Fragione et al. (2020) Fragione, G., Loeb, A., Kremer, K., & Rasio, F. A. 2020, ApJ, 897, 46, doi: 10.3847/1538-4357/ab94b2
  * Generozov & Madigan (2020) Generozov, A., & Madigan, A.-M. 2020, ApJ, 896, 137, doi: 10.3847/1538-4357/ab94bc
  * Ghez et al. (2005a) Ghez, A. M., Salim, S., Hornstein, S. D., et al. 2005a, ApJ, 620, 744, doi: 10.1086/427175
  * Ghez et al. (2005b) —. 2005b, ApJ, 620, 744
  * Ghez et al. (2004) Ghez, A. M., Wright, S. A., Matthews, K., et al. 2004, ApJ, 601, L159, doi: 10.1086/382024
  * Gillessen et al. (2009) Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075, doi: 10.1088/0004-637X/692/2/1075
  * Gillessen et al. (2012) Gillessen, S., Genzel, R., Fritz, T. K., et al. 2012, Nature, 481, 51, doi: 10.1038/nature10652
  * Gillessen et al. (2017) Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30
  * Ginat & Perets (2021a) Ginat, Y. B., & Perets, H. B. 2021a, MNRAS, 508, 190, doi: 10.1093/mnras/stab2565
  * Ginat & Perets (2021b) —. 2021b, Physical Review X, 11, 031020, doi: 10.1103/PhysRevX.11.031020
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15, doi: 10.1051/0004-6361/201833718
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2020, A&A, 636, L5, doi: 10.1051/0004-6361/202037813
  * Grishin et al. (2018) Grishin, E., Perets, H. B., & Fragione, G. 2018, Monthly Notices of the Royal Astronomical Society, 481, 4907, doi: 10.1093/mnras/sty2477
  * Grishin et al. (2016) Grishin, E., Perets, H. B., Zenati, Y., & Michaely, E. 2016, Monthly Notices of the Royal Astronomical Society, 466, 276, doi: 10.1093/mnras/stw3096
  * Gualandris & Merritt (2009) Gualandris, A., & Merritt, D. 2009, ApJ, 705, 361, doi: 10.1088/0004-637X/705/1/361
  * Gürkan & Rasio (2005) Gürkan, M. A., & Rasio, F. A. 2005, ApJ, 628, 236, doi: 10.1086/430694
  * Hamilton & Burns (1991) Hamilton, D. P., & Burns, J. A. 1991, Icarus, 92, 118, doi: https://doi.org/10.1016/0019-1035(91)90039-V
  * Hansen & Milosavljević (2003) Hansen, B. M. S., & Milosavljević, M. 2003, ApJ-Lett, 593, L77, doi: 10.1086/378182
  * Hayashi et al. (2022) Hayashi, T., Trani, A. A., & Suto, Y. 2022, The Astrophysical Journal, 939, 81, doi: 10.3847/1538-4357/ac8f48
  * Hayashi et al. (2023) —. 2023, The Astrophysical Journal, 943, 58, doi: 10.3847/1538-4357/acac1e
  * Hills (1988) Hills, J. G. 1988, Nature, 331, 687, doi: 10.1038/331687a0
  * Ivanov et al. (2005) Ivanov, P. B., Polnarev, A. G., & Saha, P. 2005, MNRAS, 358, 1361, doi: 10.1111/j.1365-2966.2005.08843.x
  * Jackson et al. (2019) Jackson, J. M., Dawson, R. I., & Zalesky, J. 2019, AJ, 157, 166, doi: 10.3847/1538-3881/ab09eb
  * Katz & Dong (2012) Katz, B., & Dong, S. 2012, ArXiv e-prints.  https://arxiv.org/abs/1211.4584
  * Kol (2021) Kol, B. 2021, Celestial Mechanics and Dynamical Astronomy, 133, 17, doi: 10.1007/s10569-021-10015-x
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Lidov (1962) Lidov, M. L. 1962, planss, 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lithwick & Naoz (2011) Lithwick, Y., & Naoz, S. 2011, ApJ, 742, 94, doi: 10.1088/0004-637X/742/2/94
  * Lu et al. (2013) Lu, J. R., Do, T., Ghez, A. M., et al. 2013, ApJ, 764, 155, doi: 10.1088/0004-637X/764/2/155
  * Lu et al. (2009) Lu, J. R., Ghez, A. M., Hornstein, S. D., et al. 2009, ApJ, 690, 1463, doi: 10.1088/0004-637X/690/2/1463
  * Luo et al. (2016) Luo, L., Katz, B., & Dong, S. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3060, doi: 10.1093/mnras/stw475
  * Maillard et al. (2004) Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, A&A, 423, 155, doi: 10.1051/0004-6361:20034147
  * Manwadkar et al. (2021) Manwadkar, V., Kol, B., Trani, A. A., & Leigh, N. W. C. 2021, MNRAS, 506, 692, doi: 10.1093/mnras/stab1689
  * Mardling & Aarseth (2001) Mardling, R. A., & Aarseth, S. J. 2001, MNRAS, 321, 398, doi: 10.1046/j.1365-8711.2001.03974.x
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Moe & Kratter (2021) Moe, M., & Kratter, K. M. 2021, MNRAS, 507, 3593, doi: 10.1093/mnras/stab2328
  * Murray & Dermott (2000a) Murray, C. D., & Dermott, S. F. 2000a, Solar System Dynamics
  * Murray & Dermott (2000b) —. 2000b, Solar System Dynamics, ed. Murray, C. D. & Dermott, S. F.
  * Mushkin & Katz (2020) Mushkin, J., & Katz, B. 2020, Monthly Notices of the Royal Astronomical Society, 498, 665, doi: 10.1093/mnras/staa2492
  * Mylläri et al. (2018) Mylläri, A., Valtonen, M., Pasechnik, A., & Mikkola, S. 2018, MNRAS, 476, 830, doi: 10.1093/mnras/sty237
  * Nagasawa & Ida (2011) Nagasawa, M., & Ida, S. 2011, ApJ, 742, 72, doi: 10.1088/0004-637X/742/2/72
  * Naoz et al. (2013a) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013a, MNRAS, 431, 2155, doi: 10.1093/mnras/stt302
  * Naoz et al. (2013b) Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013b, ApJ, 773, 187, doi: 10.1088/0004-637X/773/2/187
  * Naoz et al. (2017) Naoz, S., Li, G., Zanardi, M., de Elía, G. C., & Di Sisto, R. P. 2017, AJ, 154, 18, doi: 10.3847/1538-3881/aa6fb0
  * Naoz & Silk (2014) Naoz, S., & Silk, J. 2014, ApJ, 795, 102, doi: 10.1088/0004-637X/795/2/102
  * Naoz et al. (2020) Naoz, S., Will, C. M., Ramirez-Ruiz, E., et al. 2020, ApJ, 888, L8, doi: 10.3847/2041-8213/ab5e3b
  * Petrovich (2015) Petrovich, C. 2015, ApJ, 808, 120, doi: 10.1088/0004-637X/808/2/120
  * Poisson & Will (2014) Poisson, E., & Will, C. M. 2014, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge: Cambridge University Press)
  * Raghavan et al. (2006) Raghavan, D., Henry, T. J., Mason, B. D., et al. 2006, ApJ, 646, 523, doi: 10.1086/504823
  * Raghavan & Steprāns (2010) Raghavan, D., & Steprāns, J. 2010, ArXiv e-prints.  https://arxiv.org/abs/1010.1226
  * Rashkov & Madau (2013) Rashkov, V., & Madau, P. 2013, ArXiv e-prints.  https://arxiv.org/abs/1303.3929
  * Rauch & Hamilton (2002) Rauch, K. P., & Hamilton, D. P. 2002, in AAS/Division of Dynamical Astronomy Meeting, Vol. 33, AAS/Division of Dynamical Astronomy Meeting #33, 08.02
  * Rose et al. (2020) Rose, S. C., Naoz, S., Gautam, A. K., et al. 2020, ApJ, 904, 113, doi: 10.3847/1538-4357/abc557
  * Rose et al. (2022) Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2022, ApJ, 929, L22, doi: 10.3847/2041-8213/ac6426
  * Roy & Haddow (2003) Roy, A., & Haddow, M. 2003, Celestial Mechanics and Dynamical Astronomy, 87, 411, doi: 10.1023/B:CELE.0000006767.34371.2f
  * Sabha et al. (2012) Sabha, N., Eckart, A., Merritt, D., et al. 2012, A&A, 545, A70, doi: 10.1051/0004-6361/201219203
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Stone & Leigh (2019) Stone, N. C., & Leigh, N. W. C. 2019, Nature, 576, 406, doi: 10.1038/s41586-019-1833-8
  * The LIGO Scientific Collaboration et al. (2020a) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020a, arXiv e-prints, arXiv:2009.01075.  https://arxiv.org/abs/2009.01075
  * The LIGO Scientific Collaboration et al. (2020b) —. 2020b, arXiv e-prints, arXiv:2009.01190.  https://arxiv.org/abs/2009.01190
  * Tory et al. (2022) Tory, M., Grishin, E., & Mandel, I. 2022, Publications of the Astronomical Society of Australia, 39, e062, doi: 10.1017/pasa.2022.57
  * Vynatheya et al. (2022) Vynatheya, P., Hamers, A. S., Mardling, R. A., & Bellinger, E. P. 2022, Monthly Notices of the Royal Astronomical Society, 516, 4146, doi: 10.1093/mnras/stac2540
  * Wei et al. (2021) Wei, L., Naoz, S., Faridani, T., & Farr, W. M. 2021, ApJ, 923, 118, doi: 10.3847/1538-4357/ac2c70
  * Will (2021) Will, C. M. 2021, Phys. Rev. D, 103, 063003, doi: 10.1103/PhysRevD.103.063003
  * Yelda et al. (2014) Yelda, S., Ghez, A. M., Lu, J. R., et al. 2014, ApJ, 783, 131, doi: 10.1088/0004-637X/783/2/131
  * Zheng et al. (2020) Zheng, X., Lin, D. N. C., & Mao, S. 2020, arXiv e-prints, arXiv:2011.04653.  https://arxiv.org/abs/2011.04653