# Dynamical Fates of S-Type Planetary Systems in Embedded Cluster Environments

Elizabeth A. Ellithorpe  HL Dodge Department of Physics & Astronomy,
University of Oklahoma, Norman, OK 73019, USA  Nathan A. Kaib  HL Dodge
Department of Physics & Astronomy, University of Oklahoma, Norman, OK 73019,
USA  corresponding author<EMAIL_ADDRESS>
###### Abstract

The majority of binary star systems that host exoplanets will spend the first
portion of their lives within a star-forming cluster that may drive dynamical
evolution of the binary-planet system. We perform numerical simulations of
S-type planets, with masses and orbital architecture analogous to the solar
system’s 4 gas giants, orbiting within the influence of a $0.5\>M_{\odot}$
binary companion. The binary-planet system is integrated simultaneously with
an embedded stellar cluster environment. $\sim$10% of our planetary systems
are destabilized when perturbations from our cluster environment drive the
binary periastron toward the planets. This destabilization occurs despite all
of our systems being initialized with binary orbits that would allow stable
planets in the absence of the cluster. The planet-planet scattering triggered
in our systems typically results in the loss of lower mass planets and the
excitement of the eccentricities of surviving higher mass planets. Many of our
planetary systems that go unstable also lose their binary companions prior to
cluster dispersal and can therefore masquerade as hosts of eccentric
exoplanets that have spent their entire histories as isolated stars. The
cluster-driven binary orbital evolution in our simulations can also generate
planetary systems with misaligned spin-orbit angles. This is typically done as
the planetary system precesses as a rigid disk under the influence of an
inclined binary, and those systems with the highest spin-orbit angles should
often retain their binary companion and possess multiple surviving planets.

Keywords: stars: binaries: general, planets and satellites: dynamical
evolution and stability, galaxies: star clusters: general.

##  1 Background & Motivation

The dynamics of planetary systems serve as signposts of their formation and
subsequent evolution. Thus the study of dynamical modeling of exoplanets can
be thought of as analogous to archaeology: using both observations and
simulations of existing stable planets, we can walk back in time to discover
which physical processes govern the behavior of extrasolar worlds. Dynamics
offers a strong union of computational and observational astronomy, where an
influx of exoplanet observational data has dovetailed nicely with a surge in
advancements in numerical algorithms over the last two decades.
It has been long thought that most stars spend some of their early lives in a
cluster environment, a large scale body of gas that acts as a stellar nursery
that condenses into pre-stellar cores in the gas and dust (Lada & Lada, 2003).
Recent work (Sadavoy & Stahler, 2017) has also bolstered the theory (Kroupa,
1995) that the majority of stars are born in a stellar multiple, with their
dense pre-stellar cores mutually bound together in a dusty envelope. These
multiple bound cores, likely formed via fragmentation (Bonnell & Bate, 1994;
Bate et al., 2003; Turk et al., 2009), either stay together in a binary as
interactions with their gaseous disk decrease the separation between the pair
(Bate, 2000; Bate et al., 2003), or are pulled apart by a cluster environment
(Kroupa, 2001). While clusters can unbind so-called ‘soft binaries’ (Kroupa,
2001), dynamical interactions with a cluster are not sufficient to explain the
observed period distribution of binary stars. Thus, their nascent orbital
architecture results from the fragmentation process and interactions with a
circumstellar disk early on in their lives (Kroupa & Burkert, 2001), and
should be considered when we discuss the formation and the dynamics of planets
in a holistic manner.
While planets in a binary system have been previously studied in depth with
N-body integration, the binary that interacts with the planets is for the most
part unperturbed by outside forces (e.g. Holman & Wiegert, 1999; Haghighipour
& Raymond, 2007), or, if the binary is modeled in a cluster, planets are not
included (e.g. Kroupa & Burkert, 2001). Similarly, planetary systems as they
might interact with a cluster environment have been modeled, but often with
mock fly-bys of a star rather than a fully evolved cluster (e.g. Reche et al.,
2009; Cattolico & Capuzzo-Dolcetta, 2020). A union of these two simultaneous
phenomena, a binary’s effect on planets and a cluster’s effect on a binary, is
necessitated by the ubiquitous nature of these interactions in young systems.
With N-body simulations produced by hybrid symplectic integration scheme, we
offer a new perspective at this complex problem.
The presence of a binary companion has been considered in studies of planetary
formation since the 1960’s. Early treatments of planetary stability in a
binary system worked within the framework of the restricted 3-body problem
(Huang, 1960) and numerical estimates based on analytic models (Heppenheimer,
1978). Early work established that both the formation and the stability of
planets in a binary system were possible. Later on, numerical integration to
evolve these systems began in earnest, owing largely to developments in
symplectic integration methods (e.g. Wisdom & Holman, 1991) that allowed for
efficient integrations of entire planetary systems in the presence of a
binary. The main question that arises in such a system is in what cases does a
binary companion preclude the existence of stable planets. An early work that
used the Wisdom-Holman scheme is the exploration of the Alpha Centauri system
in Wiegert & Holman (1997), in which the authors used N-body integrations in
the Wisdom & Holman (1991) scheme to assess the stability of test particles in
and around the Alpha Centauri triplet.
The seminal works in the field of planetary stability in a binary system are
Dvorak (1986) and Holman & Wiegert (1999). Dvorak (1986) looked at integrating
systems of test particles in the restricted 3-body problem with the Lie series
method (Delva, 1985) and established an upper and lower critical orbit that
began mapping a region of stability for a planet based on the eccentricity of
the binary system. Dvorak (1986) looked at various eccentricities of the
binary companion, but in all cases the mass of the binary and the central star
were the same. Holman & Wiegert (1999) further built upon this stability
analysis via symplectically integrating binary systems with both P-type (where
the planet orbits around a binary star pair) and S-type planets (where planets
are bound to a single star with a binary companion exterior to their orbit).
They produced an analytic form of a planet’s critical orbit based on both the
eccentricity of the binary and the stellar mass ratio. Also drawing on Wisdom
& Holman (1991) is the work described in Innanen et al. (1997), in which the
solar system is integrated in the presence of a binary. The authors
exemplified the Kozai effect (Kozai, 1962) in which an inclined binary
exchanges angular momentum with the planets, and also showed that the planets
can evolve in concert as a rigid disk. This work revealed that the Kozai
effect is a pathway that can lead to increased eccentricity of planets.
Another analysis of the Kozai effect on planetary stability that takes a
different, secular approach is Fabrycky & Tremaine (2007), in which the
authors integrated the motion of planets in the presence of a binary to show
that tidal friction can circularize an orbit made eccentric by an inclined
binary companion.
The stability analysis of Holman & Wiegert (1999) remains of interest to this
day, with recent work such as Lam & Kipping (2018) further refining their
original critical orbit with a neural network and integrations performed in
the publicly available REBOUND integrator (Rein & Liu, 2012). Quarles &
Lissauer (2018) also further examined the stability of planets in the presence
of a stellar binary, focusing in particular on the planetary architectures
appropriate for orbital stability in the $\alpha$ Centauri system. Analysis of
planets in binaries continued with work in Chambers et al. (2002) in which the
authors used a mixed variable symplectic method to analyze planetary accretion
both in a system where the planetary bodies orbit a primary and are perturbed
by a distant binary and one in which they orbit around the binary system in
its entirety. While Chambers et al. (2002) utilized an integrator that allows
for the presence of a binary, further work by Beust (2003) introduced a
symplectic scheme based on the Wisdom & Holman (1991) method, called
Hierarchical Jacobian Symplectic (HJS), that allows for the integration of a
multiple system of any size (i.e., numbering more massive bodies than just a
binary) as long as there is a retained hierarchy among the masses. The
presence of a distant companion has also been considered in numerical
integrations of protoplanetary (debris) disks. Reche et al. (2009) used the
HJS scheme to integrate the HD141596 triple star system and debris disk in the
presence of stellar flybys. Beust et al. (2014) used symplectic integration to
model the Fomalhaut triplet in addition to its dust belt.

On the larger scale, N-body simulations of stellar clusters have been used to
study how stellar binaries interact with a cluster environment. Bate et al.
(2003) used high resolution simulations of a collapsing gas cloud to study the
fragmentation of gas into dense cores, and the subsequent evolution of these
young stars into binaries. They found that very tight binaries are formed by
the hardening of wider binaries via dynamical interactions. Adams et al.
(2006) and Proszkow & Adams (2009) used N-body simulations of moderately sized
stellar clusters (100-10,000 members) or embedded clusters to study the impact
of a cluster environment on planet formation, particularly in how
protoplanetary disks may be affected by stellar raditation. They found that
disruption of model solar systems in an evolved cluster environment (Adams et
al., 2006) should be relatively rare due to the paucity of encounters between
a ‘solar system’ in question and a passing star, but in younger, denser
clusters the photoevaporation of disks by a cluster environment can be
appreciable (Proszkow & Adams, 2009).
Also addressing the evolution of a stellar cluster is work by Parker & Goodwin
(2009) and Parker et al. (2009) in which N-body simulations are used to
analyze the prevalence of planets susceptible to the Kozai effect from a
binary companion, and the stability of binaries in a dense cluster environment
respectively. Parker & Goodwin (2009) found that around $20\%$ of all
exoplanets should at one point in their lives be in the presence of a binary
companion that has been sufficiently inclined by its cluster environment such
that Kozai cycles can occur. This is an intriguing finding that bolsters the
idea that a binary companion perturbed by the cluster can have an appreciable
affect on the planets, with the caveat that the authors only examined the
evolution of binary orbits, as their simulations lacked planetary bodies.
Additionally, they particularly focused on very dense clusters similar to
Orion, with a half-mass radius of only $0.1$ pc. This is not typical for
embedded clusters, which have typical half mass radii of $\sim$0.8 pc (Lada &
Lada, 2003), and therefore would have less frequent interactions between
cluster stars and a particular binary. Nevertheless, they showed that a dense
cluster environment can significantly alter the architecture of a stellar
binary which in turn can affect a protoplanetary disk or mature planet system.
Parker et al. (2009) focused on similarly dense cluster environments, but
instead explored the longevity of moderately wide ($\sim$$10^{3}\text{au}$) to
ultrawide ($>10^{4}\text{au}$) binaries. They found that cluster environments
strip away all ultra-wide binary companions, and that the denser clusters,
with half mass radii of 0.1–0.2 pc do not retain any binaries with separations
$>10^{3}\text{au}$. The less dense clusters, with half mass radii 0.4–0.8 pc,
do retain some of these moderately wide binaries. The authors noted that as
ultrawide binaries are often stripped in only a few cluster crossing times,
these very separated binaries may form in isolation An alternate channel is
that ultrawide binaries are formed during cluster dissolution (Kouwenhoven et
al., 2010).
Hao et al. (2013) takes a monte carlo approach to simulating planetary systems
in an open cluster to explore planet-planet scattering. The authors modeled
stellar flybys and found that multi-planet systems are more sensitive to an
open cluster environment than single planet systems. In the realm of exploring
planetary orbits in clusters via a model stellar fly-by is work by Malmberg et
al. (2011) and Breslau & Pfalzner (2019). Malmberg et al. (2011) showed that
fly-bys increase the chance of planetary ejection, while Breslau & Pfalzner
(2019) showed that a cluster star can create a retrograde planetary orbit.

In summary, while there is an extensive body of work on the integration of
planets in a binary system and the evolution of planets in a cluster
environment, previous work largely assumes that the binary is in isolation and
does not evolve due to external forces; work concerning cluster environments
either model interactions as a flyby or do not include a simultaneous binary
system. The work we present here is novel in its approach in that it allows
the binary companion to be altered by a cluster environment, and that we fully
integrate the passages from and evolution of the cluster rather than taking a
flyby approach.

We are largely concerned with the overall fate of binary systems and their
planets during these interactions with the cluster. We find that encounters
between the binary system and clusters stars can destabilize planets and occur
at a non-negligible rate. Moreover, our work has also revealed some intriguing
changes in orbital architecture of planetary systems that have a close
encounter with an excited binary companion. The remainder of this paper is
laid out into the following sections: Section 2 describes how we approach the
numerical integration of these systems as well as how we set up the explored
parameter space. Section 3 details our findings from these simulations
focusing on the mechanism and rate of planetary destabilizations. We then
examine the aftermath of systems that have undergone a planetary instability
in Section 4, highlighting potential tracers of instability and discussing the
broader implications of our results. Finally, we summarize our conclusions in
Section 5.

##  2 Methods

The numerical integrator we use for this work is built on the foundation of
the original Mercury package (Chambers, 1999; Chambers et al., 2002) but with
changes to allow the inclusion of multiple bound stellar mass bodies (i.e., a
binary or triplet star system) and unbound stellar mass bodies (i.e. stars in
a cluster environment). The details of the changes to the integrator are
explained in Kaib et al. (2018). We still use the democratic heliocentric
coordinates for the planetary system, but the binary is defined relative to
the center of mass of the planetary system, and the cluster stars are defined
simply with their inertial coordinates. By letting the different bodies be
treated with different coordinate systems, this version of Mercury retains the
advantages of a pure symplectic integrator while still being able to evolve
the unbound cluster stars. That is, there is not a build up of energy error,
the pseudo-Keplerian nature of the planetary system integration is retained,
and simultaneously the cluster stars and their close approaches can be
integrated efficiently in a leapfrog-like scheme.

###  2.1 Simulation Algorithm

In this integration scheme, we define 3 distinct types of bodies: planets
bound to the primary, the binary companion, and cluster stars. The cluster
star positions are treated with their inertial coordinates, measured only with
respect to the cluster’s center of mass. The binary’s position is integrated
with respect to the center of mass of the primary star system (Chambers et
al., 2002). The planets are integrated in democratic heliocentric coordinates
with respect to the primary star (Duncan et al., 1998). The choice of these
coordinates means that the cluster stars are integrated with a purely
symplectic $T+V$ leapfrog scheme. This conserves energy and angular momentum,
and is quite accurate as the timescale of cluster stars orbiting in the
cluster is much larger than the simulation time step set by the orbits of the
planetary bodies. The binary and the planets are also integrated
symplectically, although with a mixed variable symplectic Wisdom-Holman scheme
(Wisdom & Holman, 1991). Finally, when interacting bodies near each other and
non-Keplerian terms in the Hamiltonian become comparatively large, the bodies
are integrated directly with a Bulirsch-Stoer scheme (Chambers, 1999). For the
planets, this occurs when passing within a Hill radius of one another. For
cluster stars with much larger spheres of influence, the changeover to a
direct integration occurs when two stars are close enough to one another that
their free fall time is less than $500\times$ the simulation time-step of
$100$ days. For a 1 $M_{\odot}$ star, this corresponds to a close encounter
distance of $\sim$130 au.

Thus, this scheme largely preserves the angular momentum and energy
conservation of a symplectic scheme while allowing for close encounters.
Equation 1 shows how we define the position of each type of object, where
$\vec{X}_{A}$ represents the position of the primary star, $\vec{X}_{B}$ the
position of the binary, $\vec{X}_{i}$ for $1\leq i\leq N_{P}$ is the position
of a planet, and $\vec{X}_{i}$ for $N_{P}<i\leq N_{P}+N_{S}$ is the position
of a cluster star. The lowercase $\vec{x}$ represents each type of body’s
inertial coordinates relative to the origin/center of mass of the cluster.

$\begin{split}\vec{X}_{A}=\frac{m_{a}\vec{x}_{A}+m_{B}\vec{x}_{B}+\sum_{j=1}^{N_{P}}m_{j}\vec{x}_{j}}{m_{A}+m_{B}+\sum_{j=1}^{N_{P}}m_{j}}\\\
\vec{X}_{i}=\vec{x}_{i}-\vec{x}_{A}\text{ for }1\leq i\leq N_{P}\\\
\vec{X}_{B}=\vec{x}_{B}-\frac{m_{A}\vec{x}_{A}+\sum_{j=1}^{N_{P}}m_{j}\vec{x}_{j}}{m_{A}+\sum_{j=1}^{N_{P}}m_{j}}\\\
\vec{X}_{i}=\vec{x}_{i}\text{ for }N_{P}<i\leq N_{P}+N_{S}\end{split}$ (1)

As such, we define some useful position vectors for changing between reference
frames. $\vec{s}$ marks the barycenter of the planets relative to the primary
host star and $\vec{\Delta}$ marks the inertial position of the primary star
as defined with the new coordinates ($\vec{X}$’s).

$\begin{split}\vec{s}=\frac{\sum_{i-1}^{N_{P}}m_{i}\vec{X}_{i}}{m_{A}+\sum_{i-1}^{N_{P}}m_{i}}\\\
\vec{\Delta}=\vec{X}_{A}-\frac{\sum_{i=1}^{N_{P}}m_{i}\vec{X}_{i}+m_{B}(\vec{X}_{B}+\vec{s})}{m_{A}+m_{B}+\sum_{i=1}^{N_{P}}m_{i}}\end{split}$
(2)

The corresponding conjugate momenta for the bodies is as follows, again where
the lowercase terms are for the inertial case and the uppercase terms are the
momenta for the new coordinates:

$\begin{split}\vec{P}_{A}=\vec{p}_{A}+\vec{p}_{B}+\sum_{j=1}^{N_{P}}\vec{p}_{j}\\\
\vec{P}_{i}=\vec{p}_{i}-m_{i}\frac{\vec{p}_{A}+\sum_{j=1}^{N_{P}}\vec{p}_{j}}{m_{A}+\sum_{j=1}^{N_{P}}m_{j}}\text{
for }1\leq i\leq N_{P}\\\
\vec{P}_{B}=\vec{p}_{B}-m_{B}\frac{\vec{p}_{A}+\vec{p}_{B}+\sum_{j=1}^{N_{P}}\vec{p}_{j}}{m_{A}+m_{B}+\sum_{j=1}^{N_{P}}m_{j}}\\\
\vec{P}_{i}=\vec{p}_{i}\text{ for }N_{P}<i\leq N_{P}+N_{S}\end{split}$ (3)

In this symplectic scheme, we must consider how Hamilton’s equations evolve
for all members of the system. We consider mutual gravitation between the
primary star, binary star, planets and cluster stars. We also consider the
gravitational tide from the gas of the Plummer potential (Plummer, 1911) that
we assume to be binding the cluster stars together. We include the derivation
for the Plummer potential here. We treat the force from a Plummer distribution
as a constant background potential whose origin is the center of mass of the
cluster stars. As the effect of the Plummer potential is dependent on the
inertial position of each massive body within the cluster, the calculation in
our chosen (non-inertial) coordinates is fairly complicated. The expression
for the force from the sphere of gas is different for each type of body,
namely the planets, the binary, and the cluster stars. We begin with the
expression for the mass distribution of a Plummer sphere in inertial
coordinates, where $\vec{x}$ is defined relative to the cluster’s center-of-
mass.

$\rho(\vec{x})=\frac{3M_{0}}{4\pi
a_{0}}\times\biggl{(}1+\frac{|\vec{x}|^{2}}{a_{0}^{2}}\biggr{)}^{-5/2}$ (4)

This mass distribution, where $M_{0}$ is the total mass in gas of the Plummer
sphere and $a_{0}$ is a scale parameter that sets the size of the flat core
region results in a gravitational potential of:

$\Phi(\vec{x})=-\frac{GM_{0}}{(|\vec{x}|^{2}+a_{0}^{2})^{1/2}}$ (5)

The potential energy of our system can then be written as:

$\begin{split}V_{Plum}=-\sum_{i=N_{p}+1}^{N_{p}+N_{s}}\frac{GM_{0}m_{i}}{(|{\vec{X}_{i}}|^{2}+a_{0}^{2})^{1/2}}-\sum_{i=1}^{N_{p}}\frac{GM_{0}m_{i}}{(|(\vec{\Delta}+\vec{X_{i}})|^{2}+a_{0}^{2})^{1/2}}\\\
-\frac{GM_{0}m_{B}}{(|(\vec{\Delta}+\vec{s}+\vec{X}_{B})|^{2}+a_{0}^{2})^{1/2}}-\frac{GM_{0}m_{A}}{(|\vec{\Delta}|^{2}+a_{0}^{2})^{1/2}}\end{split}$
(6)

For clarification, the inertial position with which the strength of the
Plummer potential is measured in simulation coordinates is $\vec{X}_{i}$ for
the cluster stars, $\vec{\Delta}+\vec{X}_{i}$ for the planets,
$\vec{\Delta}+\vec{s}+\vec{X}_{B}$ for the binary, and $\vec{\Delta}$ for the
primary. We can then calculate the acceleration due to the Plummer sphere on
each body (indexed with $i$) for each cartesian coordinate (indexed with $u$)
via:

$m_{i}\frac{dv_{i,u}}{dt}=-\frac{\partial V_{plum}}{\partial x_{i,u}}$ (7)

####  2.1.1 Cluster Stars

The positions of the cluster stars in the integration coordinates are exactly
their inertial coordinates in this integration scheme (i.e., they are measured
relative to the origin). This makes finding the acceleration due to the
Plummer tide much more straightforward than for the planets and the binary
star pair:

$\vec{X}_{i}=\vec{x}_{i}\text{ for }i>N_{p}$ (8)

$\frac{dv_{i,u}}{dt}=-\frac{GM_{0}X_{i,u}}{(|\vec{X_{i}}|^{2}+a_{0}^{2})^{3/2}}$
(9)

We can see that the acceleration due to the Plummer sphere on the cluster
stars is simply a direct term based on the mass of the gas enclosed by the
cluster star’s current position.

####  2.1.2 Planets

Recall that $\vec{\Delta}$ is a function of the position of the planets.

$\frac{\partial\vec{\Delta}}{\partial
X_{i}}=-m_{i}\frac{1}{m_{A}+\sum_{k=1}^{N_{P}}m_{k}}$ (10)

This dependence on the planets’ positions leads to a more complicated
expression for the Plummer potential on the planetary bodies.

$\begin{split}\frac{dv_{i,u}}{dt}=\frac{GM_{0}}{m_{A}+\sum_{j=1}^{N_{P}}m_{j}}\biggl{[}\sum_{k=1}^{N_{P}}\frac{m_{k}(\Delta_{u}+X_{k,u})}{(|\vec{\Delta}+\vec{X_{k}}|^{2}+a_{0}^{2})^{3/2}}+\frac{m_{A}\Delta_{u}}{(|\vec{\Delta}|^{2}+a_{0}^{2})^{3/2}}\biggr{]}\\\
-\frac{GM_{0}(\Delta_{u}+X_{i,u})}{(|\vec{\Delta}+\vec{X_{i}}|^{2}+a_{0}^{2})^{3/2}}\end{split}$
(11)

####  2.1.3 Binary Companion

The vector $\vec{\Delta}$ is also a function of the position of the binary:

$\frac{\partial\vec{\Delta}}{\partial
X_{B}}=\frac{-m_{B}}{m_{A}+m_{B}+\sum_{i=1}^{N_{P}}m_{i}}$ (12)

Thus, as with the planets, the expression for the acceleration of the binary
star due to the Plummer potential is more complicated than what we showed
previously for the cluster stars:

$\begin{split}\frac{dv_{B,u}}{dt}=\frac{GM_{0}}{m_{A}+m_{B}+\sum_{j=1}^{N_{P}}m_{j}}\biggl{[}\sum_{k=1}^{N_{P}}\frac{m_{k}(\Delta_{u}+X_{k,u})}{(|\vec{\Delta}+\vec{X}_{k})|^{2}+a_{0}^{2})^{3/2}}\\\
+\frac{X_{B,u}+s_{u}+\Delta_{u}}{(|\vec{\Delta}+\vec{X}_{B}+\vec{s}|^{2}+a_{0}^{2})^{3/2}}+\frac{m_{A}\Delta_{u}}{(|\vec{\Delta}|^{2}+a_{0}^{2})^{3/2}}\biggr{]}-\frac{GM_{0}(X_{B,u}+s_{u}+\Delta_{u})}{(|\vec{\Delta}+\vec{X}_{B}+\vec{s}|^{2}+a_{0}^{2})^{3/2}}\end{split}$
(13)

This modified version of the Mercury simulation package is publicly available
for download at https://github.com/nathankaib/Mercury_StarCluster.

###  2.2 Cluster Environments

First, we must decide how to build our stellar clusters. Our clusters consist
of many discrete stellar mass bodies embedded within a more massive background
Plummer potential to mimic the gravitational effects of a dominant gaseous
component to the clusters. To design the parameters of the stellar cluster, we
turn to the embedded cluster catalog compiled by Lada & Lada (2003), which
Adams et al. (2006) find is mostly comprised of clusters with 100–1000 stars.

While 1000-star clusters begin to have prohibitive runtimes, we employ small,
medium, and large cluster environments with 122, 221, and 509 stars,
respectively. These have total stellar masses of 80, 160, and 341 M⊙. In
embedded clusters, the stellar mass only represents a small fraction of the
total cluster mass, as Lada & Lada (2003) estimate embedded cluster star
formation efficiencies between $\sim$10–25%. To bracket the low end of this
range, our small and medium clusters have Plummer masses of 1000 and 2000 M⊙,
yielding star formation efficiencies of 7.4%. Meanwhile, our large cluster
employs a Plummer mass of 1364 M⊙, yielding a higher star formation efficiency
of 20%. To highlight the role of the non-stellar potential, we also include a
version of our medium cluster with all stellar masses reduced by an order of
magnitude, implying a star formation efficiency of just 0.74%.

To set the Plummer radii of our potentials (and stellar distributions), we
again consult Lada & Lada (2003). Using this catalog, Adams et al. (2006)
estimate that 122-, 221-, and 509-star clusters (our small, medium, and large
clusters) have typical radii of 0.63–1.1, 0.85–1.5, and 1.3–2.3 pc,
respectively. For our small and medium cluster environments, we set our
Plummer radius to 0.62 pc, implying a half-mass radius of $\sim$0.8 pc. For
the large cluster environment, we set the Plummer radius to 1.38 pc, yielding
a half-mass radius of 1.8 pc.

Our cluster environment parameters are summarized in Table 1. The goal of this
work is to probe a significant swath of cluster parameter space relevant to
many planetary systems and not necessarily to the solar system. However, it
may be instructive to compare our parameters with constraints on the solar
birth cluster, since our work assumes planetary system architectures like our
own and it has been speculated that the early Sun may have possessed a binary
companion (Siraj & Loeb, 2020). Primarily relying on cluster perturbations to
generate the orbit of Sedna while preserving the stability of the giant
planets, Adams (2010) concludes that the Sun’s birth cluster contained
1500–7100 stars, or roughly an order of magnitude more than our clusters.
However, this assumes a cluster residence time of $\gtrsim$10 Myrs. Using the
Kuiper belt inclination induced by flybys of cluster stars, Batygin et al.
(2020) find that the product of the solar birth cluster residency and stellar
density must be below $\sim$$2\times 10^{4}$ Myr/pc3. The values for our
clusters are 1–2 orders of magnitude less than this upper limit, so it again
appears that much of the allowable parameter space for the solar birth cluster
is more populous or denser than our modeled environments.

Cluster | Plummer | Stellar | Star Formation | Plummer | N∗ | $\bar{n_{*}}$
---|---|---|---|---|---|---
| Mass | Mass | Efficiency | Radius |  |
| (M⊙) | (M⊙) | (%) | (pc) |  | (pc-3)
Small | 1000 | 80 | 7.4 | 0.62 | 122 | 47.4
Medium | 2000 | 160 | 7.4 | 0.62 | 221 | 86.0
Large | 1364 | 341 | 20 | 1.38 | 541 | 22.2
NonStellar | 2000 | 16 | 0.74 | 0.62 | 221 | 86.0
Table 1: Parameter summary of simulated cluster environments. From left to
right, columns are: environment name, mass of Plummer potential, mass of
stellar component, implied star formation efficiency, number of stars, mean
volume density of stars inside half-mass radius.

###  2.3 Simulation Setup and Progression

With cluster parameters specified, our cluster assembly and integration
prescription is largely based on the methods outlined in Levison et al.
(2010), with some minor differences, as noted below. The masses of the stars
are sampled from the Chabrier (2003) piecewise IMF, where we set an upper mass
limit of $150\>M_{\odot}$ (Weidner & Kroupa, 2004). Then, the virial speeds of
each star are calculated based on their position and we assign them a random
initial velocity vector at $8\%$ of the virial speed (Levison et al., 2010).
Once the cluster bodies are positioned following the cluster Plummer profile,
we generate different simulations, where for each simulation a random cluster
star is converted into the role of the primary, such that the starting
position of the planetary/binary system within the cluster is randomized.

A set of 4 gas giants, analogous to the orbits and masses of Jupiter, Saturn,
Uranus, and Neptune, are set around a 1 M⊙ primary star, and a coplanar 0.5 M⊙
binary companion is added as well. We choose to fix our binary mass at this
fiducial subsolar value to limit the number of parameters being simultaneously
varied. Higher companion masses will likely enhance the destabilizing power of
the binary companion, while lower ones should diminish it.

Before integration within a cluster environment, all of our planet+binary
systems are first verified via numerical integration to be inherently stable
for $10$ Myr in isolation. If they are not stable, the initial conditions are
discarded. If they do prove to be stable, they are then integrated for 10 Myrs
within the cluster environment in concert with all of the other cluster stars.
After 5 Myrs of this 10-Myr integration, the Plummer potential is
instantaneously dispersed, leaving only the stars and planets. This
corresponds to the physical dispersal of the cluster gas, and past modeling
works have employed similar dispersal timescales (Adams et al., 2006; Proszkow
& Adams, 2009; Levison et al., 2010; Brasser et al., 2012).

###  2.4 Binary Orbits

Name | Cluster | Binary | Binary | Number of
---|---|---|---|---
| Environment | Eccentricities | Semimajor Axes | Simulations
|  |  | (au) |
Small_Fixed | Small | $e=0.5$ | 300–800 | 600
Medium_Fixed | Medium | $e=0.5$ | 300–800 | 600
Large_Fixed | Large | $e=0.5$ | 300–800 | 600
Medium_Cold | Medium | $f(e)\propto\sqrt{e}$ | 300–800 | 600
Medium_Uniform | Medium | $f(e)\propto e$ | 300–800 | 600
Medium_Thermal | Medium | $f(e)\propto e^{2}$ | 300–800 | 600
Medium_Wider | Medium | $f(e)\propto e^{2}$ | 900–1300 | 500
NonStellar_Fixed | NonStellar | $e=0.5$ | 300–800 | 600
NonStellar_Wider | NonStellar | $f(e)\propto e^{2}$ | 900–1300 | 500
Table 2: Summary of simulation sets. From left to right, columns are:
simulation set name, cluster environment of simulation set (see Table 1),
binary eccentricity distribution, and range of binary semimajor axes (spanned
in 100 au intervals), and total number of simulations.

The influence of the binary on our planetary systems and its susceptibility to
perturbation from our cluster environments of course strongly depends on its
orbit. Our choices for binary orbits vary across our different sets of
simulations. To initially limit the number of parameters being varied within a
simulation set, we first perform three sets of 600 simulations in which the
binary eccentricity is fixed at 0.5. Meanwhile, in each simulation set, the
binary semimajor axes are varied in increments of 100 au between 300 and 800
au, with 100 different simulations (different starting cluster positions)
performed at each binary semimajor axis. (For semimajor axes closer than $300$
au the planets are typically immediately unstable, and some binaries beyond
$800$ au begin to dissociate before cluster gas dispersal.) These fixed
eccentricity runs allow us to better isolate the importance of binary
semimajor axis and cluster environment.

An additional three sets of simulations explicitly study the importance of the
binary’s initial eccentricity within the medium cluster environment. The three
different sets employ different distributions of binary eccentricity: uniform
in $e$, uniform in $e^{2}$ (a thermalized distribution; Ambartsumian, 1937),
and uniform in $e^{1/2}$. For each individual simulation, we again first
verify that the choice of binary eccentricity between $0$ and $1$ from the
requisite distribution resulted in a system of stable planets for $10$ Myr
without the presence of the cluster environment. This gave us an approximate
upper bound of $e=0.8$ for the widest 800 au binaries and $e=0.6$ for the
tightest 300 au binaries. Especially in the case of simulations where the
planets are on the edge of stability, we operate under the assumption that
planets can form in the presence of such an eccentric binary companion. This
can be physically explained either that we start the $10$ Myr simulation at a
snapshot where the cluster has perturbed the binary to its high starting
eccentricity and we are examining how long the planets remain stable after
further perturbations, or that when the planets are forming in a gas disk, the
eccentricity perturbations from cluster stars are damped out by the disk
material. Either way, we are embedding a system that is stable in isolation,
and examining how the cluster affects the binary from that point.

We also run another follow-on set of simulations that explore slightly larger
binary semimajor axes. As mentioned earlier, the binary dissociation times
decrease with binary semimajor axes. To destabilize the planets from an
$e=0.5$ starting eccentricity, the binary periastron has to decrease
dramatically in a widely separated binary. However, our thermalized runs begin
many binaries on higher eccentricities nearer to the point of planetary
instability. To study whether more weakly bound binaries could also trigger a
planetary system instability, this set of simulations spans binary semimajor
axes between 900 and 1300 au in semimajor axis increments of 100 au. These
binaries are also immersed in the medium cluster environment described above.

Finally, we perform two more sets of simulations in the reduced stellar mass
version of our medium cluster environment. These are meant to explore the
importance of the Plummer gas potential to our systems’ dynamics. The first
set employs binaries with a fixed eccentricity of 0.5 and semimajor axes from
300–800 au, and the second set employs a thermal binary eccentricity
distribution with semimajor axes from 900–1300 au. Our various sets of
simulations are summarized in Table 2.

##  3 Results

Name | Planetary | Binary
---|---|---
| Disruption Rate | Dissociation Rate
| (%) | (%)
Small_Fixed | 5.5 | 8.2
Medium_Fixed | 9.0 | 13.0
Large_Fixed | 7.3 | 8.8
Medium_Cold | 6.7 | 11.5
Medium_Uniform | 10.8 | 12.7
Medium_Thermal | 14.5 | 13.2
Medium_Wider | 14.4 | 22.8
NonStellar_Fixed | 2.0 | 0.3
NonStellar_Wider | 8.0 | 0.8
Table 3: Summary of instability rates in our simulation sets. From left to
right, columns are: simulation set name, percentage of destabilized planetary
systems, and percentage of unbound binary companions. Figure 1: A: The
fraction of planetary systems that experienced one or more planetary ejections
as a function of initial binary semimajor axis with fixed initial binary
eccentricity of 0.5. Different cluster mass simulations are co-added and error
bars are 1$\sigma$ Poisson uncertainties. B: The cumulative distribution of
final binary periastra for binaries beginning with $e=0.5$. Distributions are
obtained via co-adding all three different cluster mass simulations.

The rates that planetary systems are destabilized are shown in Table 3. We see
no strong dependence on the cluster masses we employ. Amongst our fixed
eccentricity binaries, we find that 5.5%, 9.0%, and 7.3% of the systems in our
small, medium, and large clusters experience planetary system instabilities
(lose at least one planet via collision or ejection). Surprisingly, we also do
not observe a strong dependence on initial binary semimajor axis. In Figure
1A, we plot the fraction of planetary systems that undergo an instability as a
function of binary semimajor axis for our fixed binary eccentricity runs. As
can be seen, there is no clear trend across binary semimajor axes of 300–800
au, with instability rates hovering between $\sim$5–9%. We can see why this is
the case in Figure 1B, where we examine the final periastron distributions for
$a=300$ au binaries and $a=800$ au binaries. Here we see that the two
distributions are radically different, retaining signatures of their identical
initial eccentricities but differing semimajor axes. However, at the low
periastron tail, the two samples have nearly identical fractions of systems
with periastra below $\sim$100 au. Coincidentally, periastron passages within
100 au are what are typically capable of destabilizing our planetary systems,
where the instability probability quickly falls off for more distant
periastron passages. Thus, the greater periastron evolution experienced by our
$a=800$ au binaries roughly offsets the effects of the lower initial periastra
of our $a=300$ au binaries.

For our variable eccentricity simulations, we see a clearer sign that
instability rates vary with different initial distributions of binary
eccentricity. The fraction of planetary instability is $6.7\%$ for the low
eccentricity distributions, $10.8\%$ for the uniform case, and $14.5\%$ for
the thermalized case, steadily increasing as the initial binary eccentricities
are biased toward higher values. We can see this dependence more obviously in
Figure 2. Here we co-add the results of our Medium_Cold, Medium_Uniform, and
Medium_Thermal runs and then separate them by binary semimajor axis. The
fraction of planetary systems that have undergone instability is then plotted
against the initial binary eccentricity normalized by the maximum binary
eccentricity used for a given binary semimajor axis. (Since we require systems
to be stable in isolation, different binary semimajor axes allow for different
maximum eccentricities.) This plot demonstrates that the probability of
planetary instability increases by a factor of $\sim$3–5 across the range of
eccentricities employed for each binary semimajor axis.

Figure 2: Planetary instability rate as a function of initial binary
eccentricity (normalized by the maximum initial eccentricity employed for a
given binary semimajor axis) in our Medium_Cold, Medium_Uniform, and
Medium_Thermal runs. Different lines and symbols correspond to different
binary semimajor axes, as detailed in the legend. Legend also lists the
maximum initial binary eccentricity for each semimajor axis.

Although we see a clear dependence of planetary instability probability on the
binary’s eccentricity, we again do not see a strong dependence on binary
semimajor axis. Even though some of our wider binaries begin to be dissociated
by the end of the cluster’s life and also need to attain a more extreme
eccentricity to disrupt their planets, this appears to be roughly offset by
the greater orbital variation these binaries incur due to cluster
perturbations. This raises the possibility that significant planetary
instability rates persist to wider binaries. Our Medium_Wider simulations
suggest that this is the case. Extending the binary semimajor axes from 800 au
to 1300 au (using a thermal binary eccentricity distribution), we see no sign
of a significant drop off in planetary instability rates, with 14.4% of
systems undergoing an instability. This is in spite of the fact that our
widest binaries ($a=1300$ au) need to typically reach a much more extreme
eccentricity of $\sim$0.93 to destabilize the planets compared to our tightest
binaries, which only need to typically exceed $e=0.62$.

A uniting factor that we can point to in the case of destabilized systems is
the number of close encounters a system has with cluster stars, as the
perturbing of the binary by the stochastic motion of the cluster can produce
planetary instability. This is what we see in the vast majority of cases:
interactions with passing cluster stars drive the binary companion to a very
eccentric orbit, allowing for close pericenter passages with the primary and
its planets. Kaib et al. (2013) found a similar result in an exploration of
wider binaries in the galactic field wherein even distant binary companions
can be perturbed to eccentric orbits with close pericenter approaches to an
inner planet system. In fact, in our simulations the frequency of close
stellar encounters seems to be the most important factor in whether or not a
system is destabilized. Unstable systems have a clear preference for multiple,
close interactions with cluster stars. This is shown in the histogram in
Figure 3, where we count the number of times cluster stars pass within
$10^{4}$ au of our binary system. A clear takeaway from this figure is that
destabilized systems typically experience far more close encounters with
stellar cluster members than stable systems. The probability of having a large
number of close stellar encounters within the cluster is largely dependent on
the positional history of our systems within the cluster environment. For each
of our simulations, we measure the distance from the primary to the cluster
center-of-mass at each time output. For systems that lose at least one planet,
this distance’s median value is over 4 times smaller compared to those that
remain stable.

Figure 4 illustrates why a large number of encounters with cluster stars
enhances the probability of a planetary instability. The example system shown
here retains its binary companion, but it evolves from an initial semi-major
axis of $a=800$ au to $a\approx 500$ au. At the same time, the binary
periastron evolves from 400 au to under 100 au, leading to strong perturbation
of the planetary orbits. During this inwards binary evolution, Jupiter and
Saturn are pulled to $\sim$$100^{\circ}$ inclination (as the binary’s orbital
inclination also increases to over 130∘), while Uranus and Neptune are ejected
from the system. This particular simulation is a good representative of many
common characteristics we observe in unstable systems: 1) the binary star
moves inward, and its close periastron approaches drive the loss of the
planets. 2) The massive inner planets are the most likely survivors of this
event. 3) The remaining bound planets are on altered orbits following the
binary’s periastron approach, having both higher inclination and eccentricity
on average than their stable counterparts. As frequent close encounters from
cluster stars are the main driver in these orbital changes, the architecture
of planetary system is largely ‘frozen in’ once the cluster becomes unbound
after $5$ Myr (Allen et al., 2007), as the frequency of stellar encounters
drops precipitously. A similar lack of stellar encounters when the stars are
unbound from the gas was seen in Adams et al. (2006).

In $73\%$ of our destabilized systems, at least one planet is left behind,
which is good news for observational signals of these disruption events. This
may allow us to make inferences about orbital changes due to binary
destabilization processes in surviving exoplanets. Figure 4 suggests that
planets surviving these binary-triggered instabilities may reside on
exceptionally eccentric and/or inclined orbits.

Figure 3: A normalized histogram containing results from all 3 suites of
simulation depicting the distribution of number of ‘close’ (within 10k au)
encounters the planetary system has with cluster stars. The systems with
unstable planets are in red, showing a tail of multiple close encounters. The
stable systems are shown in black. Figure 4: The orbital evolution of the 4
giant planets and binary companion. The left panels are the distance of the
object from the primary, where the solid black line is the semi-major axis and
the dashed red lines are the pericenter/apocenter of the body. This particular
system results in a very high predicted spin-orbital angle for the Jupiter and
Saturn mass bodies. Note the multiple ‘jumps’ as the binary companion enters a
close pericenter approach. These jumps are from encounters with cluster stars
in which angular momentum is exchanged. In this simulation both the Uranus and
Neptune mass planets become unbound and the binary companion remains bound.

We can also examine the final fate of the binary companion, and we use our
fixed binary eccentricity runs to do so. For unstable systems with higher
planetary inclinations, this effect is often tied to a similar change in the
binary orbit. Among our unstable systems that retain a binary companion and at
least one planet at the end of the simulation, the inner planet has a median
inclination of 38.4∘. If we just look at the systems with planetary
inclinations higher than this median, we find that the final median binary
inclination is 52.1∘. However, when we look at unstable systems in the bottom
half of planetary inclinations, we find a much lower final median binary
inclination of just 12.1∘.

Of course, not all binary companions remain bound until the end of our
simulations. Out of our 1800 fixed eccentricity simulations, binary disruption
is relatively rare, with 10% of binaries becoming dissociated. However,
amongst systems that underwent an instability, the rate of binary dissociation
is much higher at 56%. (This is due to the fact that exceptional numbers of
close encounters with cluster stars drive the binary evolution required to
trigger a planetary instability.) Thus, for the majority of unstable systems,
the binary that perturbed the planets is eventually lost to the cluster. Many
single-star exoplanet systems should exist in which the planetary orbital
architecture exhibits hallmarks of a violent instability, yet the
destabilizing binary companion responsible for this is no longer associated
with the system.

It is becoming recognized that the exoplanets with the highest known
eccentricities have a higher probability of possessing a stellar companion
(Venner et al., 2021). In fact, $\sim$50% of exoplanets whose eccentricities
exceed 0.8 reside in multi-star systems. Our instabilities triggered by a
quickly vacated binary companion may also explain a substantial fraction of
the $\sim$50% of such eccentric exoplanets that are not observed to currently
reside in a stellar multiple.

##  4 Discussion

###  4.1 Orbital Architecture Changes

The orbital architectures of surviving planets in unstable (lose at least one
planet) versus stable (retain all four planets) have at least two potentially
observable signposts: the eccentricity and spin-orbit angle of the surviving
planets. In the following analysis, we take the final planetary orbital
inclinations as a proxy for their spin-orbit angles. This makes two implicit
assumptions. The first is that the planets are formed on orbits whose angular
momentum is aligned with the stellar spin axis. The second is that tidal
forces do not subsequently realign the spin axes of our stars with the orbital
angular momentum vectors of our planets. Although the sample of known Hot
Jupiters with large spin-orbit angles suggests that tidal realignment can
occur after spin-orbit angle excitation, the timescale for realignment is not
well-constrained, and such timescales would be much longer in our simulated
systems whose planets orbit at multiple au (Winn et al., 2010; Hamer &
Schlaufman, 2022).

The distributions of eccentricity and the predicted spin-orbit angle of the
inner gas giants in unstable systems are statistically distinct from their
stable counterparts. High eccentricity planets, like the ones we see in our
unstable systems, are thought to have a dynamical source. First posited in
Rasio & Ford (1996), planet-planet scattering interactions can lead to
ejections of some planets and a higher eccentricity for planets that remain.
In our systems, often these surviving planets are the two inner gas giants.
The Jupiter and Saturn mass planets thus act as good markers of a past
instability, and our analysis will focus on their final orbits in each
simulation. The distributions of the final eccentricity and inclination of the
inner giants are shown in Figure 5. Each of our batches of simulations tells a
similar story. The median eccentricities of Jupiter and Saturn analogs in our
stable systems are around 0.045, but in systems that have undergone binary-
triggered instabilities this median value grows by 50–100%. In addition, a
small fraction of unstable systems are driven to more exceptional
eccentricities. Depending on the particular simulation batch, 1–3% percent of
our surviving gas giants achieve eccentricities above 0.5.

We see similar trends in the evolution of planetary inclinations and assumed
spin-orbit angle. In stable systems, the Jupiter and Saturn analogs have a
median predicted spin-orbit angle of 2–4∘. However, in instability survivors
this median value increases to 15–30∘, depending on the particular simulation
batch. Often the two inner planets will undergo a coupled evolution as if in a
rigid disk (shown in work by Innanen et al. (1997)), reaching similar
inclinations after the instability. In addition, approximately 5% of our
surviving planets in unstable systems attain spin-orbit angles in excess of
90∘.

Figure 5: A: Cumulative distributions of eccentricities of Jupiter and Saturn
analogs in stable (solid lines) and unstable (dashed lines) systems from our
Small_Fixed, Medium_Fixed, and Large_Fixed simulations. B: Cumulative
distributions of the predicted spin-orbit angles of the planets from Panel A.
C: Cumulative distributions of eccentricities of Jupiter and Saturn analogs in
stable (solid lines) and unstable (dashed lines) systems from our Medium_Cold,
Medium_Uniform, Medium_Thermal, and Medium_Wider simulations. D: Cumulative
distributions of the predicted spin-orbit angles of the planets from Panel C.

Interactions with the binary companion, whether it remains bound to the
primary or not, destabilize the planets via an exchange of orbital angular
momentum. If the planets stay bound, this interaction tends to heat up their
orbits to higher eccentricities and inclinations. However, when we examine our
unstable systems, we find that final planetary inclinations and eccentricities
are not particularly strongly correlated ($r<0.1$). This is because the most
common instability leads to the ejection of one or both ice giants, while the
Jupiter and Saturn analogs survive. The scattering of ice giants does not
necessarily strongly drive the excitation of Jupiter and Saturn analog
eccentricities. On the other hand, the binary evolution necessary to trigger
ice giant ejections also often elicits significant rigid disk precession of
the Jupiter and Saturn analogs, significantly increasing the spin-orbit angle
(Kaib et al., 2011). If we instead limit ourselves to systems that only have 1
surviving planet (nearly always Jupiter after a Saturn ejection), we find a
stronger correlation ($r=0.37$) between final inclination and eccentricity. In
addition, the median planetary eccentricity is much higher at 0.25, and this
sample also contains 2/3 of our planets that finished with eccentricities
above 0.5.

###  4.2 A Significantly Misaligned Population

If we assume that the planets in our systems begin in an aligned configuration
with the stellar spin axis, we can use their orbital inclination as a proxy
for the spin-orbit angle $\Psi$. We must then consider how the signature of
misalignment would be maintained along a hypothetical line of sight. To ensure
that the distribution of planetary spin-orbit angles for stable and unstable
populations remain distinct, we choose a random observing angle for each
system and calculate the projected spin-orbit angle $\lambda$. We then can
utilize the Kolmogorov-Smirnov test to conclude if we can discard the null
hypothesis that the projected spin-orbit angles of our stable and unstable
systems are drawn from the same parent distribution. We repeat this process of
choosing a random viewing angle for each system in the population 100 times,
and in all 100 cases we obtain a p-value from our KS test that is less than
$1\times 10^{-3}$. Therefore, even when projecting our systems along a
hypothetical viewing angle, the distributions remain statistically distinct
between stable and unstable systems. A cumulative distribution of this test
for our fixed binary eccentricity simulations is shown in Figure 6, where the
median value of the projected angle is plotted for the stable and unstable
systems and the shaded regions denote the $5^{th}$ and $95^{th}$ percentiles.
This N-body gravitational interaction pathway to misalignment is significant,
in our simulations, and bolsters the observations that high obliquity
planetary systems likely arise from planet scattering events (Winn et al.,
2010). Roughly 1/3 of our systems that experienced instability attain
projected spin-orbit angles over 30∘, while only $\sim$1/20 of our stable
systems do so.

Figure 6: Cumulative distributions of the average projected spin-orbit angle
after 10 Myrs for stable and unstable systems. Our Small_Fixed, Medium_Fixed,
and Large_Fixed simulations are co-added. The shaded regions denote the
$5^{th}$ and $95^{th}$ percentiles of projected spin-orbit angle for each bin.

While our highest eccentricity systems are very likely to be 1-planet systems,
this is not the case for our highest inclination planets that have survived
instabilities. Of the 16 unstable planetary systems driven to spin-orbit
angles over 90∘, 12 of them contain both Jupiter and Saturn analogs. In
addition, the Jupiter and Saturn analogs in our 12 multi-planet systems all
have final inclinations within $\sim$5∘ of one another (as well as modest
eccentricities of $\lesssim 0.1$), indicating that they have attained their
high inclinations via rigid disk precession (Innanen et al., 1997).

Also in contrast to planetary eccentricities, there is some correlation
between spin-orbit angle and the retention of a binary companion. If we co-add
all of our unstable systems and examine those that lost their binary, the
median predicted spin-orbit angle is 11.1∘. Meanwhile, for those unstable
systems that retain their binary, this median angle is over twice as large at
23.5∘. This suggests that the retention of the binary leads to continued
evolution of planetary inclinations (in particular, rigid disk precession)
after the instability.

Given the potentially long timescales for rigid disk precession (see Figure 4
for instance), this raises the prospect that amongst our unstable systems with
a binary companion, the spin-orbit inclinations may continue to evolve after
our 10-Myr integration times. To test this possibility, we perform follow-on
integrations of our unstable systems that possess a binary for another 100
Myrs. The results of these integrations are shown in Figure 7, where we plot
their distribution of spin-orbit angles. Here we see that systems with binary
companions continue to diverge from those that lost their companions. After
the additional 100 Myrs of integration, the median spin-orbit angle has
increased from 23.5∘ to 33.8∘. In addition, the fraction of systems with
retrograde angles has nearly doubled from 6.5% to 11.1%. (It should be noted
that our distribution of spin-orbit angles after 50 Myrs of integration does
not exhibit any statistical differences from our 100-Myr sample, suggesting
this suite of systems has reached a dynamical equilibrium by the end of its
inegration.)

This small but significant population of retrograde systems possesses distinct
properties. Mostly having achieved their spin-orbit angles via rigid disk
precession, the majority (70%) are multi-planet systems with modest
eccentricities ($\bar{e}\simeq 0.09$). Such a set of observed systems would be
very consistent with a dynamical history where a planetary system is perturbed
by a retained binary companion whose orbit is modified within a birth cluster.

Figure 7: Cumulative distributions of the predicted spin-orbit angle for
unstable systems (co-added across all simulation batches). Shown are unstable
systems that have lost their binary companion (blue) and unstable systems that
have retained their binary companion after 10 Myrs of integration (pink) and
110 Myrs of integration (green).

###  4.3 Ramifications for Exoplanets of Single Stars

As mentioned previously, planetary instabilities in our simulations are often
accompanied by the dissociation of the binary companion because both processes
are tied to strong cluster perturbations on the binary orbit. In fact,
$\sim$41% of our systems that had their binary companion stripped also
underwent a planetary orbital instability. Meanwhile, it has been argued that
most stars form with one or more stellar companions, but lose their companion
via dynamical processing (stellar encounters and tidal stripping) within their
birth cluster (e.g. Kroupa, 1995). Although only a relatively small portion of
our simulated systems undergo binary-triggered planetary instabilities (the
highest instability rate reached in our simulation sets is $\sim$15%), our
simulations do not encompass all parameter space of binaries or cluster
environments. Moreover, our simulations do not consider substructure within
our cluster potentials and make other idealized assumptions such as spherical
symmetry and instantaneous dissipation. If the actual physical environments of
clusters facilitate higher binary dissociation rates and these processes are
responsible for a large fraction of the isolated stars observed in the field,
a significant fraction (tens of percent) of planetary systems around isolated
stars may be subject to instabilities analogous to those in the simulations we
present here.

###  4.4 The Importance of Gas Potential

The bulk of the mass in our cluster environments is in the form of a Plummer
potential meant to represent the gaseous component. On the other hand, Figure
3 demonstrates that encounters with other cluster stars are critical in many
of the planetary system instabilities that our simulated binaries trigger, so
one may question whether the gas potential is important to the dynamics of
this process. To assess this, we perform two simulation sets (NonStellar_Fixed
and NonStellar_Wider) that are repeats of Medium_Fixed and Medium_Wider. Only
this time, they are performed in the NonStellar cluster environment in which
stellar masses (but not the Plummer mass) are reduced by an order of
magnitude.

In Table 3, we see that the rate of planetary instabilities is indeed lower in
both reduced stellar mass runs, but the effect is more severe if we use a
binary eccentricity distribution fixed at 0.5. Here we find that the
instability rate drops by over 75% (from 9% of all systems to 2%). In
contrast, when we employ wider binaries with a thermalized eccentricity
distribution there is a more modest decrease in instability rate of 44% (from
14.4% of all systems to 8%). The reason for this more modest change is that
our wider binary orbits are more easily perturbed and many more of them are
initialized with periastra near values necessary to destabilize our planetary
systems. Thus, more systems do not require powerful stellar encounters to
dramatically drive peristra low enough for a planetary instability, and this
can be accomplished with more modest gaseous perturbations. Thus, the
importance of the gas potential increases as one considers more eccentric and
more weakly bound binaries.

###  4.5 Comparison with Past Work

####  4.5.1 System Stability

As discussed in our introduction, there is an extensive literature studying
the stability of planets within binaries and star clusters. Many of the works
focused on the stability of planetary systems within clusters only consider
the direct effects of stellar flybys, and this limits the rate of planetary
instabilities. For instance, Adams et al. (2006) consider the rate that
stellar flybys perturb solar system architectures within embedded clusters of
100–1000 stars, and they find only $\sim$1% of Neptune-like planets around
Sun-like stars are destabilized (their Tables 3 and 7). Similarly, Proszkow &
Adams (2009) consider an even broader range of cluster parameters and
typically find destabilization rates of 0.1–1% for Neptune-like planets (their
Table 8). Meanwhile, we show here that the addition of a wide binary companion
boosts these destabilization rates by an order of magnitude within comparable
cluster environments. To achieve our binary-induced destabilization rates
around singleton stars requires much more extreme cluster environments. For
instance, Malmberg et al. (2011) find that typical open cluster environments
yield a destabilization probability of 5–15% for solar system-like
architectures, but their cluster stellar densities are 1–2 orders of magnitude
greater than ours, and their cluster lifetimes are 20 times as long.
Similarly, Hao et al. (2013) find cluster-generated stellar flybys can disrupt
99% of solar system-like architectures, but to achieve this rate requires
cluster densities $\sim$2 orders of magnitude greater than ours and 100-Myr
lifetimes. If we examine their disruption rates at 10 Myrs (their Figure 3),
it is only $\sim$2–10% with, again, a much denser cluster than any of ours.
Our simulations show that a wide binary companion acts as a potential conduit
allowing perturbations from the cluster environment to greatly boost the rate
of planetary instabilities.

Among works studying planetary stability within binaries, many do not consider
temporal evolution of the binary orbit. Holman & Wiegert (1999) developed a
widely employed empirical formula for the maximum semimajor axis, $a_{c}$, of
stable 1-planet systems in the presence of a coplanar wide binary companion
with a given mass and eccentricity:

$a_{c}=[(0.464\pm 0.006)+(-0.380\pm 0.010)\mu+(-0.631\pm 0.034)e\\\ +(0.586\pm
0.061)\mu e+(0.150\pm 0.041)e^{2}+(-0.198\pm 0.074)\mu e^{2}]a_{b}$ (14)

where $a_{b}$ and $e$ are the binary companion’s semimajor axis and
eccentricity, and $\mu$ is the ratio of the companion mass to the system mass.
This formula is approximate (e.g. Quarles & Lissauer, 2018), but if we take
Neptune’s semimajor axis of 30 au and note our fixed binary mass ratio of 1/3,
it allows us to invert Equation 1 and predict a “critical binary eccentricity”
above which our planetary system should be unstable and below which it will be
stable:

$e_{c}=\frac{.436-\sqrt{0.19-0.336(0.337-\frac{30}{a_{b}})}}{0.168}$ (15)

If we then examine our simulated initial binary eccentricities, we find that
the vast majority (97.1%) are started on orbits that will allow Neptune to
remain stable according to Equation 15. When we examine the 2.9% of binaries
whose initial eccentricities exceed $e_{c}$ we find that only 1/3 (or less
than 1% of our total systems) end up experiencing a planetary instability,
underscoring the fact that Equations 14 and 15 are approximate estimates of
stability. Thus, without considering cluster-driven binary evolution, a
stability analysis of our systems would predict nearly all of our systems to
be stable (hence their stability over an isolated integration). However, our
simulations record 395 planetary instabilities (9.6% of our integrated
systems), and 90% of these unstable systems begin with binary eccentricities
below $e_{c}$. Of course, our Medium_Uniform simulation is our only batch of
systems with an unskewed binary eccentricity distribution, but if we examine
only these systems, we find that 88% of our unstable systems begin below
$e_{c}$, and their median initial value is 0.83$e_{c}$.

Nearly half of our unstable systems lose their binary companion, but for those
that retain it, we can study how the final binary eccentricities compare to
$e_{c}$. In these systems, we find that the large majority ($\sim$80%) finish
with binary eccentricities above $e_{c}$. This demonstrates that it is
relatively rare in our simulations for binaries to evolve to a destabilizing
periastron and then evolve to a more modest periastron after triggering an
instability. Our cluster environments tend to disperse before this can occur.
Hence, final binary eccentricities relative to Equation 2 are a relatively
accurate marker of what types of planets could have been destabilized during
the cluster phase. However, these binary orbits should not be taken to
conclude that such planets could never have existed early in systems’
histories.

####  4.5.2 Binary-driven Kozai Cycles

Past work modeling the evolution of binary stars within clusters has noted
that cluster perturbations can excite the orbital inclinations of binary stars
and initiate Kozai cycles within planetary systems, possibly driving
instabilities and planetary eccentricity excitation. Parker & Goodwin (2009)
find that $\sim$20% of binaries can attain inclinations high enough ($i\gtrsim
40^{\circ}$) to excite planetary Kozai cycles. Our simulations find a
comparable fraction; 14.8% of our binary orbits are excited beyond to
$i>40^{\circ}$. However, this is not a guarantee of an instability within our
planetary systems. Of those systems whose binary orbit is excited beyond
$i\gtrsim 40^{\circ}$, 60% did not experience a planetary instability. This is
because Kozai cycles tend to be suppressed in multi-planet systems, and they
instead undergo rigid disk precession without substantial eccentricity
excitation (Innanen et al., 1997). Moreover, of our systems that undergo
instability and retain a binary, the final binary inclination is below 40∘ for
60% of systems. This underscores that Kozai cycles are not the primary driver
of instability and eccentricity excitation in our multi-planet systems.
Instead, it is the periastron evolution of the binary under cluster
perturbations (which is sometimes accompanied by significant binary
inclination evolution).

##  5 Conclusions

The simulations we present here demonstrate that wide (100 $\gtrsim
a_{b}\gtrsim$ 1000 au) binary stars likely have an early phase with active
temporal orbital evolution while they still inhabit their birth clusters. This
binary orbital evolution can in turn destabilize the planetary systems they
host, touching off episodes of planet-planet scattering. The instabilities we
document are typically triggered as cluster perturbations drive the binary
periastron toward the planets, exciting one or more planets onto crossing
orbits. Close encounters with other cluster stars are generally the cluster
perturbation that modifies the binary orbit, but the tidal potential of
cluster gas becomes increasingly relevant as we consider more eccentric and/or
more weakly bound binaries. Our simulations display a modest rate of
instabilities ranging from 5–15% of all planetary systems. However, due to
computational limitations, many regimes of the potential parameter space of
embedded star clusters are unexplored. It is possible that other cluster
environments as well as different binary and planetary system properties could
yield substantially higher rates of planetary instability.

The planetary instabilities that the binaries trigger often cause significant
alterations of the planetary architectures. Namely, the lower mass planets
(Uranus and Neptune analogs in our case) are preferentially lost, while the
surviving planets (typically Jupiter and Saturn analogs) remain with more
excited eccentricities. In addition, even though our simulations begin with
coplanar binaries, the inclinations of our planets are often altered as well.
This is because the same perturbations that alter binary periastra that
trigger instabilities also alter binary inclinations. These heightened binary
inclinations cause the planetary systems to undergo rigid disk precession
(Innanen et al., 1997), which can ultimately lead to high, even retrograde,
spin-orbit angles. Notably, it is rigid disk precession, rather than Kozai
cycles, that most often generate simulated systems with the highest spin-orbit
angles.

In systems with the highest spin-orbit angles, the binary companion that
torques the planets is usually retained, but overall, roughly half of our
unstable planetary systems also lose their binary, as the cluster
perturbations driving binary orbital evolution can lead to binary unbinding.
This raises the prospect that some exoplanets with high eccentricities and/or
spin-orbit angles around isolated field stars may be the products of cluster-
driven evolution of early binary companions that are lost early in the
systems’ histories.

##  6 Data Availability Statement

The simulation data underlying this article are available on the Harvard
Dataverse public repository.

##  7 Acknowledgements

This work was performed with support from NASA Exoplanets Research Program
grant 80NSSC19K0445, NSF grant AST-1814762, and NSF CAREER Award 1846388. Our
computing was performed at the OU Supercomputing Center for Education &
Research (OSCER) at the University of Oklahoma (OU).

## References

  * Adams (2010) Adams, F. C. 2010, ARA&A, 48, 47, doi: 10.1146/annurev-astro-081309-130830
  * Adams et al. (2006) Adams, F. C., Proszkow, E. M., Fatuzzo, M., & Myers, P. C. 2006, ApJ, 641, 504, doi: 10.1086/500393
  * Allen et al. (2007) Allen, L., Megeath, S. T., Gutermuth, R., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 361
  * Ambartsumian (1937) Ambartsumian, V. A. 1937, AZh, 14, 207
  * Bate (2000) Bate, M. R. 2000, MNRAS, 314, 33, doi: 10.1046/j.1365-8711.2000.03333.x
  * Bate et al. (2003) Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577, doi: 10.1046/j.1365-8711.2003.06210.x
  * Batygin et al. (2020) Batygin, K., Adams, F. C., Batygin, Y. K., & Petigura, E. A. 2020, AJ, 159, 101, doi: 10.3847/1538-3881/ab665d
  * Beust (2003) Beust, H. 2003, A&A, 400, 1129, doi: 10.1051/0004-6361:20030065
  * Beust et al. (2014) Beust, H., Augereau, J. C., Bonsor, A., et al. 2014, A&A, 561, A43, doi: 10.1051/0004-6361/201322229
  * Bonnell & Bate (1994) Bonnell, I. A., & Bate, M. R. 1994, MNRAS, 271, 999, doi: 10.1093/mnras/271.4.999
  * Brasser et al. (2012) Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icarus, 217, 1, doi: 10.1016/j.icarus.2011.10.012
  * Breslau & Pfalzner (2019) Breslau, A., & Pfalzner, S. 2019, A&A, 621, A101, doi: 10.1051/0004-6361/201833729
  * Cattolico & Capuzzo-Dolcetta (2020) Cattolico, R. S., & Capuzzo-Dolcetta, R. 2020, Astrophysics and Space Science, 365, doi: 10.1007/s10509-020-03885-4
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763–795, doi: 10.1086/376392
  * Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793, doi: 10.1046/j.1365-8711.1999.02379.x
  * Chambers et al. (2002) Chambers, J. E., Quintana, E. V., Duncan, M. J., & Lissauer, J. J. 2002, AJ, 123, 2884, doi: 10.1086/340074
  * Delva (1985) Delva, M. 1985, A&AS, 60, 277
  * Duncan et al. (1998) Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067, doi: 10.1086/300541
  * Dvorak (1986) Dvorak, R. 1986, A&A, 167, 379
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298, doi: 10.1086/521702
  * Haghighipour & Raymond (2007) Haghighipour, N., & Raymond, S. N. 2007, ApJ, 666, 436, doi: 10.1086/520501
  * Hamer & Schlaufman (2022) Hamer, J. H., & Schlaufman, K. C. 2022, arXiv e-prints, arXiv:2205.00040.  https://arxiv.org/abs/2205.00040
  * Hao et al. (2013) Hao, W., Kouwenhoven, M. B. N., & Spurzem, R. 2013, MNRAS, 433, 867, doi: 10.1093/mnras/stt771
  * Heppenheimer (1978) Heppenheimer, T. A. 1978, A&A, 65, 421
  * Holman & Wiegert (1999) Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621, doi: 10.1086/300695
  * Huang (1960) Huang, S.-S. 1960, PASP, 72, 106, doi: 10.1086/127489
  * Innanen et al. (1997) Innanen, K. A., Zheng, J. Q., Mikkola, S., & Valtonen, M. J. 1997, AJ, 113, 1915, doi: 10.1086/118405
  * Kaib et al. (2013) Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Nature, 493, 381, doi: 10.1038/nature11780
  * Kaib et al. (2011) Kaib, N. A., Raymond, S. N., & Duncan, M. J. 2011, ApJ, 742, L24, doi: 10.1088/2041-8205/742/2/L24
  * Kaib et al. (2018) Kaib, N. A., White, E. B., & Izidoro, A. 2018, MNRAS, 473, 470, doi: 10.1093/mnras/stx2456
  * Kouwenhoven et al. (2010) Kouwenhoven, M. B. N., Goodwin, S. P., Parker, R. J., et al. 2010, MNRAS, 404, 1835, doi: 10.1111/j.1365-2966.2010.16399.x
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Kroupa (1995) Kroupa, P. 1995, MNRAS, 277, 1491, doi: 10.1093/mnras/277.4.1491
  * Kroupa (2001) Kroupa, P. 2001, in The Formation of Binary Stars, ed. H. Zinnecker & R. Mathieu, Vol. 200, 199
  * Kroupa & Burkert (2001) Kroupa, P., & Burkert, A. 2001, ApJ, 555, 945, doi: 10.1086/321515
  * Lada & Lada (2003) Lada, C., & Lada, E. 2003, Annual Reviews of A&A
  * Lam & Kipping (2018) Lam, C., & Kipping, D. 2018, MNRAS, 476, 5692, doi: 10.1093/mnras/sty022
  * Levison et al. (2010) Levison, H. F., Duncan, M. J., Brasser, R., , & Kaufmann, D. E. 2010, Science Magazine, 329, doi: 10.1126/science.187535
  * Malmberg et al. (2011) Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859, doi: 10.1111/j.1365-2966.2010.17730.x
  * Parker & Goodwin (2009) Parker, R. J., & Goodwin, S. P. 2009, MNRAS, 397, 1041, doi: 10.1111/j.1365-2966.2009.15037.x
  * Parker et al. (2009) Parker, R. J., Goodwin, S. P., Kroupa, P., & Kouwenhoven, M. B. N. 2009, MNRAS, 397, 1577, doi: 10.1111/j.1365-2966.2009.15032.x
  * Plummer (1911) Plummer, H. C. 1911, MNRAS, 71
  * Proszkow & Adams (2009) Proszkow, E.-M., & Adams, F. C. 2009, ApJS, 185, 486, doi: 10.1088/0067-0049/185/2/486
  * Quarles & Lissauer (2018) Quarles, B., & Lissauer, J. J. 2018, AJ, 155, 130, doi: 10.3847/1538-3881/aaa966
  * Rasio & Ford (1996) Rasio, F. A., & Ford, E. B. 1996, Science, 274, 954, doi: 10.1126/science.274.5289.954
  * Reche et al. (2009) Reche, R., Beust, H., & Augereau, J. C. 2009, A&A, 493, 661, doi: 10.1051/0004-6361:200810419
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Sadavoy & Stahler (2017) Sadavoy, S. I., & Stahler, S. W. 2017, Monthly Notices of the Royal Astronomical Society, 469, 3881–3900, doi: 10.1093/mnras/stx1061
  * Siraj & Loeb (2020) Siraj, A., & Loeb, A. 2020, ApJ, 899, L24, doi: 10.3847/2041-8213/abac66
  * Turk et al. (2009) Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601, doi: 10.1126/science.1173540
  * Venner et al. (2021) Venner, A., Pearce, L. A., & Vanderburg, A. 2021, arXiv e-prints, arXiv:2111.03676.  https://arxiv.org/abs/2111.03676
  * Weidner & Kroupa (2004) Weidner, C., & Kroupa, P. 2004, Monthly Notices of the Royal Astronomical Society, 348, 187–191, doi: 10.1111/j.1365-2966.2004.07340.x
  * Wiegert & Holman (1997) Wiegert, P. A., & Holman, M. J. 1997, AJ, 113, 1445, doi: 10.1086/118360
  * Winn et al. (2010) Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010, The Astrophysical Journal, 718, L145, doi: 10.1088/2041-8205/718/2/l145
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, AJ, 102, 1528, doi: 10.1086/115978