11institutetext:  Anton Pannekoek Institute for Astronomy, University of
Amsterdam, 1090 GE Amsterdam, The Netherlands 22institutetext: Institute for
Gravitational Wave Astronomy, School of Physics and Astronomy, University of
Birmingham, Birmingham, B15 2TT, UK 33institutetext: Rudolf Peierls Centre for
Theoretical Physics, Clarendon Laboratory, Parks Road, Oxford, OX1 3PU, UK
44institutetext: CFisUC, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal 55institutetext: Instituto de Telecomunicações,
Campus Universitário de Santiago, 3810-193 Aveiro, Portugal 66institutetext:
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The
Netherlands

# Stellar triples on the edge

Comprehensive overview of the evolution of destabilised triples leading to
stellar and binary exotica

S. Toonen , , First and second author contributed equally 1122<EMAIL_ADDRESS>T. C. N. Boekholt ,⋆ , 334455<EMAIL_ADDRESS>S. Portegies Zwart 66

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context.  Hierarchical triple stars are ideal laboratories for studying the
interplay between orbital dynamics and stellar evolution. Both mass loss from
stellar winds and strong gravitational perturbations between the inner and
outer orbit cooperate to destabilise triple systems.

Aims.  Our current understanding of the evolution of unstable triple systems
is mainly built upon results from extensive binary-single scattering
experiments. However, destabilised hierarchical triples cover a different
region of phase space. Therefore, we aim to construct a comprehensive overview
of the evolutionary pathways of destabilised triple-star systems.

Methods.  Starting from generic initial conditions, we evolved an extensive
set of hierarchical triples using the code TRES, combining secular dynamics
and stellar evolution. We detected those triples that destabilise due to
stellar winds and/or gravitational perturbations. Their evolution was
continued with a direct N-body integrator coupled to stellar evolution.

Results. The majority of triples (54-69%) preserve their hierarchy throughout
their evolution, which is in contradiction with the commonly adopted picture
that unstable triples always experience a chaotic, democratic resonant
interaction. The duration of the unstable phase was found to be longer than
expected ($10^{3-4}$ crossing times, reaching up to millions), so that long-
term stellar evolution effects cannot be neglected. The most probable outcome
is dissolution of the triple into a single star and binary (42-45%). This
occurs through the commonly known democratic channel, during which the initial
hierarchy is lost and the lightest body usually escapes, but also through a
hierarchical channel, during which the tertiary is ejected in a slingshot,
independent of its mass. Collisions are common (13-24% of destabilised
triples), and they mostly involve the two original inner binary components
still on the main sequence (77-94%). This contradicts the idea that collisions
with a giant during democratic encounters dominate (only 5-12%). Together with
collisions in stable triples, we find that triple evolution is the dominant
mechanism for stellar collisions in the Milky Way. Lastly, our simulations
produce runaway and walk-away stars with speeds up to several tens of km/s,
with a maximum of a few 100km/s. We suggest that destabilised triples can
explain — or at least alleviate the tension behind — the origin of the
observed (massive) runaway stars.

Conclusions. A promising indicator for distinguishing triples that will follow
the democratic or hierarchical route, is the relative inclination between the
inner and outer orbits. Its influence can be summed up in two rules of thumb:
1) prograde triples tend to evolve towards hierarchical collisions and
ejections, and 2) retrograde triples tend to evolve towards democratic
encounters and a loss of initial hierarchy, unless the system is compact,
which experience collision preferentially. The trends found in this work
complement those found previously from binary-single scattering experiments,
and together they will help to generalise and improve our understanding on the
evolution of unstable triple systems of various origins.

###### Key Words.:

Methods: numerical – Stars: Binaries (including multiples) – Stars: kinematics
and dynamics – Stars: evolution – Stars: blue stragglers – Stars: white dwarfs
–

##  1 Introduction

While the Sun is a single star, other stars are often found as part of small
hierarchical stellar systems, consisting of two, three, or even more stars.
While stellar evolution has been studied for centuries, and binary evolution
has been studied for decades, our understanding of the evolution of triple
stars and higher-order multiples is still limited. It was Harrington (1968a,
1969) who studied the Von Zeipel-Lidov-Kozai mechanism for three stars,
describing the Hamiltonian up to third order, determining orbital stability
and demonstrating eccentricity-inclination variations in the inner orbit.
Later, the octuple term in the three-body Hamiltonian was derived and its
influence was investigated by Ford et al. (2000a) and Blaes et al. (2002a),
who also added relativistic precession in the inner orbit. It was demonstrated
by Naoz et al. (2011) that octuple order oscillations in the eccentricity and
inclination can lead to orbital flips, that is a transition from prograde to
retrograde inner orbits and vice versa. They also included tidal effects in
the equations of motion. Further progress was made through the many efforts of
the community (e.g. Kiseleva et al., 1998; Katz et al., 2011; Thompson, 2011;
Michaely & Perets, 2014; Naoz et al., 2016; Toonen et al., 2020; Hamers,
2020). It is an important step as hierarchical triples are abundant in the
Galactic field. For solar-type stars, one in four binaries has a third
companion (Raghavan et al., 2010; Tokovinin, 2014). And as these triples
evolve into interacting (mass-transferring) systems three times more often
than binaries (Toonen et al., 2020), triples should play an important role in
the formation of many close binaries. Moreover, for O- and B-type stars, the
fraction of triples is much higher (Evans, 2011; Sana et al., 2012) than for
solar-type stars. The triple fraction of massive stars even surpasses that of
massive-star binaries (Moe & Di Stefano, 2017).

The reason why triples evolve so different from binaries is that the presence
of the tertiary star can act as a catalyst for stellar interactions; three-
body dynamics (e.g. von Zeipel, 1910; Kozai, 1962; Lidov, 1962; Harrington,
1968a; Szebehely & Zare, 1977; Eggleton & Kiseleva, 1995; Naoz et al., 2013,
2016; Luo et al., 2016) lead to mass transfer in eccentric orbits (e.g.
Thompson, 2011; Shappee & Thompson, 2013; Michaely & Perets, 2014; Salas et
al., 2019; Stephan et al., 2019; Toonen et al., 2020), enhanced tidal
interactions (e.g. Mazeh & Shaham, 1979; Kiseleva et al., 1998; Fabrycky &
Tremaine, 2007; Liu et al., 2015; Bataille et al., 2018), and accelerates
stellar mergers that give rise to bright transients such as supernova type Ia,
gamma-ray bursts and gravitation wave emission (e.g. Wen, 2003; Perets &
Fabrycky, 2009; Thompson, 2011; Katz & Dong, 2012; Hamers et al., 2013;
Kimpson et al., 2016; Stephan et al., 2016; Portegies Zwart & van den Heuvel,
2016; Antonini et al., 2016, 2017a; Toonen et al., 2018). The dynamical impact
of the tertiary is strongest when the system is less hierarchical and close to
dynamically unstable. The hierarchy of an isolated stable triple can decrease
and even breakdown because of its own evolution (see Sect. 2 for an overview).
The breakdown and destabilisation of the triples is also known as the triple
evolution induced instability (TEDI, Perets & Kratter, 2012), and is the topic
of this paper.

Until recently, unstable triples were mostly considered in the context of
stellar clusters. A binary captures an interloper star temporarily, that is a
democratic resonance, after which the system dissolves again (Heggie, 1975).
In-depth studies on the energy exchange of the binary with passing by single
stars, and outcomes of democratic triple interactions, have been performed
extensively both numerically (e.g. Hills, 1975; Valtonen & Heggie, 1979; Hut &
Bahcall, 1983a; Leigh & Geller, 2012; Boekholt & Portegies Zwart, 2015;
Samsing & Ilan, 2018) and theoretically (e.g. Heggie, 1975; Monaghan, 1976;
Stone & Leigh, 2019; Ginat & Perets, 2021; Kol, 2021). They have provided us
with some statistical intuition for the outcome. For example, we know that the
least massive body is most likely to escape (but not always), thus leaving the
two more massive bodies in a binary system (e.g. Anosova, 1986; Mikkola, 1988;
Sterzik & Durisen, 1998). The eccentricity of the resulting binary is thermal
if the triple was initially close to being virial, or superthermal if it was
rather cold (Heggie, 1975). Its binding energy is independent of the mass
ratio and follows a power law with an index determined by the total angular
momentum (e.g. Boekholt & Portegies Zwart, 2015). The trajectories of the
stellar components evolve on dynamical timescales, such that stellar evolution
(SE) is typically of secondary importance.

In this paper we show that these expectations can not simply be extrapolated
to destabilised triples in the field. Due to the presence of chaos, most
individual binary-single scattering experiments tend to be unpredictable in
the sense that a small perturbation, such as a numerical error, is magnified
exponentially, resulting in diverged trajectories in phase space (Boekholt et
al., 2020), for example the triple loses memory of its initial condition.
Statistical distributions however seem to be predictable nevertheless and only
depend on the global conserved quantities, such as energy, angular momentum
and mass ratios (Valtonen & Karttunen, 2006). However, the destabilised
triples in the field have a very different type of initial condition compared
to unstable triples formed in clusters. First of all, the single third body is
on an elliptical orbit instead of a hyperbolic orbit with respect to the
binary. Secondly, the inner binary and the outer binary are usually separated
from each other such that they live on the stable side of the instability
border for hierarchical triples (Mardling, 2001). Generally, the distributions
of energy and angular momentum are thus expected to be different for
hierarchical triples compared to binary-single scatterings.

We follow-up on the work of Toonen et al. (2020) who simulated the evolution
of isolated hierarchical triples with realistic initial conditions of young
populations. Here we focus on those isolated triples that are born dynamically
stable, but become unstable as a result of the internal evolution of the stars
and by the dynamical interplay among the three stars. We perform self-
consistent simulations of the stable phase as well as the unstable phase
(Sect. 3). This includes full triple evolution calculations with TRES (Toonen
et al., 2016) for the former (similar to Toonen et al., 2020), and N-body
calculations linked with stellar evolution for the latter. We explore and
quantify the various outcomes in Sect. 4, such as stellar collisions and
runaway stars, place them in the observational context in Sect. 5, and
summarise our results in Sect. 7.

##  2 Evolutionary pathway

Except for periodic braids (Montgomery, 1998), dynamically stable triples are
hierarchical, in which two of the stars form a binary pair, and the tertiary
orbits the centre of mass of the binary with an orbital period much exceeding
that of the inner binary. The inner and the outer orbit can be approximated by
Keplerian orbits, and three-body dynamics manifests itself as an perturbation
on and interaction between these two orbits. At the lowest order of the
secular approximation (i.e. the quadrupole order) these are Lidov-Kozai cycles
(von Zeipel, 1910; Lidov, 1962; Kozai, 1962) in which the mutual inclination
and inner eccentricity vary periodically. At higher orders such as the
octupole level, more extreme inner eccentricities can be achieved, as well as
flips in the orbital orientation (see Naoz, 2016, for a comprehensive review).
The long-term evolution of hierarchical systems is governed by the interplay
between three-body dynamics and stellar evolution (Toonen et al., 2016, 2020).

The boundary between stable and unstable systems is not easily defined (or
definable), and therefore various stability criteria exist, which are based
a.o. on the concept of chaos, on the escape of one of the stellar components,
on zero velocity surfaces of the restricted three body problem, or on
numerical integrations (Georgakarakos, 2008, for a review). Here we have adopt
the stability criterion of Mardling & Aarseth (1999) including a factor that
reflects the dependence on inclination $i$ (Aarseth & Mardling, 2001):

$\displaystyle\begin{array}[]{l c l}\dfrac{a_{\rm out}}{a_{\rm in}}|_{\rm
crit}&=&\dfrac{2.8}{1-e_{\rm out}}(1-\dfrac{0.3i}{\pi})\\\ &&\\\
&&\left(\dfrac{(1.0+q_{\rm out})(1+e_{\rm out})}{\sqrt{1-e_{\rm
out}}}\right)^{2/5},\\\ \end{array}$ (4)

with $a_{\rm in}$ the orbital separation of the inner orbit, $a_{\rm out}$ the
orbital separation of the outer orbit, $e_{\rm out}$ the eccentricity of the
outer orbit, and the mass ratio $q_{\rm out}\equiv\dfrac{m_{3}}{m_{1}+m_{2}}$
with $m_{3}$ the mass of the star in the outer orbit (hereafter tertiary), and
$m_{1}$ and $m_{2}$ the masses of the two stars in the inner orbit (hereafter
primary and secondary respectively). By definition, the primary refers to the
initially most massive star, such that $m_{1}>m_{2}$ initially. Also $q_{\rm
in}\equiv\dfrac{m_{2}}{m_{1}}$. Triple systems are unstable if $\dfrac{a_{\rm
out}}{a_{\rm in}}<\dfrac{a_{\rm out}}{a_{\rm in}}|_{\rm crit}$. This stability
criterion is based on the concept of chaos and the consequence of overlapping
resonances, and later confirmed with direct n-body integrations to work well
for a wide range of parameters (Aarseth & Mardling, 2001; Aarseth, 2004; He &
Petrovich, 2018).

Eq. 4 shows that in general a hierarchical system can move over the stability
limit (or better cross into the instability region) when 1) $a_{\rm
out}/a_{\rm in}$ decreases, 2) the outer mass ratio $q_{\rm out}$ increases,
3) the outer eccentricity $e_{\rm out}$ increases, or 4) the mutual
inclination decreases. We mention five processes that can cause these orbital
parameters to change in a given triple:

  1. 1.

Three-body interactions: If we consider pure three-body dynamics, and focus on
the lowest order of the secular approximation, that is the quadrupole order,
the outer orbit does not change. This is generally not the case anymore for
higher order approximations. For example, when including the octupole term
(see review of Naoz, 2016) the outer eccentricity varies on a timescale111in
the limit of circular outer orbits of:

$\tau_{\rm oct}\sim\frac{256\sqrt{10}}{15\pi\sqrt{\epsilon_{\rm
oct}}}\,\tau_{LK},$ (5)

(Antognini, 2015) with the timescale of the Lidov-Kozai cycles (Kinoshita &
Nakai, 1999)

$t_{\rm LK}\sim\alpha\frac{P_{\rm out}^{2}}{P_{\rm
in}}\frac{m_{1}+m_{2}+m_{3}}{m_{3}}\left(1-e_{\rm out}^{2}\right)^{3/2},$ (6)

and the octupole parameter (Lithwick & Naoz, 2011; Katz et al., 2011;
Teyssandier et al., 2013; Li et al., 2014).

$\epsilon_{\rm oct}\equiv\frac{m_{1}-m_{2}}{m_{1}+m_{2}}\frac{a_{\rm
in}}{a_{\rm out}}\frac{e_{\rm out}}{1-e_{\rm out}^{2}}.$ (7)

Generally the octupole term is of importance if $|\epsilon_{\rm oct}|\gtrsim
0.001-0.01$.

  2. 2.

Stellar winds: Fast spherically symmetric winds that do not interact with the
stellar companions of the donor star have an adiabatic effect of the orbit.
From an angular momentum balance, one can show that the orbit widens as

$\frac{a^{\prime}}{a}=\frac{m_{1}+m_{2}}{m_{1}^{\prime}+m_{2}^{\prime}},$ (8)

where ’ denotes quantities after the wind mass loss. If mass is lost from the
inner binary of the triple, then

$\frac{a_{\rm out}^{\prime}}{a_{\rm in}^{\prime}}=\frac{a_{\rm out}}{a_{\rm
in}}\,\frac{m_{1}+m_{2}+m_{3}}{m_{1}^{\prime}+m_{2}^{\prime}+m_{3}^{\prime}}\,\frac{m_{1}^{\prime}+m_{2}^{\prime}}{m_{1}+m_{2}}<\frac{a_{\rm
out}}{a_{\rm in}}.$ (9)

In other words, as the fractional mass loss in the inner binary is larger than
that of the triple as a whole, the relative widening of the inner orbit is
larger than that of the outer orbit, such that the ratio $a_{\rm out}/a_{\rm
in}$ reduces in response to wind mass loss. When the two orbits approach one
another, the triple system can become unstable as described by Eq. 4. In this
way dynamically stable triple systems can become unstable due to their
internal stellar winds (e.g. Kiseleva et al., 1994; Iben & Tutukov, 1999;
Perets & Kratter, 2012).
We also note that the wind mass loss from the inner binary directly affects
the critical orbital separation ratio $\dfrac{a_{\rm out}}{a_{\rm in}}|_{\rm
crit}$ through its mass (i.e. $1+q_{\rm out}$) dependence (Eq.4). As the
stability criterion is only mildly dependent on mass, the direct effect of
mass loss on the stability of the triple is of secondary importance compared
to the orbital effect of stellar wind mass losses.

Another scenario induced by stellar wind mass loss is mass loss induced
eccentric Kozai cycles (MIEK) (Shappee & Thompson, 2013). In this case, the
triple exhibits no Von Zeipel-Lidov-Kozai cycles initially, but the mass loss
pushes the triple into a new region of parameter space where the cycles do
occur. During phases of extreme eccentricity, tidal forces and large radii
during the giant phase facilitate stellar collisions (e.g. Michaely & Perets,
2014; Naoz et al., 2016; Stephan et al., 2016, 2021). Therefore, depending on
the initial configuration of the triple, and the amount of mass loss, the
triple might follow the MIEK-track, or become dynamically unstable and even
disrupt.

  3. 3.

Mass transfer: Another internal process that can make a triple cross the
stability limit, is mass transfer in the inner binary (e.g. Kiseleva et al.,
1994; Iben & Tutukov, 1999; Freire et al., 2011; Portegies Zwart et al.,
2011). If we simplistically assume that the mass transfer is conservative, and
that stellar rotational angular momentum can be neglected compared to orbital
angular momentum, then angular momentum conservation dictates

$\frac{a^{\prime}}{a}=\left(\frac{m_{\rm d}\,m_{\rm a}}{m_{\rm
d}^{\prime}\,m_{\rm a}^{\prime}}\right)^{2},$ (10)

where $m_{\rm d}$ is the mass of the donor star and $m_{\rm a}$ that of the
accretor star. For a full description (including less stringent assumptions),
see Soberman et al. (1997) and Toonen et al. (2016). If the inner binary
widens in response to the conservative mass transfer, and approaches the non-
changing outer orbit, the triple system can cross the stability limit.

  4. 4.

Supernova explosions During the formation of a neutron star or black hole, a
kick may be imparted to the compact object, as suggested by Galactic NS
studies (Gunn & Ostriker, 1970; Cordes et al., 1993; Lyne & Lorimer, 1994;
Chatterjee et al., 2005; Hobbs et al., 2005; Becker et al., 2012; Verbunt et
al., 2017; Igoshev, 2020). Depending on the magnitude and orientation of the
kick it may unbind an orbit all together or simply change it. Besides this
direct affect on the inner and outer orbit of a triple, there is the
additional affect that a supernova explosion in the inner binary changes its
centre-of-mass (position, velocity, mass) which affects the outer orbit
(Pijloo et al., 2012). As a result the post-supernova trajectories of the
stars may result in unbound triples, unstable triples, or dynamically more
interesting triples (i.e. in which the importance of three-body dynamics has
increased) (Pijloo et al., 2012; Tauris & van den Heuvel, 2014; Lu & Naoz,
2019). An interesting prospect of a supernova explosion of the tertiary star
is the birth of a runaway binary (Gao et al., 2019).

  5. 5.

External effects External effects can be separated in two classes, the
perturbations induced by the background galactic tidal field with the
occasional relatively wide encounter with another star, and the strong
influence of nearby stars in a dense stellar cluster.

     * •

in the birth cluster The latter is most important in the earliest phase when
the triple system is still embedded in its birth clustered environment. This
phase was studied by van den Berk et al. (2007), who performed direct N-body
simulations of star clusters with primoridal triples, and binary evolution
(mass transfer in triples was not taken into account in these calculations).
Subsequent studies (Leigh & Geller, 2013; Fragione et al., 2020) studied the
effect of encounters in a clustered environment on wide outer triple orbits.

     * •

in the galactic field Due to their large cross sections, wide triples may be
disturbed by external effects even in the field. These include torques from
the Milky Way’s tide, rare, but strong encounters with passing field stars, or
a cumulative effect from many weak encounters which change the periastron of
the widest orbits through eccentricity pumping (Kaib & Raymond, 2014; Michaely
& Perets, 2016). Recently, especially the latter effect has been provoked to
lead to the formation of compact binaries and mergers (Michaely & Perets,
2016, 2019; Hamers & Thompson, 2019; Safarzadeh et al., 2020; Michaely &
Perets, 2020).

In this paper we consider triples that become unstable due to the first and
second process.

Table 1: Distributions of the initial stellar masses, rotation, and orbital
parameters. See also Fig. 1-3 in Toonen et al. (2020) for a graphical
representation of these parameters.

Parameter | Range | Model OBin | Model T14(1) | Model E09(2)
---|---|---|---|---
Mass of primary $m_{1}$ | 1-7.5$M_{\odot}$ | Kroupa IMF(3) | Kroupa IMF(3) | Kroupa IMF(3)
Mass ratio | $0-1^{*}$ | Uniform(4,5,6) | Uniform(4,5,6) | Model E09(2)
$m_{2}$, $m_{3}>0.008$$M_{\odot}$
Stellar rotation |  | Hurley rotation(7) | Hurley rotation(7) | Hurley rotation(7)
Orbital separation |  $5-5\cdot 10^{6}$$R_{\odot}$∗∗ | Uniform in log($a$)(8,9) | Log-normal∗∗∗(1,4,11) | Model E09(2)
Eccentricity | $0-1$ | Thermal distribution(10) | Thermal distribution(10) | Thermal distribution(10)
Inclination | $0-\pi$ | Uniform in cos($i$)(12) | Uniform in cos($i$)(12) | Uniform in cos($i$)(12)
Argument of pericentre | –$\pi-\pi$ | Uniform | Uniform | Uniform

222 ∗ Except for the outer mass ratio $q_{\rm out}$ in model E09; ∗∗ Formally
the maximum orbital separation for Model T14 is $10^{12}$$R_{\odot}$. Due to
the shape of the log-normal distribution, the contribution from systems with
$a_{\rm out}\gtrsim 5\cdot 10^{7}$$R_{\odot}$ is negligible; ∗∗∗ A log-normal
distribution with periods in days, mean $\mu=5$ and dispersion $\sigma=2.3$;
References: (1)Tokovinin (2014); (2)Eggleton (2009); (3) Kroupa et al. (1993);
(4)Raghavan et al. (2010); (5) Duchêne & Kraus (2013); (6) Sana et al. (2012);
(7) Hurley et al. (2000); (8)Abt (1983); (9) Kouwenhoven et al. (2007); (10)
Heggie (1975); (11)Duquennoy & Mayor (1991); (12)Tokovinin (2017);

##  3 Method

In this study we simulate the evolution of a population of stellar triples
with a two-step approach. In the first part of the triple’s evolution when the
system is dynamically stable, we use the triple population synthesis code TRES
(Toonen et al., 2016) which is based on the secular approach (von Zeipel,
1910; Lidov, 1962; Kozai, 1962; Naoz, 2016). In this phase-averaged approach,
the dynamics follow from a perturbative expansion of the Hamiltonian in terms
of $a_{rmin}/a_{\rm out}$. When the triple crosses the stability boundary
given by Eq. 4, and the secular approach is no longer valid, we switch to a
N-body based approach. Below we describe the set-up of each approach in
detail.

###  3.1 Hierarchical phase

In this work we employ the triple evolution code TRES to generate large
populations of triple systems on the zero-age main-sequence (MS), simulate
their subsequent evolution and extract those systems that become dynamically
unstable within 13.5 Gyr.

At every time-step, processes such as stellar evolution (including stellar
winds), tidal interaction, mass transfer, and three-body dynamics are
considered with the appropriate prescriptions. The (single) stellar evolution
is calculated with the SeBa-code (Portegies Zwart & Verbunt, 1996; Toonen et
al., 2012) which is based on the stellar evolutionary tracks of Hurley et al.
(2000). The tidal method used by TRES is based on the weak-friction
equilibrium tide model (Hut, 1981; Hurley et al., 2002) and is appropriate for
moderate eccentricities (Moe & Kratter, 2018) (see sect. 6 for a discussion).
Besides precession due to the secular dynamics, precession due to general
relativistic effects (Blaes et al., 2002b), by tides (Smeyers & Willems, 2001)
and by intrinsic stellar rotation (Fabrycky & Tremaine, 2007) are included in
TRES as well. Three-body dynamics is included based on the double-averaged
approach up to and including the octupole term (Harrington, 1968b; Ford et
al., 2000b; Naoz et al., 2013). TRES is incorporated into the Astrophysics
MUltipurpose Software Environment, or AMUSE (Portegies Zwart & McMillan, 2018,
see also amusecode.org). For a full description of TRES, see Toonen et al.
(2016).

To generate the initial population of triples, we apply three models (i.e.
OBin, T14 and E09) that are similar to those used in Toonen et al. (2020).
Model T14 and E09 are based on observed samples of triples constructed in
Tokovinin (2014) and Eggleton (2009) respectively, whereas model OBin is based
on our understanding of primordial binaries. We focus on primary stars in the
mass range $1\rm{M}_{\odot}\leq m_{1}<7.5\rm{M}_{\odot}$, such that the stars
experience significant stellar evolution in a Hubble time, while avoiding a
supernova explosion. An overview of the most important ingredients of the
models is given in Tbl. 1. For a detailed discussion on the models, see Toonen
et al. (2020). After drawing the initial stellar and orbital parameters from
the distributions given in Tbl. 1, we require that the triple is dynamically
stable initially, according to Eq. 4. In doing so, we neglect the dynamical
and stellar evolution on the pre-MS and its effect on the stability of the
triple system (but see e.g. Moe & Kratter, 2018). The adopted requirement
regarding dynamical stability biases the population towards larger ratios of
$a_{\rm out}/a_{\rm in}$ and away from high eccentricities (see also Fig. 1-3
in Toonen et al., 2020). The latter is in line with observations of wide
binaries (Tokovinin & Kiyaeva, 2016; Raghavan et al., 2010; Moe & Di Stefano,
2017).

To calculate Galactic event rates, we assume a constant star formation history
of 3$M_{\odot}$per year, a maximum stellar mass of 100$M_{\odot}$, and a
minimum primary mass of 0.08$M_{\odot}$. Furthermore, we assume that 15% of
stellar systems are triples, and 50% are binaries, leaving 35% of stellar
systems to be truly single stars (Tokovinin, 2014).

###  3.2  Dynamically unstable phase

Once a triple has become unstable according to the stability criterion, we
halt the integration with TRES, and continue with a direct N-body integration.
The last snapshot of TRES gives the age, masses, radii, wind mass loss rates,
and orbital elements, which serve as the initial condition for the direct
N-body code. Due to the secular nature of TRES, no mean anomaly is available.
In order to extend our sample and to investigate any dependencies on the mean
anomaly, we create 10 instances of each triple where the mean anomaly of the
outer orbit is uniformly sampled between 0 and 360 degrees, while the inner
orbit starts at apocentre.

We distinguish three types of direct simulations:

  1. 1.

pure N-body,

  2. 2.

N-body with uniform wind,

  3. 3.

N-body with stellar evolution (i.e. radius evolution and time-dependent mass
loss).

Type 1 is appropriate if the unstable phase is relatively short, such that
winds play a negligible role. It is often assumed that once a hierarchical
triple becomes unstable, the final disintegration follows after only a few
dynamical time scales. We show, however, that this assumption is somewhat
naive, as the duration of the unstable phase is on average thousands of
crossing times, even reaching up to tens of millions. The crossing time is the
time for three stars to encounter each other once during a democratic
resonance Here we define it as $t_{\rm cross}\equiv R_{v}/V$, with the
characteristic size $R_{v}=\frac{GM_{\rm tot}^{2}}{-2E_{\rm tot}}$ and
velocity $V=\sqrt{-2E_{\rm tot}/M_{\rm tot}}$, where the total mass $M_{\rm
tot}=m_{1}+m_{2}+m_{3}$ and total energy $E_{\rm tot}$ is the sum of the
kinetic and potential energy of the system.

In our type 2 ensemble, we therefore include a uniform wind mass loss rate. We
adopt the last value given by TRES, and assume it to be constant until the
mass has reached a value of $0.5\,\rm{M}_{\odot}$. Its radius, which was
constant up to that point, is then set to $0.01\,R_{\odot}$ and the object is
deemed a white dwarf. The main difference with the type I simulations is thus
the mass loss experienced by stars with winds, and the finite time frame
within which the giant stars retain their large radii and corresponding
collisional cross section.

Our type 3 ensemble can be considered as the most realistic one, as this
implements detailed mass-radius tracks obtained from the stellar evolution
code SeBa (Portegies Zwart & Verbunt, 1996; Toonen et al., 2012) as used in
TRES as well. This also allows stars which initially do not have any wind, to
become giants later on during the simulation. This was not possible for the
type 1 and 2 simulations. By comparing results from the three types of
simulations we can determine the role of stellar evolution during the unstable
phase of hierarchical triple stars.

The N-body code is the commonly used fourth-order Hermite integrator (Makino,
1991). Each time step, we update the masses and radii of the stars according
to the uniform wind mass loss rate, or the tracks from SeBa. For the type 3
ensemble, we first run the SeBa code separately, and store the mass-radius
tracks of the three stars in a table. This table is then fed to the N-body
code, which adds the wind at the time step level. Here the time step criterion
is adaptive and resolves both changes in the orbit (by taking the minimum time
scale proportional to relative distance divided by relative speed over all
pairs), and in the mass (by taking the minimum time scale proportional to mass
divided by wind mass loss rate). The direct integrations do not include
relativity, tides or other phenomena such as mass transfer or tidal
disruption. Although each of these ingredients affect the outcome of triple
evolution, it is our main aim here to study the interplay between dynamics and
stellar evolution. The results we present here also serve as a benchmark for
follow up studies on destabilised triples in which more physical ingredients
are included.

The stopping conditions for a simulation are the following:

  * •

Age has reached a Hubble time (Bound),

  * •

Unbound binary-single pair (Escape),

  * •

Bound binary-single pair separated by a parsec (Drift),

  * •

Overlap of stellar radii (Collision),

  * •

Three unbound stars (Ionisation),

  * •

CPU time exceeded maximum of 24 hours (CPU) (see app. D for motivation).

The criterion for escape is based on three constraints (Standish, 1971): the
single body is a certain distance away from the binary’s centre of mass
($\Delta r\geq 100a$, with $a$ the semimajor axis), the single body is moving
away from the binary’s centre of mass, and finally the single body has a
positive energy. If the third constraint was not met, but the separation
between the single body and the binary reached one parsec, then we also stop
the simulation and deem the triple to be a drifter. At these large separations
we can no longer treat the triple as being isolated. Drifters are expected to
occur when the single body is ejected to large separations, but would
turnaround at only very large distances (Orlov et al., 2010). Each time step
we check if any two stars overlap their radii, in which case we categorise the
triple as a collision. To analyse the collisional triples in further detail,
we store the collision components (primary, secondary, tertiary), their
stellar types at the moment of collision, and whether the non-collision
component would still be bound or not to the collision product, if replaced by
a conservative, centre of mass body. It is also theoretically possible that
the triple ionises into three single, unbound bodies. This outcome is very
rare and can only occur if the triple is already loosely bound, such that
stellar winds in the inner binary can become impulsive (Hills, 1983; Veras et
al., 2011; Toonen et al., 2017). The final stopping condition is not a
physically motivated one, but a practical one. From previous studies on the
lifetime of unstable triple systems (e.g. Orlov et al., 2010), it is known
that the lifetime distribution has an algebraic tail towards exceedingly long
lifetimes. Using direct integration, these systems will not fulfil any
physical stopping condition within a reasonable simulation time (see App. D).
In an iterative manner, we increase the maximum CPU time, and we show that if
it is set to 24 hours, we achieve a percentage of unfinished systems below 10%
in the worst case. We discuss this and other limitations of our direct model
in Sec. 6.2 and App. D.

During the simulation we keep track of whether the initial hierarchy of the
triple is preserved or not. If at any instance this is not the case, we tag
the triple as having lost its initial hierarchy, and deem it to be democratic
(Hut & Bahcall, 1983b). If the outcome of the simulation produced an escaper
or a drifter, we determine which initial components is now the single star. We
tag the triple as ‘Democratic-X’, with X denoting the single star (1= primary,
2=secondary, 3=tertiary). It is important to note however, that we found our
criterion to be rather strict as some triples are able to preserve their
initial hierarchy, even though the tertiary might temporarily have a close
encounter with either of the inner binary components during pericentre
passage. Therefore, some triples tagged as ‘Democratic-3’ could also have been
tagged as ‘Hierarchical’. The flag proved particularly useful as it shows that
the majority of triple systems tend to preserve their hierarchy, both in the
cases of a collision and a dynamical dissolution.

Figure 1: Overview of the outcome of triple evolution. About 2-4% of triples destabilise by themselves during their evolution (top panel). These triples commonly consist of three MS stars or contain one or two giant stars (middle right panel). The outcome of the destabilised phase is often an ejection of one of the stars, but collisions are common as well (middle left panel). Collisions preferably occur between MS stars (bottom panel). Fractions are based on our most advanced N-body simulations that includes stellar evolution. model OBin | model T14 | model E09 |
---|---|---|---
|  |  |
Figure 2:  Frequencies of the most abundant types of realistic triples that
become dynamically unstable during their evolution. The five main types are
shown in the outer ring, and subtypes in the inner ring as indicated by the
legend. The letters represent the stellar phases of the stars at the moment of
the dynamical instability; MS for a main-sequence star, G for a giant, WD for
a white dwarf. The first part of the name of the types represents the primary
star, the second part the secondary. The name of the subtypes have a third
part that represent the evolutionary state of the tertiary. The three pie
charts reflect the three models for the initial population of triples on the
zero-age main-sequence.

##  4 Results

###  4.1 Hierarchical evolution

model OBin
---
model E09
Figure 3: Triples that become dynamical unstable during their evolution are
born in a specific area of phase space close to the stability limit. The right
lower part of the figure is empty as this is the ’forbidden’ region of
dynamically unstable systems. The different types of triples are colour-coded
in the same way as in Fig. 2. Black small dots represent triples that did not
become dynamically unstable (see Toonen et al., 2020, for their evolution).
The x- and and y-axes show the initial distance the inner and outer orbit
would circularise to if tides were efficient.

From realistic simulations of triple evolution with TRES, we find that a few
percent of all triples become dynamically unstable due to their own internal
evolution in a Hubble time. The percentage is highest for model OBin
($4.2\%$), and lowest for T14 ($2.3\%$) and E09 ($2.5\%$). A full exploration
of the typical evolutionary pathways for stellar triples, and a comparison to
binary evolution, is presented in Toonen et al. (2020). Common triple pathways
include mass transfer between the stars of the inner binary, and mass transfer
from the tertiary onto the inner binary (Fig. 1).

The percentages mentioned above translate to a Galactic rate of 1-2 events
every 1000 yrs. We expect the event rate to increase for triples with higher
initial masses due to stronger stellar winds (especially early in the
evolution on the MS) and due to the higher triple fraction. If we follow Moe &
Di Stefano (2017) and adopt a triple, binary and single star fraction of
75$\%$, 20$\%$ and 5$\%$ (instead of $15\%$, $50\%$, and $35\%$) as valid for
O-type stars, the Galactic rate would already increase to 3-7 events per kyr.

The systems that we focus on in this study are born in dynamically stable
configurations, and cross the stability boundary simply as they evolve
(hereafter destabilised systems). There are three main causes for this; three-
body interactions, stellar winds, or a combination of the two. As stellar
winds are most potent for evolved stars333In this paper we focus on stars in
the low- and intermediate mass range, for which wind mass loss rates are small
on the MS. However, for high-mass star evolution, stellar winds play a crucial
role already on the MS., a significant fraction of triples crossing the
stability boundary contain evolved stars (Fig. 1). While the giant phases only
make up $\lesssim 10\%$ of a star’s life, 30-40% of destabilised triples have
a giant444We consider a star a giant if the star has evolved off the MS and
has not yet formed a compact object, i.e. ranging from the Hertzsprung gap to
the asymptotic giant branch. primary star at the onset of the instability.

Fig. 1 also shows that a large percentage of triples become dynamically
unstable with all stars are on the MS. This constitutes $37\%$ for model OBin,
and $60\%$ and $53\%$ for model T14 and E09 (Fig. 2). The difference in the
percentages is due to the characteristic orbital separations of each model
(Fig. 3, see also Fig. 2 in Toonen et al., 2020). As outer orbits tend to be
more compact in model T14 and E09, these models contain relatively more
destabilised triples with MS primaries (MS,MS) instead of giant primaries
(G,MS/G) or secondaries (WD,G). The latter cases occur in wide orbits ($a_{\rm
in}\gtrsim 10^{3}R_{\odot}$), such that Roche-lobe overflow before the
instability can be avoided. Hereafter when we refer to the stellar types of
the triple with ’(X,Y)’ or ’((X,Y),Z)’, X refers to the stellar type of the
primary star that is the initially most massive star of the inner binary, Y
refers to the secondary, and Z to the tertiary or outer star.

Lastly $2-6\%$ and $\lesssim 4\%$ of destabilised triples have a WD primary
star, and a MS or WD secondary star respectively (i.e. ’WD-MS’ and ’WD-WD’).
At the onset of the instability, the stellar winds are negligible in these
systems. The previous wind mass loss phase(s) has widened the inner orbit,
driving the triple closer to the stability limit. However, contrary to the
previously discussed evolution for ’G-MS/G’ and ’WD-G’ systems, the orbital
widening itself is not enough to trigger an instability. In stead, the new
configuration leads to stronger three-body effects and higher outer
eccentricities that make the system enter the instability region.

We note that for all types [(MS,MS), (G,MS/G), (WD,G), (WD,MS) and (WD,WD)]
three-body dynamics plays a role in the orbital evolution of the systems, and
therefore the diminished dynamical stability. The far majority of systems are
born in the octupole regime. The octupole parameter at birth $|\epsilon_{\rm
oct,init}|>0.001$ for $97-99\%$ and $|\epsilon_{\rm oct,init}|>0.01$ for
$76-81\%$ of all triples that become dynamically unstable. This means that it
is important to take three-body dynamics into account in the hierarchical
evolution of the triples, even for (G,MS) and (WD,G).

The delay times between the formation of the triple and the onset of the
instability range from a hundred Myr to several Gyr typically. Similarly to
Toonen et al. (2020), we also find some tertiary-driven interactions that
occur on shorter timescales. The issue is that if the delay time is a fraction
of the star-forming timescale, one may wonder if the interaction would not
have taken place before the zero-age main-sequence which we take as the start
of our simulations. If we assume the primary’s pre-MS timescale is 1$\%$
(0.1$\%$) of its MS timescale (Baraffe et al., 2002), than it concerns
$15-24\%$ ($10-17\%$) of the triples in the (MS,MS) channel. We do not find a
significant difference in the outcome of the dynamically unstable phase (Sect.
4.2) for these systems compared to the full population of triples with purely
MS components.

###  4.2 Evolution during the dynamically unstable phase

####  4.2.1 Duration

One may naively expect the dynamical phase of an unstable triple to be short,
that is of the order of several crossing times. However, the dynamical phase
of the triples considered here, which are on the boundary of dynamical
stability, typically lasts thousands or ten thousands of crossing times (Fig.
4). In other words, the dynamically unstable phase is not a short phase. As a
result other physical processes can play a role during the ’dynamical’ phase
as well, as we show in the following sections.

Previous studies of unstable triples mainly focused on encounters between a
binary and a single star in dense environments. Scatter experiments show that
unstable triples have a ‘decay’ time of order a hundred crossing times (e.g.
Mikkola & Tanikawa, 2007; Boekholt & Portegies Zwart, 2015), and an algebraic
tail towards very long lifetimes corresponding to long excursions of the
single star with respect to the binary (Orlov et al., 2010). The difference
with our destabilised triples can be understood due to their higher angular
momenta. Amongst others Boekholt & Portegies Zwart (2015) show that the
duration of the dynamical phase increases with initial angular momentum, and
also with mass ratio (more equal stellar masses). Moderate mass ratios are
expected for stable stellar triples (Tokovinin, 2014; Moe & Di Stefano, 2017).
On top of that, our hierarchical triples start off on the edge of instability,
contrary to randomly sampled triples, which can start off much deeper inside
the unstable regime.

Figure 4: Cumulative histogram of the duration of the dynamically unstable
phase in crossing times ($t_{\rm dyn}/t_{\rm cross}$). The figures show that
initially stable triples that have moved into the dynamically unstable regime
can remain there for long times. The figure represents model OBin where the
dynamics are modelled including stellar evolution. Other models show
qualitatively similar behaviour.

####  4.2.2 Outcomes

((MS, MS), MS) | ((G, MS), MS)
---|---
|
|
Figure 5: Outcomes of the dynamically unstable phase for the full population
of triples that become dynamically unstable during their evolution. Pie-charts
represent the N-body simulations with stellar evolution for model OBin. The
histograms show the outcomes for all models where different colours represent
the different dynamical methods (Sect. 3.2). Solid, dashed, and dotted line
styles represent the models for the initial population of triples OBin, T14
and E09, respectively. The results for each model are normalised to unity.

We find five possible outcomes of the dynamically unstable phase; 1) one of
the stars escapes the system, 2) one of the stars drifts away from the system,
3) there is a collision between two of the stars, 4+5) at the end of the
simulation the system is still bound and dynamically unstable. For 4) the
system has evolved for a Hubble time for 5) the maximum CPU time has been
reached (Sect. 3.2). Both escaping and drifter stars are ejected from the
system, but the former with high velocities such that the triple has become
unbound, and the latter with low velocities such that the triple is
technically still bound. The reason to consider drifters as an outcome is that
a separation was reached of one parsec at which point perturbations from the
Galactic environment cannot be neglected (Jiang & Tremaine, 2010). In the
following sections we discuss the characteristics of the various outcomes in
detail.

All outcomes are prevalent (Fig. 1). However their relative frequency varies
with stellar types and initial conditions (Fig. 5). The left column of Fig. 5
shows that unstable ((MS,MS),MS) triples frequently eject a stellar component
and dissolve into a separate binary and single star. Collisions between MS
stars are also common. Few unstable ((MS,MS),MS) survive for a Hubble time.
There is little difference between the simulations with the pure N-body method
and those with a constant wind mass loss rate for destabilised ((MS,MS),MS)
triples as the typical wind mass loss rate for MS stars are small in the
adopted mass range. The N-body simulations with stellar evolution show an
enhancement in the fraction of escapers due to the evolution of the stars off
the MS during the dynamical interaction.

The typical outcomes for unstable ((G,MS),MS) triples are very different from
that of ((MS,MS),MS) triples. Furthermore, there are large differences
depending on which method we use to model the dynamically unstable phase. This
demonstrates that it is important to take stellar evolution into account
during the dynamically unstable phase itself. The stellar winds from the giant
primary not only cause the system’s hierarchy to weaken during the
hierarchical phase, but continue to effect the dynamics of the system after
entering the dynamically unstable regime. The evolution during this phase is
therefore not purely dynamical. The effect can be seen in the relative large
fraction of ejections and small fraction of long-lived bound triples in the
N-body simulations that include stellar winds and evolution compared to that
of the pure N-body method (bottom right panel of Fig. 5).

Another striking difference in the outcomes of unstable ((G,MS),MS) triples is
the number of collisions (bottom right panel of Fig. 5). In the simulations
with the pure N-body method, where the stellar parameters are assumed to be
frozen in time throughout the dynamically unstable phase, collisions are
remarkable frequent compared to the other methods. The typical timescale of
the dynamically unstable phase is, however, longer than the lifetime of the
giant star. In conclusion, for these systems it is therefore vital to take the
change in stellar mass and radius into account.

####  4.2.3 Ejections: escape and drift

The destabilisation of a triple typically leads to the ejection of one of its
stars; either because a star escapes or drifts away. In our most sophisticated
dynamical models that include stellar evolution, an escape occurs for 42-45%
of destabilised triples and a drift for 12-22% (Fig. 1, see Tbl. 2 for other
models).

Figure 6: Cumulative histogram of the inclination at destabilisation for
triples that will eventually eject one of the stellar components. Triples with
prograde orbits tend to experience hierarchical encounters, whereas retrograde
orbits tend to lead to democratic encounters.

To investigate the dynamical interactions, we differentiate between triples
for which the hierarchy is maintained, and those that undergo democratic
resonances. We find a link between this and the orbital orientation at the
onset of the instability. Triples with prograde555We refer to triples as being
prograde if they have a relative inclination between the inner and outer orbit
smaller than 90∘, and as retrograde otherwise.  orbits typically experience
hierarchical encounters and eject their tertiary star (Fig. 6). Those with
retrograde orbits typically experience democratic encounters, during which
most often the secondary star is ejected (Democratic-2), see also Tbl. 3). We
expect that this is the case since retrograde orbits are more stable compared
to prograde orbits when all other parameters are kept the same (Eq. 4).
Therefore, in order for a retrograde triple to cross the stability limit, the
ratio of $a_{\rm out}/a_{\rm in}$ is therefore smaller compared to the
prograde case, which makes is easier for a democratic encounter to ensue.
Additionally, democratic encounters are more common when winds and stellar
evolution are taken into account in the dynamical modelling compared to the
pure N-body simulations. Stellar winds push the system deeper in the
instability region.

Figure 7: Which star is ejected from the destabilised triple? Is it the lowest
mass star of the system as is often assumed? The figure shows a cumulative
histogram of the mass of the ejected star over that of the lightest component
of the inner binary. If this ratio is smaller than one (left of black dashed
line), then the ejected star is the lightest component. Only for Democratic-2
encounters is it commonly the case that the lightest component is ejected.

Next we investigate which star is ejected from the triple. For their
velocities, and the formation of runaway stars, see Sect. 5.1. Naively, one
may expect that this would be the lowest mass star predominantly. Fig. 7 shows
this is not the case for destabilised triples. Hierarchical encounters do not
typically lead to the ejection of the lowest-mass star, neither do
Democratic-1 and Democratic-3 encounters. The naive picture only holds for the
Democratic-2 encounters; the lowest mass star is ejected in more than 85% of
the Democratic-2 encounters. The expectation of ejecting the lowest mass star
can be traced back to the seminal work of Hills (1975), who performed
numerical scattering experiments between binaries and field stars with inner
components of 3 solar mass each, and a tertiary of one solar mass. On the
other hand, for truly ergodic, equal mass triples, the probability to eject
any of the components is a third (Spitzer, 1987). The destabilised triples
considered here do not have extreme mass ratios $m_{\rm
ejected}/$min$(m_{1,bin},m_{2,bin})$; 50% of triples have mass ratios between
$\sim$0.6-1.7.

The ejected star is typically on the MS. At the onset of the instability this
comprises 79-87% of the systems. At the moment of dissolution the fraction has
slightly decreased to 72-82% as stars evolve of the MS. For the majority of
other systems, the ejected star was initially on a giant branch during the
destabilisation (i.e. for 12-16% of systems). However, at the moment of
ejection it is very rare that the ejected star is still a giant star. At the
moment of ejection, 17-28% of ejected stars has already reached the WD stage.
The formation rate of single WDs through this channel in the Milky Way is
$(1-4)\times 10^{-4}$ per year, where as the birth rate of WDs from single
stellar evolution is approximately a few 0.1 per year (Verbeek et al., 2013;
Toonen et al., 2017). As the ejected stars are typically the original
secondaries and tertiaries of the triples, the ejected WDs do not follow the
same mass distribution as WDs from single stellar evolution.

Next, we focus on the remaining binary and how the orbit changes due to the
dynamical interaction. In the simplest scenario of purely dynamical ejections
without stellar exchanges (i.e. hierarchical and Democratic-3 encounters), we
expect the orbital separation to shrink due to energy conservation. There are
two exception two this scenario. Firstly, if stellar winds play an important
role, than the orbital separation also experience widening. This is clearly
visible as the difference between the top and bottom panel of Fig. 8.
Secondly, if stellar exchanges take place than the orbit can widen as well.
This can be seen in the top panel of Fig. 8 in the red and green markers that
lie above the black line. This can be understood in the following way:
consider a triple that dissolves into a single star and a binary with a
specific post-dissolution binding energy $E_{bin}$ and $m_{2}<m_{3}$. In a
hierarchical interaction $E_{bin}\equiv\frac{m_{1}m_{2}}{2a_{bin,hier}}$,
where as for a Democratic-2 encounter
$E_{bin}\equiv\frac{m_{1}m_{3}}{2a_{bin,dem-2}}$. From $m_{2}<m_{3}$, it
follows that if the low-mass secondary is ejected the orbital separation
$a_{bin,dem-2}>a_{bin,hier}$. Additionally, the post-interaction
eccentricities of the remaining binary are shown in Fig. 9. Where as
hierarchical interactions, as well as Democratic-3 encounters, lead to
eccentricity distributions that are more or less thermal (black line),
Democratic-2 encounters give rise to a distinct eccentricity distribution that
is significantly sub-thermal, that is larger fraction at smaller
eccentricities. We observe that when the lowest-mass body escapes, that the
remaining binary is not only more massive, but also statistically wider.
Combined with the smaller eccentricities, these trends all tend to maximise
the angular momentum of the binary. This implies that when the lightest body
is ejected after a democratic resonance, it does so on a low angular momentum
orbit.

Finally, destabilised triples with compact inner orbits ($a_{\rm in}\lesssim
10^{4}R_{\odot}$), tend to have prograde orbits, and therefore keep their
hierarchy and eject the tertiary star (Fig. 8). Given that the orbits are
compact, the ejected star escapes from the triple; very few triples experience
a drift, nor do they remain bound for a Hubble time. The retrogradely orbiting
triples do not typically lead to democratic interactions and ejections, but to
collisions instead (see also App.B).

N-body simulations
---
N-body + stellar evolution
Figure 8: Change in orbital separation for destabilised triples that unfold
into a binary and single star. The x-axis shows the orbital separation of the
inner orbit at the onset of the instability, and the y-axis shows the post-
interaction orbital separation of the binary. The results of the pure N-body
simulations are shown in the top panel, those of simulations that include
stellar evolution are shown on the bottom. The black dash-dotted line
represents no change in the orbital separation. The figure shows that besides
adiabatic wind mass losses leading to orbital widening, ejecting a star from
the inner orbit also tends to lead to orbital widening (see text). The two
panels represent model OBin. Model T14 and E09 show qualitatively similar
trends.  Figure 9: Cumulative histogram of the eccentricity of the remaining
binary after the ejection. Democratic-2 encounters lead to lower
eccentricities then hierarchical, Democratic-1, and Democratic-3 encounters.
The black dashed and dash-dotted line represents a thermal and superthermal
distribution, respectively.

In a small number of simulations, there was no remaining binary. Both orbits
are dissolved during the interaction. This is not possible in purely dynamical
situations due to energy conservation. However, it is known that impulsive
winds can unbind a binary (Veras et al., 2011; Toonen et al., 2017; Boekholt &
Toonen, In preparation), or in this case a triple. Not surprisingly, the
particular triples have wide orbits ($a_{\rm in}\sim 10^{5}-10^{6}R_{\odot}$),
and the stars are close to their apocentre. As a result their orbital
velocities are relatively low (and the impact of the winds relatively large).
In the dynamical simulations with wind or stellar evolution, it comprises
0.04% - 0.027% of all destabilised triples. The rate is the highest in model
OBin, as in this model we start the simulations with relatively many wide
triples.

####  4.2.4 Collisions

One of the exciting outcomes of a triple instability, is the collision of two
stars. Typically the collision involves the two stars that initially were in
the inner orbit, i.e the primary and secondary star. During the dynamically
unstable phase leading up to the collision, the original hierarchy of the
system is routinely preserved. 87-95% of destabilised colliding systems never
experience a democratic resonance. And even for those triples in our ensemble
which are flagged as having lost their hierarchy, stellar exchanges are rare
(i.e. 1.6-2.4% of all collisions). The inner orbits exhibits extreme behaviour
with median eccentricities at the time of collision of more than 0.996 in all
models. This is consistent with results from Bhaskar et al. (2021a), who find
that in the regime of similar mass stars ($m_{2}/m\geq 0.1$) and a semi-major
axis ratio $a_{\rm{in}}/a_{\rm{out}}\geq 0.3$, the maximum eccentricity
exceeds those estimated from secular theory, reaching up to 0.99 and beyond.

Even though the hierarchy is usually maintained in the dynamically unstable
phase, we see fluctuations in the configuration of the outer orbit. Not only
the eccentricity of the outer orbit varies (also in cases without stellar
winds), as expected in the eccentric Lidov-Kozai regime, but also the orbital
separation varies. This illustrates that the double-averaging approximation,
as used in TRES during the hierarchical phase, breaks down in the dynamically
unstable phase, justifying the transition to N-body simulations. Accordingly,
there are variations in the mutual inclination of the orbits, commonly
including orbital flips. It is the breakdown of the orbits that drives the
collisions in an extraordinary efficient way.

The fraction of destabilised triples that experience a collision ranges
between 13 and 35% (Tbl. 2). The pure N-body models give the highest fractions
(32-35%). The models that take stellar evolution into account during the
dynamically unstable phase provide a much lower fraction, that is 13-23%
depending on the initial distributions of triples (Fig. 5). This can be easily
understood as the radius of the evolving star decreases as it becomes a WD,
which decreases the cross section for the collision. In fact, the fraction of
triples with collisions in ((MS,MS),MS) systems does not change with the
various dynamical modelling methods, whereas the fraction for ((G,MS),MS)
systems decreases with a factor 3.5-18 when stellar evolution is taken into
account (Fig. 5). The systems that avoid the collision typically have wide
orbits ($a_{\rm in}\gtrsim 5\times 10^{4}R_{\odot}$, Fig. 10) and dissolve
into a binary and a single star through a stellar escape.

Figure 10: The origin of destabilised triples and their outcomes. Where as
drifters and long-lived triples originate from triples with large initial
orbital separations, colliding stars and CPU-limited systems come from compact
initial triples.

With our Galactic model as described in Sect. 3.1, we estimate a collision
rate of 2-3 per $10^{4}$ yr based on the N-body simulations with stellar
evolution. The majority of these involve collisions between MS stars (Fig. 1);
76-94% depending on the initial conditions. Collisions between a giant and a
MS star occur in 2-6%. More common are merger between WD and MS stars (4-15%)
that mostly come from ((G,MS),MS) systems whose giant primary has evolved into
a WD by the time of the collision. We also find low rates of giant-WD and WD-
WD collisions ($\sim$1% and $\sim$0.5-2% respectively). Assuming mass is
completely conserved during the merger, the merger remnant usually remains
bound to the tertiary star. Less than 0.5% of collisions in all models lead to
an unbound pair. For a discussion on the remaining binary, see Sect. 5.

####  4.2.5 Long-lived bound unstable triples

At the end of our simulations, a small fraction of our destabilised triples is
still bound, that is 7-23%. The vast majority of these systems has preserved
their initial hierarchy; only in 0.2-0.5% of cases was the system flagged as
undergoing a democratic encounter. However, this did not lead to the
disruption of the triple.

Similar to the previous section on collisions, we find that the double
averaged method breaks down for these systems. Besides strong variations in
the inner and outer eccentricity, we also notice variations in the outer semi-
major axis (also in cases without stellar winds). As a result of these
variations, the triple can move back and forwards over the stability limit
spending some time just inside the instability region, and some time just
outside (top panel of Fig. 11). In the pure dynamical simulations, 90% of
systems are within a factor 1.13-1.19 of $(a_{\rm out}/a_{\rm in})|_{\rm
crit}$ (Eq. 4). In the dynamical simulations with stellar winds and evolution,
this increases to 1.7-2.1 and 2.4-2.8 for simulations including stellar winds
and evolution, respectively. When including stellar winds and evolution, the
additional wind mass loss can push the system deeper in the instability region
(when mass is lost in the inner ’orbit’) or out of the instability region
(when the ’tertiary’ star loses mass). The latter is reminiscent of the
secular evolution freeze out studied by Michaely & Perets (2014).

The necessity of taking stellar evolution into account in the modelling of
these systems is also apparent in the bottom panels of Fig. 5. The fraction of
bound triples after a Hubble time decreases correspondingly as the stellar
winds push the system further into the instability region (Sect. 2).

The long-lived triples have wide orbits initially, that is initial inner
orbital separations $\gtrsim 10^{4}R_{\odot}$ (Fig. 10). Even though a wide
orbit implies a relatively long crossing time, the median of the evolved time
was still $10^{5}$ crossing times (Fig. 4). This is significantly more than in
the case of ejections or collisions, which indicates that the long-lived
triples have remained near the instability border with only mild variations in
its orbital elements.

At the onset of the dynamically unstable phase, the long-lived triples
typically have a stellar component on the giant branch that is (G,MS/G) and
(WD,G). For example, this is the case for 76%, 81%, and 64% of the triples in
model OBin, T14, and E09 respectively for dynamical simulations including
stellar evolution. At the end of the simulation, when a Hubble time has past,
these giant components have become white dwarfs. If the other stellar
components had masses above 1$\rm{M}_{\odot}$ initially, they too move off the
MS and become white dwarfs by the end of the simulation. This is the case for
60-80% of the ((G,MS,MS)) that remain bound in a Hubble time. Furthermore, all
destabilised triples with (WD,WD) inner binaries remain bound within a Hubble
time.

---
Figure 11: For those destabilised triples that remain bound, where are those
systems compared to the stability limit? This is quantified by the ratio of
$a_{\rm out}/a_{\rm in}|_{\rm crit}$ (that is the stability criterion of Eq. 4
over the current configuration (that is $a_{\rm out}/a_{\rm in}$). When the
ratio is larger than 1, the system is considered to be in the dynamically
unstable regime. If the ratio is smaller than 1, the system is dynamically
stable. On the top we show the systems that remain bound after a Hubble time,
on the bottom those systems that remain bound after simulating the dynamically
unstable phase for 24 hours. Solid, dashed, and dotted line styles represent
the models for the initial population of triples OBin, T14 and E09,
respectively.

####  4.2.6 CPU-time limited triples

Lastly, we discuss the evolution of those destabilised triples for which the
simulation of the dynamically unstable phase was terminated after simulating
for 24 hours. The fraction of systems that is CPU-limited depends on the
initial population of triples, and to a lesser degree on the modelling of the
dynamically unstable phase. For model OBin it comprises 2.1-3.9% of
destabilised triples, 8.8-10.1% for model T14, and 4.9-6.1% for model E09. We
have limited the simulation of each system in such a way that the total
fraction less than 10% for all models. In App. D, we demonstrate that in order
to significantly reduce the fraction of CPU-time limited triples it would
require over $10^{3}$ CPU-hours, which is beyond the scope of this project.
Our rate estimates of the most common outcomes of destabilised triple
evolution are therefore accurate to within an order of magnitude, while that
of less frequent channels should be considered a lower limit.

Many crossing times are simulated for the CPU-time limited triples, that is
typically over $10^{6}$ (Fig. 4), however the systems do not change in a
catastrophic way. The systems predominantly keep their hierarchy throughout
the simulation (Tbl. 2) and stay close to the stability limit (bottom panel of
Fig. 11). 75% of systems are within a factor 1.10-1.24 of the stability
criterion (Eq. 4) for all models except one. The factor increases to 1.88 for
the model that combines N-body dynamics with stellar evolution for the initial
triple population of T14. In this model we find more systems that have moved
into the dynamically stable regime (bottom panel of Fig. 11). In addition,
similarly to colliding systems and long-lived systems, the double averaged
method breaks down for the cpu-limited systems, which demonstrates the
necessity and justifies the transition from the secular simulations of TRES to
the N-body calculations.

CPU-limited triples have initial inner orbital separations between
$10^{2}-10^{4}R_{\odot}$ and outer orbital separations
$10^{3}-10^{5}R_{\odot}$ (Fig. 10). As more such triples are created in model
E09 compared to model OBin, the fraction of CPU-limited systems is higher in
the former. In addition, the initial eccentricities of the outer orbits are
remarkably high with a median value of 0.89-0.94 across all models. Given
these relatively tight inner orbits, it is not surprising that the majority of
CPU-limited triples are ((MS,MS),MS) systems (Tbl. 2). The small fraction of
((G,MS),MS) triples decreases further when stellar winds and evolution are
taken into account in the dynamical simulations.

##  5 Observational exotica

###  5.1 Runaway stars

The destabilisation of a triple often leads to the ejection of one of the
stars, and the formation of runaway and walkaway stars. These are stars that
are moving through space with abnormally high velocities compared to their
surroundings. Runaway stars are stars with velocities above 30(-40) km/s
(Blaauw, 1961; De Donder et al., 1997; Portegies Zwart, 2000; Dray et al.,
2005; Eldridge et al., 2011; Boubert & Evans, 2018; Renzo et al., 2019;
Bischoff et al., 2020), whereas their slower cousins were nicknamed ’walkaway
stars’ by (Renzo et al., 2019). Runaway stars are thought to form as
companions to massive stars collapsing in a supernova explosion (Blaauw,
1961), or dynamically ejected from clusters in multiple body encounters
(Poveda et al., 1967). In this paper we consider a third scenario, namely the
destabilisation of an isolated triple. Fig. 12 shows the velocity distribution
of the ejected stars in our simulation relative to the remaining binary
system. The velocity distribution peaks around a few km/s with a tail towards
large velocities. Drifters mainly contribute to velocities below 0.1km$/$s,
and escaping stars above it. Velocities up to several tens km/s are reached,
placing these systems firmly in the regime of runaway stars. It is analogous
to the escape velocities derived from scattering experiments of three-body
interactions with zero total angular momentum (see e.g. Sterzik & Durisen,
1995, 1998), or from N-body simulations for entire star clusters (Fujii &
Portegies Zwart, 2011). As these interactions are typically hierarchical
events in our simulations, the ejected star is typically the original tertiary
of the system.

The magnitude of the escape velocity is related to the orbital velocity before
the destabilisation. For this reason, the fastest escaping stars originate
from the most compact initial systems, and are on the MS during the ejection.
We estimate the terminal velocity $v_{\rm terminal}$ by considering the
average amount of energy that the tertiary extracts from the inner binary
(Appendix C). Assuming Solar masses and a minimum distance between the stars
of 10$R_{\odot}$ (5$R_{\odot}$), the maximum terminal velocity would be 67-157
km/s (95-222 km/s). For more massive stars, the maximum velocity increases
further. For 10$M_{\odot}$ stars and a minimum distance of 10$R_{\odot}$
(50$R_{\odot}$), the maximum terminal velocity reaches 95-222 km/s (213-231
km/s). In comparison, the expected velocities of runaway stars from the binary
(supernova) channel are of the order of 20 km/s (Renzo et al., 2019), as the
orbital velocities decreases substantially during the pre-supernova mass
transfer phase (i.e. orbital widening). As the observed runaway stars are
typically O- and B-type stars666There is a strong observational selection
effect to luminous and therefore massive stars. , and massive stars are
typically in compact orbits with inner periods of a few days ($a_{\rm in}\sim
10-50R_{\odot}$ Sana et al., 2012; Moe & Di Stefano, 2017), we suggest that
destabilised massive triples have the potential to solve the mystery behind
the origin of massive runaway stars.

Lastly, the ejection of stars from destabilised triples (at all velocities),
releases stars into the field that masquerade as ’normal’ single stars.
However, as the ejected stars are formed as part of a multiple, their mass
distribution is likely different from that of a single star formed from a
single cloud. In this way, the destabilisation of triples, is an interesting
formation mechanism for single brown dwarfs , analogous to the hydrodynamical
star formation calculations of unstable multiples (Bate et al., 2002; Delgado-
Donate et al., 2004)

Figure 12:  Histogram of the escape velocities of the escaping stars relative
to the binary for different models. Different colours represent the different
dynamical modelling methods, different line styles represent the different
models for the initial population of triples; OBin (solid), T14 (dashed), E09
(dotted).

###  5.2 MS-MS collisions

Collisions between MS stars frequently occur in dense Galactic clusters (Hills
& Day, 1976; Portegies Zwart et al., 1997; Hurley et al., 2001, 2005; Umbreit
et al., 2008). The rate of such collisions was estimated by Perets & Kratter
(2012) to be of the order of $4\times 10^{-6}(N_{\rm GC}/20)$ per year, where
$N_{\rm GC}\approx 20$ is the number of massive and dense clusters in the
Milky Way. The Galactic rate of MS-MS collisions from destabilised triples is
$(2-2.2)\times 10^{-4}$ per year based on our simulations including SE during
the N-body calculations. If you include collisions in stable and destabilised
triples, the full rate of MS-MS collisions in triples would be even higher.
Stable triples can be driven to contact through Lidov-Kozai cycles and the
related high eccentricities without destabilising the system. In our
simulations with TRES we label these systems highly-eccentric mass-
transferring systems (Fig. 1). Their eccentricities at the onset of the mass
transfer phase ranges from 0.4-1.0 with a median around 0.85-0.9 (see Fig. 7
in Toonen et al., 2020), whereas the eccentricities of the colliding
destabilised triples are more extreme with a median value over 0.996. Based on
the results of Toonen et al. (2020), the Galactic rate would be several
$10^{-3}$ events per year. This is in good agreement with the estimate of He &
Petrovich (2018) for the full MS-MS collision rate in triples, that is about
$10^{-2}$ per year.

Kaib & Raymond (2014) proposed a third channel for producing stellar
collisions. Their channel involves wide binaries that are perturbed by passing
stars and/or the Galactic tide. They estimate a Galactic rate of
$(1.3-10)\times 10^{-4}$ per year. This is similar to the collision rate we
estimate for destabilised triples. However, if we include collisions in stable
triples, then triple evolution is the dominant mechanism to lead to collisions
in the Milky Way.

In a cluster, MS-MS collisions can be recognised as (massive) blue stragglers
(Freitag & Benz, 2005; Lombardi et al., 1995, 1996; Sills et al., 1997, 2001,
2002, 2005). Hydrodynamical simulations show that during the merger little
mass is lost (Lombardi et al., 2002; Dale & Davies, 2006), typically of the
order of a few percent (Glebbeek & Pols, 2008). Some matter may be lost after
the merger in order to spin down the star if the blue straggler is rapidly
rotating, as expected for an off-axis collision (Sills et al., 2001, 2005). If
the remnant maintains its high rotation throughout the main-sequence phase,
rotational mixing extends the lifetime of the blue straggler, as well as
making it bluer and brighter (Sills et al., 2002; Glebbeek & Pols, 2008). In
case of head-on collisions and non- or slowly rotating remnants, substantial
convective regions do not develop during the thermal relaxation phase, and no
significant mixing occurs after the collision (Sills et al., 2001; Glebbeek &
Pols, 2008). Nevertheless, an extension of the lifetime of the remnant is
expected due to hydrodynamical mixing (e.g. Dale & Davies, 2006). The amount
of rejuvenation depends on how much hydrogen is carried into the core during
the merger, which in turn depends on the entropy profiles of the colliding
stars (see e.g. Lombardi et al., 1996, 2002; Glebbeek & Pols, 2008).
Unfortunately, collisional remnants from field triple evolution can not be
detected easily in a Hertzsprung-Russell (HR) diagram as is the case of
cluster blue stragglers, which show two quite distinct tracks (Ferraro et al.,
2009) that indicate a different origin (Portegies Zwart, 2019). Triple star
evolution can, in such environments, also lead to a curious twin blue
stragglers, such as WOCS ID 7782 (Portegies Zwart & Leigh, 2019). In
destabilised triples, the apparent age difference between the collision
product and its outer tertiary star means the system will appear as a non-
coeval binary. Moreover, the rotation of the remnant, in line with the orbital
motion of the former inner binary, is typically inclined with respect to its
orbital motion (Fig. 13).

Figure 13: Distribution of inclinations at the time of the collisions. Green
lines refer to the model where SE is included in the N-body simulations.
Solid, dashed, and dotted line styles represent the initial population of
triples of model OBin, T14 and E09, respectively. The black dashdotdotted line
represents a random distribution in the cosine of the inclination.

The merger event can also be detected directly (Soker & Tylenda, 2003, 2006;
Smith, 2011), as in the case of V1309 Sco (Tylenda et al., 2011), V838 Mon
(Tylenda & Soker, 2006), and even extragalactic events such as M85 OT2006-1
(Kulkarni et al., 2007). An observational link with triples has already been
made, as Kamiński et al. (2021) recently showed that V838 Mon was a binary
merger in a triple system. These intermediate luminosity optical transients
(ILOTs) have a wide range in energies that span the range from classical novae
to supernovae. This is likely related to the large range of possible accretion
masses ( e.g. Soker & Kashi, 2011). In particular, the subset of slowly-
evolving red transients, also known as red luminous novae (RLN), are thought
to originate from stellar mergers in common-envelope events (Ivanova et al.,
2013) with giant donors (see next section Blagorodnova et al., 2017, 2021).

###  5.3 Collisions involving giants

The first estimate of the collision rate in destabilised triples was done by
Perets & Kratter (2012). They recognised the importance of stellar winds to
drive the instability, and therefore the involvement of evolved stars in the
TEDI. They estimate a Galactic rate of collisions involving stellar giants
(mostly AGB stars) of $10^{-4}$ per year. Our estimate of giant-MS collisions
in destabilised triples is about an order of magnitude lower, that is
$(0.4-1.8)\times 10^{-5}$ per year. The difference can likely be attributed to
the evolution of the triples up to the destabilisation, in other words the
initial conditions at the moment of the destabilisation. Whereas they assume a
circular random distribution of inclinations, our systems avoid inclinations
around 90∘. Systems with inclinations close to 90∘ experience strong
gravitational perturbations and eccentricity variations due to three-body
dynamics, and would likely lead to collisions during the dynamically unstable
phase. However, if three-body dynamics is taken into account in the
hierarchical phase, as is the case with TRES, the systems interact prior to
the destabilisation leading to a mass transfer phase or a collision earlier in
the evolution (e.g. on the MS). As noted in Sect. 4.1 for the hierarchical
phase, the far majority of our destabilised triples are in and even born in
the octupole regime, and three-body dynamics is an important contributor to
the formation of the destabilisation. Despite the lower rate compared to
Perets & Kratter (2012), the events are still at least as common as collisions
in Galactic clusters.

Depending on the impact parameter of the collision, we envision three outcomes
(Bailey & Davies, 1999): Either 1) the MS star merges with the giant, or 2) a
compact binary is formed removing the envelope in a common-envelope like event
(for a review on common-envelope evolution see Ivanova et al., 2013), or 3)
the MS star passes through the envelope removing (part of) the envelope on its
way out. The former would lead to a giant star with an anomalously large
envelope mass for its core and an anomalous internal rotation pattern, which
could be detectable with asteroseismology. The latter scenario, may give rise
to a sub-subgiant or a red straggler (Geller et al., 2017a; Leiner et al.,
2017). Like blue stragglers, these stars are found in atypical areas of the
HR-diagram of clusters. But whereas blue stragglers lie bluewards of the turn-
off, sub-subgiants and red stragglers are observed to be redder than normal
main-sequence stars. Sub-subgiants are fainter than normal (sub)giants, and
red stragglers lie redward of the red giant branch (although often considered
to be part of the sub-subgiant family). Leiner et al. (2017) showed that a
rapid partial removal of the envelope mass leads to a rapid decrease in
luminosity of the giant and the formation of a sub-subgiant. However, close
encounters between a binary and a passing star are not frequent enough to
explain the observed systems (Leiner et al., 2017; Geller et al., 2017b), and
so it would be interesting to see if destabilised triples could solve this
mystery. Interestingly, Geller et al. (2017a) found that the number of sub-
subgiants per unit mass increases towards lower-mass clusters, which is
consistent with hierarchical triples being more prevalent in open clusters
than in the dense globular clusters (Leigh & Geller, 2013).

In the second and third case, when most or all of the envelope mass is
stripped from the giant, a planetary nebulae (PNe) may form, assuming that the
matter is ionised sufficiently by the central object. Given the highly
eccentric nature of the collisions and the close proximity of the tertiary
star, the resulting PN can deviate from spherical symmetry, and even axial-
and point-symmetry. The origin of the PN morphologies are sought in magnetic
fields, binary motion and stellar rotation amongst others (e.g. Soker & Hadar,
2002; De Marco, 2009; De Marco & Soker, 2011; Zijlstra, 2015; Jones & Boffin,
2017), however Bear & Soker (2017) estimate that one in six PN are too messy
to be accounted for by the processes mentioned. Soker (2016) estimates that
one in eight PN by triple evolutionary pathways, like the one studied in this
paper. A direct connection between PN and triples was even made for PN Sh 2-71
(Jones et al., 2019).

###  5.4 Collisions involving white dwarfs

After MS-MS collisions, the most typical collision involves a WD. Our Galactic
rate estimates are $(0.72-3.9)\times 10^{-5}$ per year for WD-MS collisions,
$\lesssim 3.4\times 10^{-6}$ per year for collisions between a WD and an
evolved star, and $\lesssim 1.7\times 10^{-6}$ per year for WD-WD collisions.
As the WDs in this channel are post-AGB objects, they are composed of carbon-
oxygen (CO) and/or neon, and have masses above 0.5M⊙.

During a collision between a WD and a MS, the MS experiences strong tidal
torques that can lead to the tidal disruption of the star (for
0.1-1$M_{\odot}$ MS stars), and a significant disruption for more massive
stars (Regev & Shara, 1987). Hydrodynamical simulations show that a massive
disk forms around the WD, liberating the binding energy of the MS (Shara &
Regev, 1986; Regev & Shara, 1987; Soker et al., 1987; Rozyczka et al., 1989).
Michaely & Shara (2021) estimate a bolometric luminosity of
$10^{7}-10^{9}L_{\odot}$, which makes these events similar to kilonova in
brightness. After the merger a giant-like star forms (Rozyczka et al., 1989;
Hurley et al., 2002). Similar to what is discussed in Sect. 5.2, the newly
formed giant likely has a anomalous core-to-envelope mass ratio, internal
rotation pattern, and position in the HR-diagram.

In the case of a collision between two WDs, hydrodynamical simulations have
shown that the shock can trigger nuclear reactions and the formation of a
short-lived R Coronae Borealis star (Webbink, 1984; Lorén-Aguilar et al.,
2010; Schwab, 2019, 2021), or lead to the total disruption of the star in a
supernova-like explosion (Rosswog et al., 2009; Raskin et al., 2009; Kushnir
et al., 2013). Expected luminosities range from supernova Type Ia-like
(Rosswog et al., 2009; Kushnir et al., 2013, and even matching spectra and
lightcurves) to underluminous supernovae (Raskin et al., 2009). Lorén-Aguilar
et al. (2010) also showed that even in the case of a post-collision remnant,
the fallback luminosities are close to $10^{48}$erg/s.

###  5.5 Eta Carinae

An interesting system to mention in the context of collisions and triples is
Eta Carinae (Damineli et al., 1997). Eta Carinae is currently a binary system
with two massive stars in a highly eccentric orbit, and is well known for its
Great Eruption in the middle of the nineteenth century (de Vaucouleurs &
Eggen, 1952). A link with triples was made by Portegies Zwart & van den Heuvel
(2016) who proposed that the Great Eruption was caused by the merger of the
stars in the inner binary due to the Lidov-Kozai cycles imposed by the
tertiary companion (see also Toonen et al., 2016). Later on Smith et al.
(2018) suggested a formation channel involving a destabilised triple. In their
scenario the instability is driven by the expansion of the inner orbit due to
mass transfer. During the dynamically unstable phase, the original secondary
and tertiary collide, causing the Great Eruption. Detailed simulations of the
scenario have been performed by Hirai et al. (2021). These consist of
hydrodynamical simulations of the merger, as well as pure N-body simulations
of the dynamically unstable phase for one triple with 1000 initialisations of
the initial orbital phase of the inner binary. Out of their runs resulting in
a merger, only 8-11% experienced a swap of the companions. While our
simulations do not consider mass-transfer-driven instabilities, they show that
collisions occur between the original primary and secondary (and not the
tertiary) with few exceptions. Despite the similar mass ratio between our
simulations and the proposed scenario for Eta Carinae, we only find collisions
involving the tertiary in 1.6-2.4% of all collisions. This does not change
significantly when focusing only on ((MS,MS),MS) systems (i.e. 1.2-1.5%) or
compact orbits (i.e. 1.8-2.0% for $a_{\rm in}<10^{4}R_{\odot}$). This would
imply that for the scenario to work for every Eta Carinae there should be 50
similar systems that have experienced a collision in the inner binary.

###  5.6 Wide binary formation

The majority of destabilised triples come out of the dynamically unstable
phase as a wide binary system, either through ejecting a star or through a
collision of two of the stars. That the latter is possible, in particular that
the tertiary can remain bound, is demonstrated by OW Geminorum (Eggleton,
2002), and see Eggleton (2013) for three other cases. This binary with a 3.45
yr period and eccentricity of 0.53 consists of two sub-giants with
surprisingly different masses of 5.5$M_{\odot}$ and 3.8$M_{\odot}$ (Griffin &
Duquennoy, 1993; Terrell et al., 2003; Galan et al., 2008). As the MS-lifetime
is a strong function of stellar mass, two single stars with such different
masses should not be on the giant branch simultaneously. Galan et al. (2008)
estimates the lower mass star to be at least twice as old (200 million years)
as its companion. OW Geminorum posed a problem for standard stellar evolution
theory if stellar interactions are not provoked. Indeed, Eggleton (2002)
proposed a triple-star solution; OW Geminorum formed through a merger in a
triple with a compact (2d) inner binary with a $\sim$4$M_{\odot}$ and
$\sim$2$M_{\odot}$ star that was orbited by a $\sim$4$M_{\odot}$ tertiary. If
OW Geminorum formed through the destabilisation of a triple instead, we
envision the following set-up:
\- masses the same as suggested by Eggleton (2002)
\- assuming the outer orbit has not changed much during the unstable phase,
the inner period at destabilisation would have been 50-100 days (Eq. 4). This
is wide enough for tides not to quench the three-body dynamics and (to a
secondary) degree for the stars to evolve and develop winds. Both processes
contribute to the breakdown of the stability.

---
Figure 14: Properties of the (isolated) binaries that are formed from
destabilised triples. Solid, dashed, and dotted line styles represent the
N-body + SE models for the initial population of triples OBin, T14 and E09,
respectively. The black dashed line represents a thermal distribution of
eccentricities.

The typical orbital separations and eccentricities of the remaining binaries
are shown in Fig. 14. The eccentricity distributions are sub-thermal, and most
strongly so for the post-collision binaries. The orbital separations of the
newly found binaries are large, and span a broad range from
$10^{2}-10^{3}R_{\odot}$ upwards. The top panel of Fig. 14 shows that the
adopted initial orbital separation distribution has a strong imprint on the
demographics of the binaries.

As discussed above, post-collision binaries can be recognised as they harbor
peculiar stars (e.g. blue, yellow and red stragglers) or by the apparent age
difference (as for OW Geminorum). Blue stragglers are indeed typically found
in wide binaries. Mathieu & Geller (2009) found that 76% of blue stragglers in
NGC 188 are part of a binary with typical orbits around a $10^{3}$d (or at
least $10^{3}$d as longer periods are hard to detect observationally).
Moreover, the orbits have eccentricities between 0 and 0.8. This is surprising
as the classical mass-transfer scenario in isolated binaries would have
circularised the orbits (but see Bonačić Marinović et al., 2008; Geller &
Mathieu, 2011; Dermine et al., 2013). The eccentricities indicate an origin in
multiple stars. One scenario involves Lidov-Kozai cycles in combination with
tidal friction (Perets & Fabrycky, 2009; Naoz & Fabrycky, 2014). Here we show
that also collisions in destabilised triples contribute in their formation,
albeit in the limit where Roche lobe overflow is neglected.

###  5.7 Compact binary formation

In this paper we focus on isolated, hierarchical stellar triples that become
dynamically unstable as they evolve. We have shown that the duration of the
dynamically unstable phase is much longer than what is generally expected from
binary-single resonant interactions or triples drawn from Plummer spheres.
Although triples on the edge of stability are considered unstable according to
various stability criteria, they tend to be long lived reaching up to
thousands of crossing times. And so, we argue that it is important to take
stellar evolution into account during the unstable phase. Besides stellar
evolution, dissipative processes may also play a role. Examples are stellar
tides, gravitational wave emission or general relativistic precession. These
processes act most efficiently during strong encounters or as a cumulative
effect during repeated close passages (Kaib & Raymond, 2014; Michaely &
Perets, 2016, 2019; Hamers & Thompson, 2019; Michaely & Perets, 2020; Michaely
& Shara, 2021). As a result, we expect dissipative processes to be most
important during democratic encounters (see e.g. Ginat & Perets, 2021). In
general, for ordinary triples in which dissipative processes play an important
role in the dynamics, a compact inner binary forms that is kinematically
decoupled from the tertiary (Mazeh & Shaham, 1979; Kiseleva et al., 1998;
Fabrycky & Tremaine, 2007; Liu et al., 2015; Antonini et al., 2017a; Rodriguez
& Antonini, 2018a; Bataille et al., 2018). Such a compact orbit may eventually
lead to a merger of the two stars in the inner binary. How effective
dissipative processes are for destabilised triples, which is outside of the
scope of this paper, can be verified by using a direct N-body code, such as
TIDYMESS (Boekholt & Correia, In preparation), which combines orbital
evolution with stellar spins and gravitational tides between each pair of
bodies in a self-consistent manner. It would likely lead to the formation of
compact inner binaries at the expense of collisions and stellar escapes. The
further evolution of the system may once again lead to a dynamical
instability. As the stars in the newly-formed compact inner binary evolve,
mass transfer may develop in the inner binary. If the mass transfer is stable
and the orbit widens consequently, another dynamical instability may commence.
On the other hand, if the stars in the compact inner binary are compact
objects, a gravitational wave source can form that is detectable by LISA or
even LIGO/Virgo (e.g. Antognini et al., 2014; Silsbee & Tremaine, 2017;
Antonini et al., 2017b; Hoang et al., 2018; Rodriguez & Antonini, 2018b;
Bonetti et al., 2019). If any eccentricity remains in the orbit at the time of
detection, the system would stand out compared to those formed through
isolated binary evolution. The compact binary with compact objects will
eventually merge and give rise to electromagnetic transients such as supernova
type Ia, gamma-ray bursts or kilonova (Webbink, 1984; Metzger et al., 2010;
Pakmor et al., 2012; Shen et al., 2018).

##  6 Discussion

In this study we have constructed three different models for the initial
population of triples that differ from one another based on their mass, mass
ratio and period distribution. For all models we have assumed that the mutual
inclinations of the orbits follow a circular uniform distribution. The
observations of Tokovinin (2017) support this assumption for outer projected
separations $\gtrsim 10^{5}R_{\odot}$. For low-mass triples with outer
projected separations $\lesssim 10^{4}R_{\odot}$ there is more alignment.
Which destabilised triples would be affected by this and how? The mentioned
range of affected orbital separations corresponds to the regime where the
destabilised triples mainly consist of three MS stars (Fig. 3), which become
unstable due to three body dynamics (i.e. the octupole term). Whereas the
amplitude of Lidov-Kozai cycles are strongly dependent on the initial
inclination (limiting the cycles to $39.2-140.8^{\circ}$ in the test-particle
approximation), the eccentric Lidov-Kozai mechanism is accessible for a larger
range of inclinations (Naoz & Fabrycky, 2014; Naoz, 2016). Moreover, the
inclination distribution of the destabilised triples in our simulations avoid
inclinations around 90∘ in the first place, as these typically lead to mass
transfer (Toonen et al., 2020, see e.g. their Fig.15). And so, if in real life
compact triples have more alignment, it would mean our synthetic compact
destabilised triples would have retrograde orbits less frequently than
currently modelled. These are the systems that predominantly experience
collisions, where as their prograde counterparts predominantly dissolve into a
binary and an escaping star.

###  6.1 Limitations of the model

The simulations performed in this paper consist of two parts. The hierarchical
phase is simulated with TRES, which is based on the double averaged method, up
to the stability limit. The subsequent evolution is simulated with the N-body
approach. In general the stability limit that we adopted (Eq. 4, Mardling &
Aarseth, 1999; Aarseth & Mardling, 2001) is considered to be
conservative777That is, it may take longer than a Hubble time for the
instability to manifest itself (Mylläri et al., 2018; He & Petrovich, 2018).
This is further corroborated by the fractions of bound systems left in our
various N-body simulations and the large number of crossing times simulated
(Fig. 4). However, it is good that we transition from TRES to the N-body
simulations, as we notice the break down of the double-averaging method in
sections 4.2.4-4.2.6. This behaviour leading to variations in the outer
orbital separations, as well as higher maximum eccentricities we would not
have caught with the double-averaging approach with quadropole and octupole
terms (see e.g. Antonini & Perets, 2012; Katz & Dong, 2012; Luo et al., 2016;
Grishin et al., 2018; Bhaskar et al., 2021b).

This leads one to wonder: what if we would have stopped the TRES simulations
earlier? The experiment in appendix B shines light on this matter, in
particular model B & D. These models show the difference in outcome for pure
N-body simulations of retrograde triples on the stability limit (model B) or
slightly above it (model D). Whereas over 60% of triples in model B experience
a catastrophic event (i.e. collision or ejection), this is only the case for
2% in model D, for which the far majority of triples is long-lived and remain
bound (Tbl. 4). This behaviour confirms that the triples are indeed not
dynamically unstable in model D. If the double averaged method including the
octupole term underestimates the maximum eccentricity for a given triple as it
approaches the stability limit (Antonini & Perets, 2012; Katz & Dong, 2012;
Luo et al., 2016; Grishin et al., 2018; Bhaskar et al., 2021b), then the
triple would have destabilised earlier in the evolution compared to what is
assumed in this work. This is mostly important for the ((MS,MS),MS) systems,
as the instability for systems with giants are largely driven by their stellar
winds. However, if the destabilisation of a ((MS,MS),MS) system occurred
earlier in its evolution, the outcome of the dynamically unstable phase would
not be very different.

###  6.2 Limitations of the direct model

Our direct integrations include both the orbital and the mass-radius evolution
of the three stars. These two physical ingredients cause the destabilisation
of hierarchical triple star systems. When two stars overlap, we consider them
to have collided, that is the sticky sphere approximation. Our direct model
has some limitations, as we do not include relativistic effects, tides, or
mass transfer during close passages within the Roche limit. These effects
could potentially play a key role in the formation of close binaries and other
phenomena.

The interplay between von Zeipel-Lidov-Kozai cycles and relativistic
precession can lead to various eccentricity behaviours. If the time scale for
relativistic precession becomes significantly shorter than the time scale for
von Zeipel-Lidov-Kozai cycles, then it is possible for the maximum attained
eccentricity to be quenched. However, in the case when the two time scales
resonate, eccentricity excitation can be re-triggered (Ford et al., 2000a;
Naoz et al., 2013; Will, 2017; Lim & Rodriguez, 2020; Hansen & Naoz, 2020).
Depending on the radii of the objects in the inner binary, a change in the
maximum attained eccentricity could directly influence our measured collision
rate. However, for triples on the edge of stability, for which
$a_{\rm{out}}\sim 3a_{\rm{in}}$, the inner orbit is strongly perturbed by the
tertiary star. On the other hand, the relativistic precession time scale of
the inner orbit does not depend on the outer orbit. Hence, for destabilised
triples in which all orbital elements of the inner and outer binary vary in
time, it is unclear whether relativistic precession acts to quench
eccentricity oscillations in the inner binary, or rather acts as a
perturbation to a chaotic system, leading to loss of memory of its initial
condition. The latter scenario on the effect of relativity in a non-
hierarchical, chaotic triple system was recently studied by Boekholt et al.
(2021). They find that relativistic perturbations to the Newtonian force, grow
exponentially due to the exponential sensitivity in the chaotic three-body
problem. This results in a sensitive dependence in the lifetime of the triple,
but the statistical outcome remains unchanged, that is a dissolution or head-
on collision. Only when the velocities start to approach the speed of light,
do gravitational wave captures become prominent outcomes, which mostly applies
to compact stellar remnants.

Tidal dissipation tends to circularise eccentric binaries into an orbit whose
semi-major axis is about twice the initial pericentre separation. The damping
of the eccentricity prevents direct collisions at pericentre in the sticky
sphere approximation, and instead leads to the production of close binaries,
or compact triples if the tertiary remains. On the other hand, for
destabilised triple systems which lose their hierarchy and become democratic,
it is still uncertain how tides affect the outcome. The equilibrium tide model
used in our secular calculations, assumes that the tidal lag is always small
compared to the orbital time scale (e.g. Mignard, 1979; Hut, 1981). For high
eccentricity, this will not always hold during close pericentre passages,
where the orbital time scale can be very short. In this case, non-linear tides
is a better description (such as the Maxwell model, e.g. Correia et al.
(2014)). Furthermore, if the orbital time scale resonates with internal
oscillation modes of the star as in dynamical tide models (e.g. Fuller & Lai,
2012), this triggers a strong enhancement of the tidal dissipation. This would
lead us to expect that the equilibrium tide model underestimates the amount of
tidal dissipation during very close pericentre passages. As a consequence,
high eccentricities in the inner orbit of a triple are less effectively
damped, thereby prolonging the possibility of a collision or dynamical
dissolution. Our estimated fractions for these two outcomes can thus be
considered as upper limits, and a benchmark test for future studies with more
sophisticated tidal models.

Another limitation of our numerical setup is the simplified handling of close
encounters and collisions. We consider two stars to have collided if their
mutual distance is less than the sum of their radii, that is they overlap. If
their separation is only slightly larger however, tidal effects and Roche lobe
overflow events occur instead. For two solar mass stars in a circular orbit,
this happens within a separation of about 1.25 times the sum of their radii,
while for a giant and a white dwarf a separation of up to 2.5 times the sum of
their radii suffices. In our simulations these effects are not taken into
account, which could potentially introduce a bias. This can be tested by
including physical models for tidal effects and eccentric mass transfer
(episodes), whose influence on the eccentricity fluctuations of the inner
binary is to be compared to the chaotic fluctuations caused by the tertiary.
The inner binaries’ semi-major axis could also be affected if orbital energy
is dissipated by tides or if there is mass transfer and/or mass loss. However,
(episodic) mass transfer events are obviously challenging to model efficiently
in a population synthesis study, but apart from that, the workings of mass
transfer in very eccentric inner binaries is still poorly understood. The
inclusion of relativity, tides, and other effects in the direct simulation of
destabilised triples is left for a future study, whose benchmark will be the
current study. Here, our aim is to determine the influence of different
stellar evolution models on the outcome of destabilised triple systems.

##  7 Conclusions

We have presented a numerical study on the evolution of destabilised,
hierarchical triple-star systems, that is ‘triples on the edge’. The current
study follows-up on Toonen et al. (2020) who studied the typical evolution of
hierarchical triples, and focuses on those triples that start out in a stable
configuration (Aarseth & Mardling, 2001), but become dynamically unstable due
to stellar wind mass loss and/or strong dynamical perturbations between the
inner and outer orbit. The initial hierarchical phase is modelled in the
secular regime with the multi-physics code TRES, while the unstable phase is
treated by direct N-body integration. To determine the influence of stellar
evolutionary effects during the unstable phase, we consider 1) a model with
non-evolving stars, 2) a model with constant mass loss rates, and 3) a
realistic model of stellar evolution following the code SeBa (Portegies Zwart
& Verbunt, 1996; Toonen et al., 2012). A comprehensive overview of our main
results regarding the evolution of destabilised triples is given below.

  * •

Hierarchy: 2-4% of all low and intermediate mass triples become destabilised.
Of these destabilised triples, 80-83% preserve their initial hierarchy
throughout the evolution in the case when stellar winds are neglected, while
54-69% do so when stellar evolution is included. This is in contradiction with
the commonly adopted picture that unstable triples always experience a
chaotic, democratic resonant interaction.

  * •

Lifetime: The duration of the unstable phase has a median of $10^{3-4}$
crossing times, reaching up to millions. This is in contradiction with the
assumption that the lifetime of unstable triples is relatively short, so that
long-term effects can be neglected. The extended duration allows for other
physical effects, such as stellar evolution, to play a crucial role in the
evolution.

  * •

Lifetime: In fact for a number of systems, the lifetime is sufficiently long,
in excess of $10^{6}-10^{7}$ crossing times, such that it exceeds our maximum
CPU time of 24 hours. It comprises 2-4% of destabilised triples for model
OBin, 9-10% for model T14, and 5-6% for model E09. Our rate estimates can
therefore be regarded as lower limits, while for the most common outcomes of
destabilised triple evolution they accurate to within the order of magnitude.

  * •

Ejections: The most probable outcome is a dissolution of the triple into a
single star and binary: 42-45% leads to unbound escapers, and 12-22% to slow
drifters.

  * •

Ejections: Two ejection modes are: 1) the commonly known democratic ejection
during which the initial hierarchy is lost and the lightest body tends to
escape (57-70% for escapers, 24-44% for drifters), and 2) hierarchical
slingshot of the tertiary irrespective of its mass (30-43% for escapers,
56-76% for drifters).

  * •

Runaway stars: Our simulations produce runaway and walkaway stars with speeds
up to several tens km/s. In appendix C, we calculate a maximum velocity of one
to a few hundred km/s. The main factor determining the ejection speed is the
compactness (or total energy) of the initial triple. We suggest that
destabilised massive triples have the potential to solve the mystery behind
the origin of the observed (massive) runaway stars.

  * •

Collisions: Collisions are common and occur in 13-24% of destabilised triples.
Almost all of them occur between the original components of the inner binary,
only in 1.6-2.4% is the tertiary star involved. This is consistent with the
fact that the initial hierarchy is preserved until the moment of the collision
(88-95%), and contradicts the idea that collisions during democratic
encounters dominate (only 5-12%).

  * •

Collisions: Collisions typically involve main-sequence stars (77-94%). This is
in contradiction with the expectation that giant stars are usually involved,
since due to their enhanced radii they form large targets during democratic
resonant encounters. The relatively short time-frame of the giant phase, and
the random walk of the inner eccentricity, result in collisions being
effectively avoided. Once the giant has evolved to a much smaller white dwarf,
collisions can still take place due to the much longer time frame afterwards.
Indeed, we find a higher fraction of collisions between a white dwarf and
main-sequence star (4-14%) than between a giant and main-sequence star (2-6%).

  * •

Collision rates: We estimate the following Galactic event rates for
collisions, which originated from destabilised, hierarchical triple systems:
$(2-2.2)\times 10^{-4}\,\rm{yr^{-1}}$ for MS-MS stars, $(0.72-3.9)\times
10^{-5}\,\rm{yr^{-1}}$ for WD-MS stars, $(0.4-1.8)\times
10^{-5}\,\rm{yr^{-1}}$ for G-MS stars, $3.4\times 10^{-6}\,\rm{yr^{-1}}$ for
G-WD stars, and $1.7\times 10^{-6}\,\rm{yr^{-1}}$ for WD-WD pairs. When
including collisions in stable triples, we find that triple evolution is the
dominant mechanism for stellar collisions in the Milky Way.

The results listed above complement previous results on the dynamical
evolution of unstable triples from binary-single scattering experiments.
Destabilised, hierarchical triple systems can evolve not only along the
democratic track, but also along the hierarchical track. If the evolutionary
path takes the triple through a democratic resonance, then the initial
hierarchy is lost and the main trends found in extensive binary-single
scattering experiments in the literature apply. However, a large fraction of
destabilised hierarchical triples preserve their hierarchy, resulting in
different ejection and collision scenarios. When studying triples on the edge
of stability, both these pathways are prevalent and affect the outcome
statistics.

Based on our statistical analysis, we present two rules of thumb. These can
provide a general prediction for the fate of a specific triple, based on the
initial inclination of the system. Specifically, they concern whether a triple
will destabilise along the hierarchical or the democratic track. Obviously
detailed and accurate predictions are better made with the N-body approach
described in Sect. 3, but -taken with care- the rules of thumb provide a new
and fast insight in the problem of destabilised triples.

  * •

First rule of thumb: prograde triples888See footnote 5 tend to evolve towards
hierarchical collisions and ejections, that is the tertiary star is ejected
irrespective of its mass.

  * •

Second rule of thumb: triples with retrograde orbits tend to evolve towards
democratic encounters and loss of initial hierarchy. This commonly leads to
ejection of the lowest-mass body. However, compact retrograde systems
experience hierarchical collisions.

These statistical rules indicate that the inclination plays a key role in the
energy exchange between the inner and outer orbit of the triple. Since these
are triples on the edge of stability, the pericentre distance of the outer
orbit is comparable to the inner binary separation. As a consequence, the
inner binaries’ orbital speed and the tertiary’s pericentre speed can be
comparable in magnitude. Since the two velocity vectors point in approximately
the same direction for prograde triples, the relative velocity will be small,
which in turn maximises the interaction time in which energy can be exchanged.
In the retrograde case, the velocity vectors point in opposite directions, and
the interaction time is minimzed. Based on this simplified analysis, we thus
expect tertiaries to be able to escape through a ‘hierarchical slingshot’
mostly in prograde triples. Further studies are required to reveal the role of
inclination, and other potential indicators, in shaping the dynamical
evolution of destabilised, hierarchical triple star systems and corresponding
observables.

###### Acknowledgements.

We thank Nathan Leigh for the discussions on stellar collisions. This project
was supported by funds from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program under grant
agreement No 638435 (GalNUC). ST acknowledge support from the Netherlands
Research Council NWO (VENI 639.041.645 grants and VIDI 203.061 grants). TB
acknowledges support from ENGAGE SKA RI, grant POCI-01-0145-FEDER-022217,
funded by COMPETE 2020 and FCT, Portugal. This work is funded by FCT/MEC
through national funds and when applicable co-funded by FEDER – PT2020
partnership agreement under the project UID/EEA/50008/2019. The calculations
were performed using the LGMII (NWO grant #621.016.701).
Data availability. A data reproduction package for this paper can be found at
the following DOI: 10.5281/zenodo.6319686

## References

  * Aarseth (2004) Aarseth, S. J. 2004, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 21, Revista Mexicana de Astronomia y Astrofisica Conference Series, ed. C. Allen & C. Scarfe, 156–162
  * Aarseth & Mardling (2001) Aarseth, S. J. & Mardling, R. A. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 229, Evolution of Binary and Multiple Star Systems, ed. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, 77
  * Abt (1983) Abt, H. A. 1983, ARA&A, 21, 343
  * Anosova (1986) Anosova, J. P. 1986, Ap&SS, 124, 217
  * Antognini et al. (2014) Antognini, J. M., Shappee, B. J., Thompson, T. A., & Amaro-Seoane, P. 2014, MNRAS, 439, 1079
  * Antognini (2015) Antognini, J. M. O. 2015, MNRAS, 452, 3610
  * Antonini et al. (2016) Antonini, F., Chatterjee, S., Rodriguez, C. L., et al. 2016, ApJ, 816, 65
  * Antonini & Perets (2012) Antonini, F. & Perets, H. B. 2012, ApJ, 757, 27
  * Antonini et al. (2017a) Antonini, F., Toonen, S., & Hamers, A. S. 2017a, ApJ, 841, 77
  * Antonini et al. (2017b) Antonini, F., Toonen, S., & Hamers, A. S. 2017b, ApJ, 841, 77
  * Bailey & Davies (1999) Bailey, V. C. & Davies, M. B. 1999, MNRAS, 308, 257
  * Baraffe et al. (2002) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563
  * Bataille et al. (2018) Bataille, M., Libert, A.-S., & Correia, A. C. M. 2018, MNRAS, 479, 4749
  * Bate et al. (2002) Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 332, L65
  * Bear & Soker (2017) Bear, E. & Soker, N. 2017, ApJ, 837, L10
  * Becker et al. (2012) Becker, W., Prinz, T., Winkler, P. F., & Petre, R. 2012, ApJ, 755, 141
  * Bhaskar et al. (2021a) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2021a, AJ, 161, 48
  * Bhaskar et al. (2021b) Bhaskar, H., Li, G., Hadden, S., Payne, M. J., & Holman, M. J. 2021b, AJ, 161, 48
  * Bischoff et al. (2020) Bischoff, R., Mugrauer, M., Torres, G., et al. 2020, Astronomische Nachrichten, 341, 908
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Blaes et al. (2002a) Blaes, O., Lee, M. H., & Socrates, A. 2002a, ApJ, 578, 775
  * Blaes et al. (2002b) Blaes, O., Lee, M. H., & Socrates, A. 2002b, ApJ, 578, 775
  * Blagorodnova et al. (2021) Blagorodnova, N., Klencki, J., Pejcha, O., et al. 2021, A&A, 653, A134
  * Blagorodnova et al. (2017) Blagorodnova, N., Kotak, R., Polshaw, J., et al. 2017, ApJ, 834, 107
  * Boekholt & Portegies Zwart (2015) Boekholt, T. & Portegies Zwart, S. 2015, Computational Astrophysics and Cosmology, 2, 2
  * Boekholt & Correia (In preparation) Boekholt, T. C. N. & Correia, A., C. M. In preparation
  * Boekholt et al. (2021) Boekholt, T. C. N., Moerman, A., & Portegies Zwart, S. F. 2021, Phys. Rev. D, 104, 083020
  * Boekholt et al. (2020) Boekholt, T. C. N., Portegies Zwart, S. F., & Valtonen, M. 2020, MNRAS, 493, 3932
  * Boekholt & Toonen (In preparation) Boekholt, T. C. N. & Toonen, S. In preparation
  * Bonačić Marinović et al. (2008) Bonačić Marinović, A. A., Glebbeek, E., & Pols, O. R. 2008, A&A, 480, 797
  * Bonetti et al. (2019) Bonetti, M., Sesana, A., Haardt, F., Barausse, E., & Colpi, M. 2019, MNRAS, 486, 4044
  * Boubert & Evans (2018) Boubert, D. & Evans, N. W. 2018, MNRAS, 477, 5261
  * Chatterjee et al. (2005) Chatterjee, S., Vlemmings, W. H. T., Brisken, W. F., et al. 2005, ApJ, 630, L61
  * Cordes et al. (1993) Cordes, J. M., Romani, R. W., & Lundgren, S. C. 1993, Nature, 362, 133
  * Correia et al. (2014) Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50
  * Dale & Davies (2006) Dale, J. E. & Davies, M. B. 2006, MNRAS, 366, 1424
  * Damineli et al. (1997) Damineli, A., Conti, P. S., & Lopes, D. F. 1997, New A, 2, 107
  * De Donder et al. (1997) De Donder, E., Vanbeveren, D., & van Bever, J. 1997, A&A, 318, 812
  * De Marco (2009) De Marco, O. 2009, PASP, 121, 316
  * De Marco & Soker (2011) De Marco, O. & Soker, N. 2011, PASP, 123, 402
  * de Vaucouleurs & Eggen (1952) de Vaucouleurs, G. & Eggen, O. J. 1952, PASP, 64, 185
  * Delgado-Donate et al. (2004) Delgado-Donate, E. J., Clarke, C. J., Bate, M. R., & Hodgkin, S. T. 2004, MNRAS, 351, 617
  * Dermine et al. (2013) Dermine, T., Izzard, R. G., Jorissen, A., & Van Winckel, H. 2013, A&A, 551, A50
  * Dray et al. (2005) Dray, L. M., Dale, J. E., Beer, M. E., Napiwotzki, R., & King, A. R. 2005, MNRAS, 364, 59
  * Duchêne & Kraus (2013) Duchêne, G. & Kraus, A. 2013, ARA&A, 51, 269
  * Duquennoy & Mayor (1991) Duquennoy, A. & Mayor, M. 1991, A&A, 248, 485
  * Eggleton (2013) Eggleton, P. 2013, in Giants of Eclipse, Vol. 45, 10201
  * Eggleton & Kiseleva (1995) Eggleton, P. & Kiseleva, L. 1995, ApJ, 455, 640
  * Eggleton (2002) Eggleton, P. P. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 279, Exotic Stars as Challenges to Evolution, ed. C. A. Tout & W. van Hamme, 37
  * Eggleton (2009) Eggleton, P. P. 2009, MNRAS, 399, 1471
  * Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, MNRAS, 414, 3501
  * Evans (2011) Evans, N. R. 2011, Bulletin de la Societe Royale des Sciences de Liege, 80, 663
  * Fabrycky & Tremaine (2007) Fabrycky, D. & Tremaine, S. 2007, ApJ, 669, 1298
  * Ferraro et al. (2009) Ferraro, F. R., Beccari, G., Dalessandro, E., et al. 2009, Nature, 462, 1028
  * Ford et al. (2000a) Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000a, ApJ, 535, 385
  * Ford et al. (2000b) Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000b, ApJ, 535, 385
  * Fragione et al. (2020) Fragione, G., Martinez, M. A. S., Kremer, K., et al. 2020, ApJ, 900, 16
  * Freire et al. (2011) Freire, P. C. C., Bassa, C. G., Wex, N., et al. 2011, MNRAS, 412, 2763
  * Freitag & Benz (2005) Freitag, M. & Benz, W. 2005, MNRAS, 358, 1133
  * Fujii & Portegies Zwart (2011) Fujii, M. S. & Portegies Zwart, S. 2011, Science, 334, 1380
  * Fuller & Lai (2012) Fuller, J. & Lai, D. 2012, MNRAS, 420, 3126
  * Galan et al. (2008) Galan, C., Mikolajewski, M., Tomov, T., et al. 2008, The Observatory, 128, 298
  * Gao et al. (2019) Gao, Y., Li, J., & Jia, S. 2019, MNRAS, 487, 3178
  * Geller et al. (2017a) Geller, A. M., Leiner, E. M., Bellini, A., et al. 2017a, ApJ, 840, 66
  * Geller et al. (2017b) Geller, A. M., Leiner, E. M., Chatterjee, S., et al. 2017b, ApJ, 842, 1
  * Geller & Mathieu (2011) Geller, A. M. & Mathieu, R. D. 2011, Nature, 478, 356
  * Georgakarakos (2008) Georgakarakos, N. 2008, Celestial Mechanics and Dynamical Astronomy, 100, 151
  * Ginat & Perets (2021) Ginat, Y. B. & Perets, H. B. 2021, Physical Review X, 11, 031020
  * Glebbeek & Pols (2008) Glebbeek, E. & Pols, O. R. 2008, A&A, 488, 1017
  * Griffin & Duquennoy (1993) Griffin, R. F. & Duquennoy, A. 1993, The Observatory, 113, 53
  * Grishin et al. (2018) Grishin, E., Perets, H. B., & Fragione, G. 2018, MNRAS, 481, 4907
  * Gunn & Ostriker (1970) Gunn, J. E. & Ostriker, J. P. 1970, ApJ, 160, 979
  * Hamers (2020) Hamers, A. S. 2020, MNRAS, 494, 5298
  * Hamers et al. (2013) Hamers, A. S., Pols, O. R., Claeys, J. S. W., & Nelemans, G. 2013, MNRAS, 430, 2262
  * Hamers & Thompson (2019) Hamers, A. S. & Thompson, T. A. 2019, ApJ, 882, 24
  * Hansen & Naoz (2020) Hansen, B. M. S. & Naoz, S. 2020, MNRAS, 499, 1682
  * Harrington (1968a) Harrington, R. S. 1968a, AJ, 73, 190
  * Harrington (1968b) Harrington, R. S. 1968b, AJ, 73, 190
  * Harrington (1969) Harrington, R. S. 1969, Celestial Mechanics, 1, 200
  * He & Petrovich (2018) He, M. Y. & Petrovich, C. 2018, MNRAS, 474, 20
  * Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
  * Hills (1975) Hills, J. G. 1975, AJ, 80, 809
  * Hills (1983) Hills, J. G. 1983, ApJ, 267, 322
  * Hills & Day (1976) Hills, J. G. & Day, C. A. 1976, Astrophys. Lett., 17, 87
  * Hirai et al. (2021) Hirai, R., Podsiadlowski, P., Owocki, S. P., Schneider, F. R. N., & Smith, N. 2021, MNRAS, 503, 4276
  * Hoang et al. (2018) Hoang, B.-M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, ApJ, 856, 140
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Hurley et al. (2005) Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
  * Hurley et al. (2001) Hurley, J. R., Tout, C. A., Aarseth, S. J., & Pols, O. R. 2001, MNRAS, 323, 630
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Hut & Bahcall (1983a) Hut, P. & Bahcall, J. N. 1983a, ApJ, 268, 319
  * Hut & Bahcall (1983b) Hut, P. & Bahcall, J. N. 1983b, ApJ, 268, 319
  * Iben & Tutukov (1999) Iben, Icko, J. & Tutukov, A. V. 1999, ApJ, 511, 324
  * Igoshev (2020) Igoshev, A. P. 2020, MNRAS, 494, 3663
  * Ivanova et al. (2013) Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59
  * Jiang & Tremaine (2010) Jiang, Y.-F. & Tremaine, S. 2010, MNRAS, 401, 977
  * Jones & Boffin (2017) Jones, D. & Boffin, H. M. J. 2017, Nature Astronomy, 1, 0117
  * Jones et al. (2019) Jones, D., Pejcha, O., & Corradi, R. L. M. 2019, MNRAS, 489, 2195
  * Kaib & Raymond (2014) Kaib, N. A. & Raymond, S. N. 2014, ApJ, 782, 60
  * Kamiński et al. (2021) Kamiński, T., Tylenda, R., Kiljan, A., et al. 2021, A&A, 655, A32
  * Katz & Dong (2012) Katz, B. & Dong, S. 2012, ArXiv e-prints [arXiv:1211.4584]
  * Katz et al. (2011) Katz, B., Dong, S., & Malhotra, R. 2011, Physical Review Letters, 107, 181101
  * Kimpson et al. (2016) Kimpson, T. O., Spera, M., Mapelli, M., & Ziosi, B. M. 2016, MNRAS, 463, 2443
  * Kinoshita & Nakai (1999) Kinoshita, H. & Nakai, H. 1999, Celestial Mechanics and Dynamical Astronomy, 75, 125
  * Kiseleva et al. (1998) Kiseleva, L. G., Eggleton, P. P., & Mikkola, S. 1998, MNRAS, 300, 292
  * Kiseleva et al. (1994) Kiseleva, L. G., Eggleton, P. P., & Orlov, V. V. 1994, MNRAS, 270, 936
  * Kol (2021) Kol, B. 2021, Celestial Mechanics and Dynamical Astronomy, 133, 17
  * Kouwenhoven et al. (2007) Kouwenhoven, M. B. N., Brown, A. G. A., Portegies Zwart, S. F., & Kaper, L. 2007, A&A, 474, 77
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Kroupa et al. (1993) Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545
  * Kulkarni et al. (2007) Kulkarni, S. R., Ofek, E. O., Rau, A., et al. 2007, Nature, 447, 458
  * Kushnir et al. (2013) Kushnir, D., Katz, B., Dong, S., Livne, E., & Fernández, R. 2013, ApJ, 778, L37
  * Leigh & Geller (2012) Leigh, N. & Geller, A. M. 2012, MNRAS, 425, 2369
  * Leigh & Geller (2013) Leigh, N. W. C. & Geller, A. M. 2013, MNRAS, 432, 2474
  * Leiner et al. (2017) Leiner, E., Mathieu, R. D., & Geller, A. M. 2017, ApJ, 840, 67
  * Li et al. (2014) Li, G., Naoz, S., Holman, M., & Loeb, A. 2014, ApJ, 791, 86
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
  * Lim & Rodriguez (2020) Lim, H. & Rodriguez, C. L. 2020, Phys. Rev. D, 102, 064033
  * Lithwick & Naoz (2011) Lithwick, Y. & Naoz, S. 2011, ApJ, 742, 94
  * Liu et al. (2015) Liu, B., Muñoz, D. J., & Lai, D. 2015, MNRAS, 447, 747
  * Lombardi et al. (1995) Lombardi, James, C. J., Rasio, F. A., & Shapiro, S. L. 1995, ApJ, 445, L117
  * Lombardi et al. (1996) Lombardi, J. C. J., Rasio, F. A., & Shapiro, S. L. 1996, ApJ, 468, 797
  * Lombardi et al. (2002) Lombardi, J. C. J., Warren, J. S., Rasio, F. A., Sills, A., & Warren, A. R. 2002, ApJ, 568, 939
  * Lorén-Aguilar et al. (2010) Lorén-Aguilar, P., Isern, J., & García-Berro, E. 2010, MNRAS, 406, 2749
  * Lu & Naoz (2019) Lu, C. X. & Naoz, S. 2019, MNRAS, 484, 1506
  * Luo et al. (2016) Luo, L., Katz, B., & Dong, S. 2016, MNRAS, 458, 3060
  * Lyne & Lorimer (1994) Lyne, A. G. & Lorimer, D. R. 1994, Nature, 369, 127
  * Makino (1991) Makino, J. 1991, ApJ, 369, 200
  * Mardling & Aarseth (1999) Mardling, R. & Aarseth, S. 1999, in NATO Advanced Science Institutes (ASI) Series C, Vol. 522, NATO Advanced Science Institutes (ASI) Series C, ed. B. A. Steves & A. E. Roy, 385
  * Mardling (2001) Mardling, R. A. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 229, Evolution of Binary and Multiple Star Systems, ed. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, 101
  * Mathieu & Geller (2009) Mathieu, R. D. & Geller, A. M. 2009, Nature, 462, 1032
  * Mazeh & Shaham (1979) Mazeh, T. & Shaham, J. 1979, A&A, 77, 145
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650
  * Michaely & Perets (2014) Michaely, E. & Perets, H. B. 2014, ApJ, 794, 122
  * Michaely & Perets (2016) Michaely, E. & Perets, H. B. 2016, MNRAS, 458, 4188
  * Michaely & Perets (2019) Michaely, E. & Perets, H. B. 2019, ApJ, 887, L36
  * Michaely & Perets (2020) Michaely, E. & Perets, H. B. 2020, MNRAS, 498, 4924
  * Michaely & Shara (2021) Michaely, E. & Shara, M. M. 2021, MNRAS, 502, 4540
  * Mignard (1979) Mignard, F. 1979, Moon and Planets, 20, 301
  * Mikkola (1988) Mikkola, S. 1988, On the Effects of Unequal Masses in the Statistics of Three-And Four-Body Interactions, ed. M. J. Valtonen, 261
  * Mikkola & Tanikawa (2007) Mikkola, S. & Tanikawa, K. 2007, MNRAS, 379, L21
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Moe & Kratter (2018) Moe, M. & Kratter, K. M. 2018, ApJ, 854, 44
  * Monaghan (1976) Monaghan, J. J. 1976, MNRAS, 176, 63
  * Montgomery (1998) Montgomery, R. 1998, Nonlinearity, 11, 363
  * Mylläri et al. (2018) Mylläri, A., Valtonen, M., Pasechnik, A., & Mikkola, S. 2018, MNRAS, 476, 830
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441
  * Naoz & Fabrycky (2014) Naoz, S. & Fabrycky, D. C. 2014, ApJ, 793, 137
  * Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155
  * Naoz et al. (2016) Naoz, S., Fragos, T., Geller, A., Stephan, A. P., & Rasio, F. A. 2016, ApJ, 822, L24
  * Orlov et al. (2010) Orlov, V. V., Rubinov, A. V., & Shevchenko, I. I. 2010, MNRAS, 408, 1623
  * Pakmor et al. (2012) Pakmor, R., Kromer, M., Taubenberger, S., et al. 2012, ApJ, 747, L10
  * Perets & Fabrycky (2009) Perets, H. B. & Fabrycky, D. C. 2009, ApJ, 697, 1048
  * Perets & Kratter (2012) Perets, H. B. & Kratter, K. M. 2012, ApJ, 760, 99
  * Pijloo et al. (2012) Pijloo, J. T., Caputo, D. P., & Portegies Zwart, S. F. 2012, MNRAS, 424, 2914
  * Portegies Zwart (2019) Portegies Zwart, S. 2019, A&A, 621, L10
  * Portegies Zwart & Leigh (2019) Portegies Zwart, S. & Leigh, N. W. C. 2019, ApJ, 876, L33
  * Portegies Zwart & McMillan (2018) Portegies Zwart, S. & McMillan, S. 2018, Astrophysical Recipes; The art of AMUSE
  * Portegies Zwart et al. (2011) Portegies Zwart, S., van den Heuvel, E. P. J., van Leeuwen, J., & Nelemans, G. 2011, ApJ, 734, 55
  * Portegies Zwart (2000) Portegies Zwart, S. F. 2000, ApJ, 544, 437
  * Portegies Zwart et al. (1997) Portegies Zwart, S. F., Hut, P., McMillan, S. L. W., & Verbunt, F. 1997, A&A, 328, 143
  * Portegies Zwart & van den Heuvel (2016) Portegies Zwart, S. F. & van den Heuvel, E. P. J. 2016, MNRAS, 456, 3401
  * Portegies Zwart & Verbunt (1996) Portegies Zwart, S. F. & Verbunt, F. 1996, A&A, 309, 179
  * Poveda et al. (1967) Poveda, A., Ruiz, J., & Allen, C. 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4, 86
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
  * Raskin et al. (2009) Raskin, C., Timmes, F. X., Scannapieco, E., Diehl, S., & Fryer, C. 2009, MNRAS, 399, L156
  * Regev & Shara (1987) Regev, O. & Shara, M. M. 1987, MNRAS, 227, 967
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66
  * Rodriguez & Antonini (2018a) Rodriguez, C. L. & Antonini, F. 2018a, ApJ, 863, 7
  * Rodriguez & Antonini (2018b) Rodriguez, C. L. & Antonini, F. 2018b, ApJ, 863, 7
  * Rosswog et al. (2009) Rosswog, S., Kasen, D., Guillochon, J., & Ramirez-Ruiz, E. 2009, ApJ, 705, L128
  * Rozyczka et al. (1989) Rozyczka, M., Yorke, H. W., Bodenheimer, P., Mueller, E., & Hashimoto, M. 1989, A&A, 208, 69
  * Safarzadeh et al. (2020) Safarzadeh, M., Hamers, A. S., Loeb, A., & Berger, E. 2020, ApJ, 888, L3
  * Salas et al. (2019) Salas, J. M., Naoz, S., Morris, M. R., & Stephan, A. P. 2019, MNRAS, 487, 3029
  * Samsing & Ilan (2018) Samsing, J. & Ilan, T. 2018, MNRAS, 476, 1548
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Schwab (2019) Schwab, J. 2019, ApJ, 885, 27
  * Schwab (2021) Schwab, J. 2021, ApJ, 906, 53
  * Shappee & Thompson (2013) Shappee, B. J. & Thompson, T. A. 2013, ApJ, 766, 64
  * Shara & Regev (1986) Shara, M. M. & Regev, O. 1986, ApJ, 306, 543
  * Shen et al. (2018) Shen, K. J., Kasen, D., Miles, B. J., & Townsley, D. M. 2018, ApJ, 854, 52
  * Sills et al. (2005) Sills, A., Adams, T., & Davies, M. B. 2005, MNRAS, 358, 716
  * Sills et al. (2002) Sills, A., Adams, T., Davies, M. B., & Bate, M. R. 2002, MNRAS, 332, 49
  * Sills et al. (2001) Sills, A., Faber, J. A., Lombardi, James C., J., Rasio, F. A., & Warren, A. R. 2001, ApJ, 548, 323
  * Sills et al. (1997) Sills, A., Lombardi, James C., J., Bailyn, C. D., et al. 1997, ApJ, 487, 290
  * Silsbee & Tremaine (2017) Silsbee, K. & Tremaine, S. 2017, ApJ, 836, 39
  * Smeyers & Willems (2001) Smeyers, P. & Willems, B. 2001, A&A, 373, 173
  * Smith (2011) Smith, N. 2011, MNRAS, 415, 2020
  * Smith et al. (2018) Smith, N., Andrews, J. E., Rest, A., et al. 2018, MNRAS, 480, 1466