# A Criterion for the Onset of Chaos in Compact, Eccentric Multiplanet Systems

Daniel Tamayo Department of Astrophysical Sciences, Princeton University, 4
Ivy Lane, Princeton, NJ 08540, USA  Norman Murray  Canadian Institute for
Theoretical Astrophysics, University of Toronto, 60 St. George Street,
Toronto, ON M5S 3H8, Canada  Scott Tremaine  Canadian Institute for
Theoretical Astrophysics, University of Toronto, 60 St. George Street,
Toronto, ON M5S 3H8, Canada  Institute for Advanced Study, Einstein Drive,
Princeton, NJ 08540, USA  Joshua Winn  Department of Astrophysical Sciences,
Princeton University, 4 Ivy Lane, Princeton, NJ 08540, USA

###### Abstract

We derive a semi-analytic criterion for the presence of chaos in compact,
eccentric multiplanet systems. Beyond a minimum semimajor-axis separation,
below which the dynamics are chaotic at all eccentricities, we show that (i)
the onset of chaos is determined by the overlap of two-body mean motion
resonances (MMRs), like it is in two-planet systems; (ii) secular evolution
causes the MMR widths to expand and contract adiabatically, so that the
chaotic boundary is established where MMRs overlap at their greatest width.
For closely spaced two-planet systems, a near-symmetry strongly suppresses
this secular modulation, explaining why the chaotic boundaries for two-planet
systems are qualitatively different from cases with more than two planets. We
use these results to derive an improved angular-momentum-deficit (AMD)
stability criterion, i.e., the critical system AMD below which stability
should be guaranteed. This introduces an additional factor to the expression
from Laskar & Petit (2017) that is exponential in the interplanetary
separations, which corrects the AMD threshold toward lower eccentricities by a
factor of several for tightly packed configurations. We make routines for
evaluating the chaotic boundary available to the community through the open-
source SPOCK package.

††thanks: NHFP Sagan Fellow<EMAIL_ADDRESS>
##  1 Introduction

Radial-velocity surveys (e.g., Udry et al., 2019) and transit surveys
(Lissauer et al., 2011) have both revealed a substantial population of compact
multiplanet systems. Various authors have defined the term “compact” in
different ways (e.g., Lissauer et al., 2011; Volk & Gladman, 2015). We will
use a dynamically motivated definition: a system (or subsystem) of planets is
compact if the period ratio of any adjacent pair of planets is between 1 and
$2$. Within this period-ratio range, first-order mean-motion resonances
(hereafter MMRs), with period ratios of $j:j-1$, strongly modify the dynamical
structure of phase space; at wider separations (and low orbital
eccentricities), the strongest resonances are typically secular (between the
orbits’ apsidal and nodal precession frequencies). The restriction in the
range of period ratios also allows us to simplify many algebraic expressions
using the approximation $|\Delta a|/a\ll 1$, where $a$ refers to the semimajor
axis. By this definition, Venus, Earth and Mars are a compact trio, and Uranus
and Neptune are a compact pair. In the exoplanet context, $\approx\,40\%$ of
3+ planet systems have at least one pair of adjacent planets that is compact.

Exoplanets typically complete $10^{10}\mathchar 45\relax 10^{12}$ orbits
between the time they form and the time we observe them. Given these vast
spans of time, a central question for connecting observations to planet-
formation models is whether planetary populations look the same after 10 Gyr
as they do following disk dispersal at ages of a few Myr. Several studies have
used N-body experiments to argue that multiplanet systems evolve significantly
over their lifetimes, starting with larger numbers of planets and subsequently
destabilizing and colliding to form progressively longer and longer lived
configurations with fewer and fewer planets (e.g., Pu & Wu, 2015; Volk &
Gladman, 2015). However, our theoretical understanding of the dynamics driving
these instabilities remains incomplete.

The ability to predict the stability of arbitrary orbital configurations would
also help to constrain uncertain planetary masses and orbital parameters. This
argument is based on the “temporal Copernican principle” — a corollary of
which is that an observer examining a planetary system at an arbitrary moment
in the system’s lifetime should not expect to detect it just prior to a major
cataclysm (e.g., Gott, 1993). In other words, one would expect a randomly
discovered system to have comparable past and future lifetimes. One can
therefore reject values of planet masses and orbital parameters for which the
system’s estimated future lifetime would be much shorter than its current age,
making it possible to constrain planetary masses and orbital parameters that
are challenging to measure directly (Tamayo et al., 2021). This approach has
occasionally been implemented by directly checking for stability over
timescales of $\sim 10^{9}-10^{10}$ orbits through brute-force numerical
integrations for some of the most important exoplanet discoveries (e.g.,
Steffen et al., 2013; Tamayo et al., 2015, 2017; Quarles et al., 2017; Wang et
al., 2018; Rosenthal et al., 2019). However, the Gyr ages of typical systems
mean that $N$-body integrations are orders of magnitude too computationally
expensive for a full exploration of parameter space, even for a single system.

Many authors, starting with Chambers et al. (1996), have used N-body
integrations to explore the dependence of instability times on the initial
conditions and parameters of multiplanet systems (Yoshinaga et al., 1999;
Marzari & Weidenschilling, 2002; Faber & Quillen, 2007; Zhou et al., 2007;
Smith & Lissauer, 2009; Funk et al., 2010; Pu & Wu, 2015; Obertas et al.,
2017; Gratia & Lissauer, 2019; Lissauer & Gavino, 2021). While these numerical
experiments have provided insight into the dynamical mechanisms that drive
these instabilities, the practical restrictions that were variously imposed on
the allowed parameter space (e.g., equal spacing, equal orbital
eccentricities, equal masses, etc.) lead to inaccurate instability time
predictions for real systems, often by orders of magnitude (Cranmer et al.,
2021).

Tamayo et al. (2016) and Tamayo et al. (2020) tackled this high-dimensional
prediction problem for compact multiplanet systems using machine learning
techniques. They demonstrated that their Stability of Planetary Orbital
Configurations Klassifier (SPOCK) can accurately assess the stability of
compact systems over $10^{9}$ orbits based on an analysis of a much shorter
integration. This results in computation times up to $10^{5}$ times faster
than direct N-body integrations. Recently, Cranmer et al. (2021) generalized
this result with a deep learning model that can not only classify a system’s
stability, but also predict the time to instability out to $10^{9}$ orbits.

In this paper we pursue the problem of predicting the long-term stability of
compact planetary systems analytically. Such efforts also have a long history,
and it has become clear that distinct dynamical mechanisms limit the stability
of planetary systems in different regions of parameter space. In order to put
our work into context and to introduce several of the concepts and results we
will use below, we begin by summarizing some of the previous key results.

##  2 Background

A central challenge in understanding the stability of compact $N$-planet
system is that there are $7N$ parameters to consider (6 orbital elements and 1
mass per planet). Analytical work has therefore progressed by starting from
simple cases and adding parameters sequentially. One powerful simplification
we will use throughout this paper is to consider only the planar case. As we
motivate numerically below, this is an excellent approximation for mutual
inclinations less than a few degrees, typical of transiting multiplanet
systems (e.g., Fabrycky et al., 2014).

###  2.1 Two-Planet Systems

The first analytic investigations considered the case of two bodies on
initially circular orbits. Wisdom (1980) showed that, at the smallest
separations, chaos is driven by the overlap of first-order ($j:j-1$)
MMRs111Throughout this paper, we will use the term ‘order’ to describe the
power of the eccentricity $e$ carried by a particular term. The strength of a
$j:j-k$ two-body MMR is proportional to $e^{k}$ (e.g., Murray & Dermott,
1999). First-order ($k=1$) MMRs are thus the strongest two-body MMRs for small
eccentricities.. The intervals between the period ratios of adjacent first-
order MMRs become smaller and smaller as $j$ increases (the 2:1 is further in
period ratio from the 3:2 than the 3:2 is from the 4:3, etc.). At the same
time, the widths of these ever more tightly spaced resonances increase due to
the decreasing interplanetary separation (even at zero eccentricity). Wisdom
(1980) showed that this implies a threshold interplanetary separation below
which the MMR widths overlap, leading to chaotic motion,

$\Bigg{(}\frac{a_{2}-a_{1}}{a_{1}}\Bigg{)}_{\rm{\\!\\!crit}}\approx
1.46\,\mu^{2/7},$ (1)

where the $a_{i}$ are the semimajor axes, and $\mu$ is the planet-star mass
ratio. This critical spacing, originally derived in the restricted case of a
planet, star, and test particle, can be generalized to the case of two planets
with non-zero masses by replacing the mass ratio $\mu$ by the sum of the
planetary mass ratios (Deck et al., 2013),

$\mu=\frac{m_{1}+m_{2}}{M_{\star}},$ (2)

where the $m_{i}$ represent the planetary masses, and $M_{\star}$ the stellar
mass.

Wisdom (1980) also showed that MMR widths grow with increasing orbital
eccentricities. Thus, for given masses and separations beyond the critical
separation for cirular orbits (Eq. 1), there exists a critical eccentricity
beyond which first-order MMRs will overlap and drive chaos. However, while
higher order MMRs have negligible widths in the initially circular case
(Ferraz-Mello, 2007) and thus can be ignored, Hadden & Lithwick (2018) show
that they must be considered at higher eccentricities where their widths can
be substantial.

The results of Hadden & Lithwick (2018) are based on two key developments that
facilitate a simple analytic criterion for resonance overlap.

##### The relative eccentricity

It has long been known that the dynamical evolution near a first-order MMR
does not depend on both planets’ eccentricity vectors independently, but only
on a particular linear combination of the two (Sessin & Ferraz-Mello, 1984;
Wisdom, 1986). Hadden (2019) showed that this result also holds to excellent
approximation for higher order MMRs in the compact limit $|\Delta a/a|\ll 1$.

In particular, at the closest separations—where the equations of motion reduce
to Hill’s equations (Hénon & Petit, 1986)—there are two additional conserved
quantities (e.g., Goldreich & Tremaine, 1981), which at low eccentricities can
be thought of as a “center-of-mass” eccentricity vector. Defining eccentricity
vectors in the complex plane $\boldsymbol{e}_{i}=e_{i}\exp{(i\varpi_{i}})$
with magnitude given by the eccentricity $e_{i}$ and direction given by the
longitude of pericenter $\varpi_{i}$ (henceforth vectors will be denoted in
boldface), the additional conserved quantity in the Hill problem is, to
leading order in the eccentricities,
$\boldsymbol{e_{com}}=(m_{1}\boldsymbol{e_{1}}+m_{2}\boldsymbol{e_{2}})/(m_{1}+m_{2})$.
Therefore, in the Hill problem the dynamics do not depend on the two
eccentricity vectors independently but rather only on their vector difference
$\boldsymbol{e_{12}}\equiv\boldsymbol{e_{2}}-\boldsymbol{e_{1}}\equiv
e_{12}\exp(i\varpi_{12}).$ (Hénon & Petit, 1986).

The expansion of the disturbing function of the planet interaction potential
includes several “sub-resonances” corresponding to a $j:j-k$ MMR, which are
proportional to
$e_{1}^{l}e_{2}^{k-l}\cos\left[j\lambda_{2}-(j-k)\lambda_{1}-(k-l)\varpi_{2}-l\varpi_{1}\right]$
for $l=0,\dots,k$. As argued by Hadden (2019), the fact that in the Hill limit
the disturbing function can only depend on the eccentricities through the
combinations $e_{12}$ and $\varpi_{12}$ implies that the sum over these sub-
resonances must be expressible as a single resonant term proportional to
$e_{12}^{k}\cos\left[j\lambda_{2}-(j-k)\lambda_{1}-k\varpi_{12}\right]$.

While these results strictly apply only in the Hill limit, Hadden (2019)
showed that they persist to a good approximation222Hadden (2019) derive
slightly modified definitions of the center-of-mass eccentricity and
eccentricity difference that depend on the interplanetary separation and
provide an even better approximation. However, in order to combine the widths
of different MMRs, which lie at different interplanetary separations, we
follow Hadden & Lithwick (2018) in ignoring these corrections. far beyond the
regime where Hill’s equations apply, and in fact hold reasonably well for all
period ratios $\lesssim 2$.

##### The orbit-crossing eccentricity

In general, MMR strengths depend on both the interplanetary separation and the
orbital eccentricities. In a second important development, Hadden & Lithwick
(2018) showed that in the compact limit ($|\Delta a/a|\ll 1$), the dependence
on the semimajor axes can be absorbed by expressing eccentricities as a
fraction of their orbit-crossing value $e_{12}^{\rm{cross}}$

$e_{12}^{\rm{cross}}=\frac{a_{2}-a_{1}}{a_{2}}=1-\Bigg{(}\frac{j-k}{j}\Bigg{)}^{\\!\\!2/3}\approx\frac{2k}{3j},$
(3)

We therefore define the normalized, relative eccentricity vector333This is
equivalent to $\mathcal{Z}_{12}/\mathcal{Z}_{12}^{\rm{cross}}$ in the notation
of Hadden & Lithwick (2018), though we choose
$e_{12}=\sqrt{2}\mathcal{Z}_{12}$ and
$e_{12}^{\rm{cross}}=\sqrt{2}\mathcal{Z}_{12}^{\rm{cross}}$ for simplicity of
notation and geometric interpretation.

$\boldsymbol{\tilde{e}_{12}}\equiv\frac{\boldsymbol{e_{2}}-\boldsymbol{e_{1}}}{e_{12}^{\rm{cross}}}.$
(4)

With this definition, $\boldsymbol{\tilde{e}_{12}}=1$ is the collision
criterion for two arbitrarily oriented, coplanar orbits in the compact limit
(Kholshevnikov & Vassiliev, 1999; Laskar & Petit, 2017).

In the pendulum approximation, which is valid if the eccentricity is not too
close to zero (e.g., Murray & Dermott, 1999), the half-width $\delta a$ of a
$k^{\mathrm{th}}$ order MMR is (Hadden & Lithwick, 2018)

$\frac{\delta a}{a}\approx
4\sqrt{\frac{h_{k}}{3}}\mu^{\frac{1}{2}}\tilde{e}^{\frac{k}{2}},$ (5)

plus terms of higher order in the eccentricity. Here $h_{k}$ is a constant
specific to each set of $k$th order MMRs (and is of order unity for modest
values of $k$), and $\mu$ and $\tilde{e}$ are defined in Eqs. 2 and 4,
respectively.

Figure 1: Different chaotic boundaries imposed by various dynamical
mechanisms, and studies deriving those criteria. Yellow regions are chaotic,
while green regions are not The $x$-axis is the dynamical spacing between
adjacent planets (Eq. 6), while the $y$-axis is the linear combination of
eccentricity vectors that dominantly sets the resonant dynamics (Eq. 4). The
left panel shows the two-planet case, with schematic MMR widths that start
overlapping at the red boundary (Eq. 8). The right panel shows an equal-mass,
equal-spacing, equal-eccentricity multiplanet case, where 3-body MMRs and
secular evolution additionally can play a major role.

Hadden & Lithwick (2018) consider the aggregate sum of all MMR widths, and
calculate a covering fraction or optical depth $\tau$ of MMRs within a fixed
range of interplanetary separations. A key parameter, which we will call the
“dynamical separation” $S_{12}$ between a pair of planets, is given by

$S_{12}=\frac{a_{2}-a_{1}}{R_{1/4}},\>{\rm where}\>R_{1/4}\equiv
a_{2}\,\mu^{1/4}$ (6)

and $\mu$ is defined by Equation (2)444Another way to portray $S$ is as a
normalized crossing eccentricity, $e^{\rm{cross}}/\mu^{1/4}$. We note that $S$
does not increase without bound for wide interplanetary separations, but
instead approaches a constant value of $\mu^{-1/4}$.. They derive a covering
fraction of two-body MMRs between the pair of planets,

$\tau_{12}\approx\Bigg{(}\frac{1.8}{S_{12}}\Bigg{)}^{\\!2}|\ln(\tilde{e}_{12})|^{-3/2}.$
(7)

Solving for the critical eccentricity at which $\tau$ equals unity and MMRs
are overlapped, Hadden & Lithwick (2018) obtain a chaotic boundary for two
planets

$\tilde{e}_{12}^{\rm{crit}}\approx\exp{\Bigg{[}-\Bigg{(}\frac{1.8}{S_{12}}\Bigg{)}^{\\!\\!4/3\,}\Bigg{]}}\;\;(2\>\rm{planets}).$
(8)

In the left panel of Fig. 1, we schematically plot the widths of various MMRs
in gray, which grow with increasing eccentricity. The boundary at which the
covering fraction of MMRs reaches unity is plotted in red (Hadden & Lithwick,
2018). In summary, the chaos in two-planet systems is always driven by the
overlap of two-body MMRs. At the closest separations, first-order MMRs overlap
even at zero eccentricity (Wisdom, 1980). At wider separations, MMRs only
overlap at non-zero eccentricity, and higher order MMRs need to be considered
(Hadden & Lithwick, 2018).

We can also gain some intuition from the following approximations (discussed
and compared to the above expressions in Appendix B), which are valid at tight
separations where $S\lesssim 1.8$,

$\displaystyle\tau_{12}$ $\displaystyle=$
$\displaystyle\Bigg{(}\frac{2.9}{S_{12}}\Bigg{)}^{\\!\\!2}\tilde{e}_{12}\>\>\>\>\>\>\>\>\>\>\>\>\>\>(\tilde{e}\lesssim
0.4)$ $\displaystyle\tilde{e}_{12}^{\rm{crit}}$ $\displaystyle=$
$\displaystyle\Bigg{(}\frac{S_{12}}{2.9}\Bigg{)}^{\\!\\!2}\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>(\tilde{e}\lesssim
0.4),$ (9)

From these equations and Eq. 6, we see that $\tau$ scales approximately
quadratically with the fractional interplanetary separation, linearly with the
eccentricity, and as the square root of the mass ratio (and has a steeper
scaling with eccentricity for $\tilde{e}>0.4$). Because the stability
criterion depends more strongly on eccentricity than mass ratio, when the
requirement of long-term stability is applied to exoplanet systems, the
resulting eccentricity constraints are typically stronger than the constraints
on masses (Tamayo et al., 2021).

In addition, the three-body problem allows for a rigorous criterion whereby
close encounters are forbidden beyond a critical value of a dimensionless
parameter combining the total angular momentum and energy (Marchal & Bozis,
1982). This general result was applied to planetary systems by Gladman (1993),
where it is commonly referred to as the Hill criterion. At low eccentricities
(see Petit et al. 2018 for the eccentric two-planet case), close encounters
are forbidden for interplanetary separations $|\Delta a|>2\sqrt{3}\,R_{H}$,
where $R_{H}$ denotes the mutual Hill radius $R_{H}=a_{1}(\mu/3)^{1/3}$.
Partly because of the Hill criterion, dynamical separations in compact systems
are traditionally expressed as a multiple of the mutual Hill radius. However,
the Hill criterion is not valid for systems with more than two planets, while
resonance overlap criteria can be generalized to higher multiplicities.

We may therefore expect $S$ (Eq. 6) to be a more useful measure of dynamical
separation than the Hill separation in multiplanet systems. In practice, the
differences between the scaling of the Hill criterion, $|\Delta
a|\sim\mu^{1/3}$, the resonance-overlap criterion for circular orbits,
$|\Delta a|\sim\mu^{2/7}$, and the dynamical separation, $|\Delta
a|\sim\mu^{1/4}$, are small enough to be challenging to distinguish in
$N$-body experiments. However, comparing the ends of the planet-mass
distribution, the differences become significant. A value of $S=2.5$ that
limits $\tilde{e}$ to be below $\approx 0.5$ (Eq. 8) corresponds to a Hill
separation of $\approx 10R_{H}$ for a pair of Earth-mass planets around the
Sun, while it translates to $\approx 6R_{H}$ for two Jupiter-mass planets.

###  2.2 Initially Circular Multiplanet Systems

The instability boundary in two-planet systems is sharp, with configurations
largely divisible into those that go unstable on timescales comparable to the
synodic period between planet conjunctions, and those that remain stable over
typical planetary system lifetimes (Gladman, 1993). By contrast, a large
number of numerical results starting with those of Chambers et al. (1996) have
shown that initially circular systems of three or more planets exhibit a much
broader dynamic range of instability timescales.

Quillen (2011), and more recently Petit et al. (2020), explain this behavior
as a consequence of the overlap of three-body resonances. In the range of
period ratios between first and second order MMRs (these MMRs have finite
widths at zero eccentricity and such configurations often exhibit much shorter
instability times in N-body integrations, e.g., Obertas et al. 2017), this is
a compelling explanation. Both secular resonances and higher order MMRs have
vanishing widths at zero eccentricity (Ferraz-Mello, 2007) — and therefore the
overlap of those resonances cannot be driving the instabilities observed in
numerical integrations.

So far we have been considering a perturbative expansion in powers of the
(typically small) eccentricities, but at leading (first) order in the mass
ratio $\mu$. In this limit, the strongest resonances are first-order two-body
MMRs, with strengths $\propto\mu\tilde{e}$ (e.g., Murray & Dermott, 1999) and
widths $\propto\mu^{1/2}\tilde{e}^{1/2}$ (Eq. 5). At the next power in the
masses $\propto\mu^{2}$, however, there appear three-body MMRs involving the
orbital parameters of trios of planets. These can again be separated into
three-body MMRs of different orders in the eccentricity. Unlike the case of
two-body MMRs (excluding co-orbital resonances), there exist three-body MMRs
that are zeroth-order in (i.e., independent of) the eccentricity (e.g.,
Quillen, 2011). While these three-body MMRs are therefore weaker by a factor
of $\mu$ than two-body MMRs, their independence of the eccentricity renders
them important for circular orbits. Indeed, Quillen (2011) showed that the
finite widths of such three-body MMRs at zero eccentricity allow them to
overlap and drive chaos at interplanetary separations that are significantly
wider than the two-planet criterion (Eq. 1) of Wisdom (1980). Building on that
work by accounting for variable chaotic diffusion rates in different regions
of the three-body resonance network, Petit et al. (2020) were able to
accurately predict the instability times measured in numerical integrations of
such systems.

While the instability times in regions where three-body resonances overlap
span a broad range of timescales, they are almost always smaller than the
typical main-sequence lifetimes ($\sim\,10^{11}$ orbits) of the known stars
hosting compact multiplanet systems. The critical separation inside which
three-body resonances overlap therefore provides a sufficient criterion for
instability in typical planetary systems. Petit et al. (2020) derive this
critical separation for the general case in which the masses and
interplanetary spacings between a trio of planets are unequal. For the simple
case of equal masses and separations (EMS), the results of Quillen (2011) and
Petit et al. (2020) can be written as555There is a typographical error in Eq.
61 of Petit et al. 2020 (A. Petit, personal communication).

$S_{12\rm{crit}}\approx 1.68\>\>(\rm{3\mathchar 45\relax
body\>MMR\>overlap,\>EMS}).$ (10)

In the case of EMS systems, one can write the 3-body MMR overlap criterion in
terms of the dynamical separation $S$ between adjacent planets (with $\mu$
given by twice the planet-star mass ratio through Eq. 2), but in the general
case the criterion depends on both of the separations between adjacent
planets, and $\mu$ is replaced by a generalized mass parameter which is not
simply the sum of the mass ratios (see Eqs. 45 and 53 of Petit et al., 2020).

Figure 2: Distribution of dynamical spacings (Eq. 6) from the transiting
multiplanet sample of Weiss et al. (2018), with masses estimated using a mass-
radius relationship. Histograms are separated into systems with two known
planets (blue) and those with three or more known planets (red).

Figure 2 shows the distribution of dynamical spacings between adjacent planets
for the 909 planets in 355 multi-planet systems analyzed by Weiss et al.
(2018) as part of the California Kepler Survey (Petigura et al., 2017).
Calculating these spacings requires adopting planet masses, which Weiss et al.
(2018) estimate using the mass-radius relations from Weiss & Marcy (2014). We
plot separate histograms for systems with two known planets and ones with
higher planet multiplicities. Fewer than 2% of adjacent planet pairs in 3+
planet systems have $S<1.68$. More careful analysis of these most closely
packed systems using the general expressions of Petit et al. (2020) for
unequal spacings and masses reveal that three-body resonances do not overlap
in these systems.

This is unsurprising given that such configurations would have destabilized
and rearranged long ago. However, the observed distribution of dynamical
spacings seems suppressed even beyond $S=1.68$, perhaps hinting that many more
widely spaced configurations can nevertheless be unstable. In particular,
while initially circular orbits may be stable beyond $S=1.68$, eccentric
orbits need not be. This expectation is consistent with numerical experiments
showing that stability is a strong function of eccentricities (e.g., Zhou et
al., 2007; Gratia & Lissauer, 2019).

This motivates trying to understand the general case, particularly given that
we expect planets to have eccentric orbits. Terrestrial planet formation
models predict eccentric orbits (e.g., Hansen & Murray, 2012; Dawson et al.,
2016; Tremaine, 2015), and transit timing variation (Lithwick et al., 2012;
Jontof-Hutter et al., 2014; Hadden & Lithwick, 2016, 2017) and transit
duration measurements (Van Eylen & Albrecht, 2015; Xie et al., 2016; Van Eylen
et al., 2019) confirm this expectation.

###  2.3 Eccentric Multiplanet Systems (this paper)

At the outset of this work, it was not clear to us whether the stability of
most multiplanet configurations at wider separations remains limited by the
overlap of two-body MMRs, as in the two-planet case, or by eccentricity-
dependent, higher order 3-body MMRs. A clue to the underlying dynamics is that
the machine learning model SPOCK (Tamayo et al., 2020) was only trained on
three-planet systems in which one pair was in (or near) a strong 2-body MMR,
and yet the trained model made accurate predictions even for randomly sampled
compact systems. As argued by Tamayo et al. (2020), this suggests that the
chaotic dynamics in eccentric compact systems is dominated by two-body MMRs
and not 3-body MMRs.

We took an approach that is qualitatively similar to that of Petit et al.
(2017), who exploited the wide separation of timescales between the fast
resonant dynamics and the much slower secular modulation of the eccentricities
$\tilde{e}$. In particular, this slow modulation causes the eccentricity-
dependent MMR widths (Eq. 5) to adiabatically grow and shrink over a secular
cycle (e.g., Wisdom, 1985), which in turn affects whether or not MMRs overlap.
However, our work differs in three main ways.

The first and simplest difference is that we consider the secular modulation
of a different MMR boundary. Petit et al. (2017) consider whether the secular
evolution can raise the eccentricities to the point at which first-order MMRs
overlap (or to $\approx\,e_{\rm{cross}}$ in cases where the latter is smaller;
Laskar & Petit 2017). However, Hadden & Lithwick (2018) showed that the
overlap of MMRs of all orders provides a more accurate match to the chaotic
boundary. We therefore instead consider the modulation of the two-planet
boundary of Hadden & Lithwick (2018), which, as we will see below, often
yields chaotic boundaries that are a factor of two or more lower in the
eccentricities than the estimates of Petit et al. (2017) for compact systems.

Second, Petit et al. (2017) and Laskar & Petit (2017) consider the long-term
redistribution of the system’s angular momentum deficit (AMD; Laskar, 2000)

${\rm
AMD}=\sum_{j=1}^{N}m_{j}\sqrt{GM_{\star}a_{j}}\Bigg{(}1-\sqrt{1-e_{j}^{2}}\cos{i_{j}}\Bigg{)},$
(11)

where $G$ is the gravitational constant, the $m_{j}$ are the planetary masses,
the $e_{j}$ the orbital eccentricities, the $i_{j}$ the orbital inclinations,
and the sum runs over all $N$ planets. In particular, they consider the worst-
case scenario, in which the system exhibits secular chaos which transfers all
of the system’s AMD to the $\tilde{e}$ of a single pair of planets. The
advantage of this estimate based on the worst-case scenario is simplicity: for
given planetary masses and semimajor axes, it generates a single system-wide
critical AMD value, below which the system should always be stable. However,
it is not clear how to interpret the criterion in the limit in which one body
is massless (in which case its eccentricity would always be able to reach
unity) or in the limit of a distant planetary or stellar companion (which can
raise the AMD to arbitrarily large values and would allow the eccentricities
of all the inner bodies to reach unity).

Additionally, if most compact systems were not significantly affected by
secular chaos, then the threshold eccentricities that guarantee AMD stability,
while technically correct, would be too low, hedging against a possibility
that seldom occurs in practice. While secular chaos in our more widely
separated solar system drives major changes in Mercury’s orbit on timescales
of $10^{10}$–$10^{11}$ orbits (Laskar & Gastineau, 2009), and secular chaos
has been demonstrated for some widely separated giant planets (Wu & Lithwick,
2011) as well as some compact multiplanet configurations (Volk & Malhotra,
2020), it is also plausible that for most low-mass compact multiplanet
systems, the requirement to avoid 2-body MMR overlap constrains the
eccentricities to low enough values that secular resonances do not overlap and
secular chaos is not important.

Rather than trying to address this difficult question, we restrict ourselves
in this paper to deriving a chaos criterion for two-body MMR overlap in
compact multiplanet systems in the presence of eccentricity modulation by the
shorter-term regular, i.e., non-chaotic, secular dynamics. This will allow us
to check our derived boundaries against the results from computationally
tractable $N$-body integrations and to establish whether the shortest chaotic
timescales are set by two-body MMR overlap, or whether other resonances, e.g.,
higher order 3-body MMRs, need to be considered.

The third difference between the present work and Petit et al. (2017) and
Laskar & Petit (2017) is that we consider not only the most tightly packed
planet pair but all of the planet pairs. We show that in compact systems it is
necessary to consider two-body MMRs between different pairs of planets in the
system simultaneously.

We begin in Sec. 3 by exploring the secular modulation of $\tilde{e}$ in the
two-planet case, and showing why this effect is negligible in the case studied
by Hadden & Lithwick (2018). We consider multiplanet systems in Sec. 4, first
analyzing the effect of two-body MMRs with multiple bodies in the system (Sec.
4.1), and then isolating the effects of the secular dynamics (Sec. 4.3). We
put everything together in Sec. 5, and conclude in Sec. 6.

##  3 Two Planets

To approach the problem analytically, we exploit the hierarchy of timescales
in the problem. Because the secular timescales are significantly longer than
both the orbital periods and the MMR libration timescales, we average over the
latter two fast timescales to isolate the secular evolution. To lowest order
in the eccentricities, this is the classical Laplace-Lagrange problem, in
which the semimajor axes remain constant. While we ignore MMR contributions to
the precession rates, which can be significant near and inside MMRs (Wisdom,
1985; Sansottera & Libert, 2019), we will find that this simplest
approximation is sufficient in most cases for determining the chaotic
boundary.

The planar Laplace-Lagrange problem with two planets has four dynamical
variables corresponding to each planet’s orbital eccentricity and pericenter
orientation, with the total energy and angular momentum as conserved
quantities. Thus the two-planet Laplace-Lagrange system has two phase-space
dimensions or one degree of freedom. In the compact limit, there is an
additional conserved quantity $\boldsymbol{e_{com}}$, which removes the
remaining freedom for the orbits to exchange angular momentum, so the
eccentricities remain fixed (i.e., one of the Laplace-Lagrange
eigenfrequencies goes to zero). While this conserved quantity only exists in
the limit of the Hill problem, i.e., at the closest separations, the vestiges
of this symmetry persist across the period ratio range where strong MMRs exist
($P_{2}/P_{1}\lesssim 2$), strongly suppressing secular eccentricity
oscillations — as we show below. We argue that this is one of the fundamental
reasons why the chaotic boundaries for two-planet and three-planet systems are
qualitatively different from one another.

###  3.1 Coordinates and Notation

We develop our Hamiltonian treatment in canonical Poincaré variables, using
democratic heliocentric coordinates (see Duncan et al., 1998; Hernandez &
Dehnen, 2017). In that case the actions and their conjugate angles are

$\displaystyle\Lambda_{i}=m_{i}\sqrt{GM_{\star}a_{i}}$
$\displaystyle\leftrightarrow$ $\displaystyle\lambda_{i}$ (12)
$\displaystyle\Gamma_{i}=\Lambda_{i}\Bigg{(}1-\sqrt{1-e_{i}^{2}}\Bigg{)}$
$\displaystyle\leftrightarrow$ $\displaystyle\gamma_{i}=-\varpi_{i}$ (13)

where the $\lambda_{i}$ are the mean longitudes, and the $\varpi_{i}$ are the
longitudes of pericenter. Because we average over the fast MMR oscillations
the problem is purely secular, we may assume that the semimajor axes and their
associated actions $\Lambda_{i}$ remain fixed, while the $\Gamma_{i}$ evolve
with time. We also define complex variables

$\boldsymbol{G_{i}}\equiv\sqrt{2\Gamma_{i}}e^{i\gamma_{i}}.$ (14)

###  3.2 Laplace-Lagrange Hamiltonian

The secular Hamiltonian for two planets is

$\mathcal{H}=-\frac{Gm_{1}m_{2}}{a_{20}}\mathcal{R}_{sec},$ (15)

where, to leading order in the eccentricities (e.g., Murray & Dermott, 1999),

$\mathcal{R}_{sec}=\frac{\alpha}{8}\Bigg{[}b_{3/2}^{1}(\alpha)(e_{1}^{2}+e_{2}^{2})-2b_{3/2}^{2}(\alpha)e_{1}e_{2}\cos(\varpi_{2}-\varpi_{1})\Bigg{]},$
(16)

where $\alpha=a_{1}/a_{2}$ is the semimajor axis ratio, and the $b_{s}^{j}$
are Laplace coefficients evaluated at $\alpha$.

The secular disturbing function can be written in matrix notation in terms of
the complex eccentricity vectors $\boldsymbol{e_{i}}$ as (e.g., Hadden, 2019),

$\mathcal{R}_{sec}=\frac{\alpha}{8}\begin{pmatrix}\boldsymbol{e_{1}}^{*}&\boldsymbol{e_{2}}^{*}\end{pmatrix}\begin{pmatrix}b_{3/2}^{1}(\alpha)&-b_{3/2}^{2}(\alpha)\\\
-b_{3/2}^{2}(\alpha)&b_{3/2}^{1}(\alpha)\end{pmatrix}\begin{pmatrix}\boldsymbol{e_{1}}\\\
\boldsymbol{e_{2}}\end{pmatrix}.$ (17)

For $e_{1},e_{2}\ll 1$ the eccentricity vectors are connected to the canonical
variables through

$\begin{pmatrix}\boldsymbol{e_{1}}^{*}\\\
\boldsymbol{e_{2}}^{*}\end{pmatrix}\simeq\begin{pmatrix}\sqrt{\frac{1}{\Lambda_{1}}}&0\\\
0&\sqrt{\frac{1}{\Lambda_{2}}}\end{pmatrix}\begin{pmatrix}\boldsymbol{G_{1}}\\\
\boldsymbol{G_{2}}\end{pmatrix}.$ (18)

Plugging Eq. 18 into Eq. 17 yields the secular Hamiltonian

$\mathcal{H}=\frac{1}{2}\begin{pmatrix}\boldsymbol{G_{1}}&\boldsymbol{G_{2}}\end{pmatrix}\cdot{\mathcal{M}}\cdot\begin{pmatrix}\boldsymbol{G_{1}^{*}}\\\
\boldsymbol{G_{2}^{*}}\end{pmatrix},$ (19)

where, after some algebra,

$\mathcal{M}=-\frac{n_{2}}{4M_{\star}}\begin{pmatrix}\alpha^{1/2}\,b_{3/2}^{1}(\alpha)\,m_{2}&-\alpha^{3/4}\,b_{3/2}^{2}(\alpha)\,\sqrt{m_{1}m_{2}}\\\
-\alpha^{3/4}\,b_{3/2}^{2}(\alpha)\,\sqrt{m_{1}m_{2}}&\alpha\,b_{3/2}^{1}(\alpha)\,m_{1}\end{pmatrix}$
(20)

and $n_{2}$ is the outer planet’s mean motion $\sqrt{GM_{\star}/a_{2}^{3}}$.
We retain factors of $\alpha$ in the secular evolution to later consider the
secular influence of distant perturbers on tightly spaced planets (Sec.
4.3.1). For cases considered below with more than two planets, we evaluate the
matrix $\mathcal{M}$ numerically using the open-source celmech
package666https://github.com/shadden/celmech..

The Laplace-Lagrange equations of motion in canonical variables read

$\frac{d}{dt}\begin{pmatrix}\boldsymbol{G_{1}}\\\
\boldsymbol{G_{2}}\end{pmatrix}=-i\mathcal{M}\cdot\begin{pmatrix}\boldsymbol{G_{1}}\\\
\boldsymbol{G_{2}}\end{pmatrix}.$ (21)

Since the matrix $\mathcal{M}$ is real and symmetric, this set of linear
differential equations can be solved by finding a rotated basis in which the
matrix $\mathcal{D}=\mathcal{R}\cdot\mathcal{M}\cdot\mathcal{R}^{T}$ is
diagonal. The variables rotated into this basis

$\begin{pmatrix}\boldsymbol{A_{1}}\\\
\boldsymbol{A_{2}}\end{pmatrix}=\mathcal{R}\cdot\begin{pmatrix}\boldsymbol{G_{1}}\\\
\boldsymbol{G_{2}}\end{pmatrix}$ (22)

are the Laplace-Lagrange modes, whose amplitudes $A_{m}$ remain fixed, and
rotate with frequencies given by the corresponding eigenvalues of
$\mathcal{D}$. The evolution of the $\boldsymbol{G_{i}}$ can thus be
decomposed into linear combinations of these modes, which each rotate
uniformly with constant amplitude.

Physically, the components of the matrix $\mathcal{R}_{mp}$ represent how much
of mode $\boldsymbol{A_{m}}$ is in the planet with index $p$. These matrix
elements are set by the planetary masses and semimajor axes, while the initial
eccentricities and pericenter orientations set the conserved mode amplitudes
$A_{m}$ as well as the initial phases (Eq. 22).

Assuming there are no secular resonances between these frequencies, the
maximum value of $G=|\boldsymbol{G}|$ for the planet with index $p$ can then
be simply estimated as the sum of the contributions to the planet from each of
the Laplace-Lagrange mode amplitudes,

$\Big{(}G_{p}\Big{)}_{\rm{max}}=\sum_{m=1}^{N}|\mathcal{R}_{pm}|A_{m}.$ (23)

###  3.3 A Canonical Rotation

To estimate the widths of MMRs, we specifically require the eccentricity
combination
$\boldsymbol{\tilde{e}_{12}}\propto\boldsymbol{e_{2}}-\boldsymbol{e_{1}}$ (Eq.
4). We can study the evolution of this particular combination of
eccentricities through a mass-dependent canonical rotation (Wisdom, 1986;
Batygin & Morbidelli, 2013; Hadden, 2019). Defining

$\begin{pmatrix}\boldsymbol{F_{1}}\\\
\boldsymbol{F_{2}}\end{pmatrix}=\mathcal{R}_{\psi}\cdot\begin{pmatrix}\boldsymbol{G_{1}}\\\
\boldsymbol{G_{2}}\end{pmatrix},$ (24)

where $\mathcal{R}$ represents a (4-dimensional) counterclockwise rotation
through an angle
$\psi=\pi+\sin^{-1}(\Lambda_{1}/(\Lambda_{1}+\Lambda_{2}))^{1/2}$:

$\mathcal{R}_{\psi}=\begin{pmatrix}\cos{\psi}&-\sin{\psi}\\\
\sin{\psi}&\cos{\psi}\end{pmatrix}=-\begin{pmatrix}\sqrt{\frac{\Lambda_{2}}{\Lambda_{1}+\Lambda_{2}}}&-\sqrt{\frac{\Lambda_{1}}{\Lambda_{1}+\Lambda_{2}}}\\\
\sqrt{\frac{\Lambda_{1}}{\Lambda_{1}+\Lambda_{2}}}&\sqrt{\frac{\Lambda_{2}}{\Lambda_{1}+\Lambda_{2}}}\end{pmatrix}$
(25)

One can check through direct substitution from Eq. 18 into Eq 24 that

$\displaystyle\boldsymbol{\tilde{e}_{12}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\Lambda_{1}+\Lambda_{2}}{\Lambda_{1}\Lambda_{2}}}\frac{\boldsymbol{F_{1}}^{*}}{e_{12}^{\rm{cross}}}$
$\displaystyle\boldsymbol{e_{com}}$ $\displaystyle\approx$
$\displaystyle\frac{\Lambda_{1}\boldsymbol{e_{1}}+\Lambda_{2}\boldsymbol{e_{2}}}{\Lambda_{1}+\Lambda_{2}}=-\frac{\boldsymbol{F^{*}_{2}}}{\sqrt{\Lambda_{1}+\Lambda_{2}}},$
(26)

where the approximate equality in the second line holds when $\alpha\approx 1$
in the compact limit. Thus, the new complex canonical coordinates
$\boldsymbol{F}$ are scaled versions of the center-of-mass eccentricity and
the $\boldsymbol{\tilde{e}_{12}}$ we require to estimate the MMR widths.

This yields a Hamiltonian in the new variables

$\displaystyle\mathcal{H^{\prime}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\boldsymbol{F_{1}}&\boldsymbol{F_{2}}\end{pmatrix}\cdot{\mathcal{M}^{\prime}}\cdot\begin{pmatrix}\boldsymbol{F_{1}^{*}}\\\
\boldsymbol{F_{2}^{*}}\end{pmatrix}$ (27) $\displaystyle\mathcal{M}^{\prime}$
$\displaystyle=$
$\displaystyle\mathcal{R}_{\psi}\cdot{\mathcal{M}}\cdot\mathcal{R}_{\psi}^{T};$

we will call $\mathcal{M}^{\prime}$ the rotated Laplace-Lagrange matrix.

###  3.4 A Nearly Conserved Quantity for 2 Planets

We now show that the rotation $\mathcal{R}_{\psi}$ from Poincaré variables to
$\boldsymbol{F}$ is close to the matrix that diagonalizes $\mathcal{M}$. This
implies that the combinations $\boldsymbol{\tilde{e}_{12}}$ and
$\boldsymbol{e_{\rm{com}}}$ are approximately the conserved secular eigenmodes
of the system.

After some algebraic manipulation, using Eqs. 20 and 25, the rotated Laplace-
Lagrange matrix can be written as

$\mathcal{M}^{\prime}=-\begin{pmatrix}\omega_{12}&0\\\
0&0\end{pmatrix}+d\omega_{12}\mathcal{M}_{d},$ (28)

where, since we will show that $d$ is a small parameter, $\omega_{12}$ is
approximately the eigenvalue of $\boldsymbol{F_{1}}$, and thus close to the
precession rate of $\boldsymbol{\tilde{e}_{12}}$ (Eq 26),

$\omega_{12}=\Bigg{(}\frac{\alpha\,m_{1}+\alpha^{1/2}\,m_{2}}{M_{\star}}\Bigg{)}\frac{b_{3/2}^{1}(\alpha)}{4}\,n_{2}.$
(29)

The parameter $d$ is the fractional difference between the two Laplace
coefficients appearing in $\mathcal{M}$,

$d=\frac{b_{3/2}^{1}(\alpha)-b_{3/2}^{2}(\alpha)}{b_{3/2}^{1}(\alpha)}$ (30)

and

$\mathcal{M}_{d}=\begin{pmatrix}2\frac{\Lambda_{1}\Lambda_{2}}{(\Lambda_{1}+\Lambda_{2})^{2}}&\frac{\sqrt{\Lambda_{1}\Lambda_{2}}}{\Lambda_{1}+\Lambda_{2}}\frac{\Lambda_{1}-\Lambda_{2}}{\Lambda_{1}+\Lambda_{2}}\\\
\frac{\sqrt{\Lambda_{1}\Lambda_{2}}}{\Lambda_{1}+\Lambda_{2}}\frac{\Lambda_{1}-\Lambda_{2}}{\Lambda_{1}+\Lambda_{2}}&-2\frac{\Lambda_{1}\Lambda_{2}}{(\Lambda_{1}+\Lambda_{2})^{2}}\end{pmatrix}$
(31)

has elements that are $\mathcal{O}(1)$ or smaller.

Using an approximation introduced by Goldreich & Tremaine (1981), we find

$b_{3/2}^{j}(\alpha)\approx\frac{2|j|}{\pi(1-\alpha)}K_{1}(|j|(1-\alpha)),$
(32)

where $K_{1}$ is a modified Bessel function of the second kind. Given that our
expressions only involve $j=1$ and $j=2$, in the compact limit the argument
$j(1-\alpha)\ll 1$, so we can use the asymptotic expansion near zero for
$K_{1}(z)\approx z^{-1}$, yielding

$b_{3/2}^{j}(\alpha)\approx\frac{2}{\pi(1-\alpha)^{2}},\>j(1-\alpha)\ll 1,$
(33)

independent of $j$. In fact, we find Eq. 33 provides a closer approximation to
$b_{3/2}^{2}$ than Eq. 32 for $\alpha>0.75$, and for $b_{3/2}^{1}$ for
$\alpha>0.33$, within 10% of the exact value. In this approximation, $d$
vanishes, so $d\lesssim 0.1$ for $\alpha>0.75$, and drops to zero as $\alpha$
approaches unity (for reference $\alpha>0.63$ inside the 2:1 MMR).

The vanishing second eigenvalue in Eq. 31 implies that both the magnitude and
direction of $\boldsymbol{F_{2}}\propto\boldsymbol{e_{\rm{com}}}$ are
conserved in the limit where $d$ is negligible, as we know must be the case in
the compact limit from the discussion in the introduction to this section. At
finite separations, the conservation is not perfect, but the oscillations in
$\boldsymbol{F_{1}}$ are attenuated by the factor $d<0.1$, which suppresses
variations in the widths of MMRs.

This near symmetry allowed Hadden & Lithwick (2018) to consider the resonant
dynamics of two-planet systems ignoring any secular evolution. However,
additional planets break this symmetry and render the secular evolution
important.

##  4 Multiplanet Systems

Two separate effects must be considered for multiplanet systems. As mentioned
above, one must calculate how much the secular eccentricity oscillations will
modulate the stability boundary. Second, one must account for MMRs with the
additional planets in the system. To better illustrate the dynamics, in what
follows, we will first isolate each of these effects. Then, in Sec. 5, we will
put them together.

###  4.1 MMRs with Additional Planets

The approach by Hadden & Lithwick (2018) of calculating an optical depth of
MMRs between a pair of planets is straightforward to generalize to the case
with additional planets. For a given planet, we now show that the simplest
approximation—simply adding the MMR optical depths (Eq. 7) with all other
planets—provides a reliable estimate of the chaotic boundary. We begin by
considering the case of three planets, and will argue in Sec. 4.2 that
considering adjacent trios of planets in higher multiplicity systems is a
useful simplification. A new complication for three-planet systems is that the
MMR optical depths of the middle planet with each of its adjacent neighbors
depend on both $\tilde{e}_{12}$ and $\tilde{e}_{23}$, and these will generally
evolve with time due to the secular dynamics.

We can nevertheless isolate the interaction of MMRs between different pairs of
planets and visualize the results by considering a three-planet system
initialized in a single Laplace-Lagrange eigenmode. Exciting only a single
eigenmode implies that, at the level of our leading-order approximation, all
eccentricities, and thus $\tau_{12}$ and $\tau_{23}$, remain constant in time.

We seem to be left with four dimensions: the two dynamical separations between
adjacent planets (Eq. 6), and $\tilde{e}_{12}$ and $\tilde{e}_{23}$ (Eq. 4).
However, the ratios between the various orbital eccentricities (and thus
$\tilde{e}_{12}/\tilde{e}_{23}$) are actually fixed by the eigenmode
components, and only depend on the planetary masses and orbital periods (e.g.,
Murray & Dermott, 1999). If we additionally fix the initial orbital periods of
the innermost and outermost planets, then choosing the middle planet’s orbital
period $P_{2}$ sets both $S_{12}$ and $S_{23}$. The system can then specified
by two parameters $(S_{12},\tilde{e}_{12})$, just as in the two-planet problem
considered by Hadden & Lithwick (2018), and we can similarly calculate a one-
dimensional optical depth, summing the contributions from each pair of
adjacent planets.

In order to have both adjacent planet pairs contribute comparably to the MMR
optical depth, we consider a system of three $10M_{\oplus}$ planets. We first
consider a tightly spaced case, with the inner planet period fixed at
$P_{1}=1$, and the outer period fixed at $P_{3}=2.033$ (with the particular
value chosen to avoid a strong MMR between in the innermost and outermost
planet). We then vary the middle planet’s orbital period in the range
$P_{2}=[1.33,1.5]$, which approximately symmetrically covers the region
between the 4:3 and 3:2 MMR with both pairs of planets. For each simulation,
we solve for the Laplace Lagrange eigenmodes at the given masses and
separations. We then initialize the system purely in the eigenmode with the
fastest eigenfrequency, with the appropriate amplitude to yield the value of
$e_{12}/e_{\rm{cross}}$ plotted on the y-axis in Fig. 3. The eigenmode
components allow us to compute the appropriate (and different) initial
eccentricity vector for the third planet.

For each configuration we run a numerical integration using the WHFast
integrator in the REBOUND N-body package, using a fixed timestep of $5\%$ of
the pericenter passage timescale, in order to accurately integrate eccentric
orbits (Wisdom, 2015). In each case we calculate the Laplace-Lagrange
eigenfrequencies. In compact configurations, one of the Laplace-Lagrange
eigenfrequency is much slower than the other ones, reflecting the near
conservation of the system’s center-of-mass eccentricity. We therefore choose
to run the integrations for three times the period of the next slowest
eigenfrequency of the system, which corresponds to $\sim 10^{4}$ orbital
periods for the masses and separations we consider. We then measure the Mean
Exponential Growth of Nearby Orbits (MEGNO) chaos indicator (Cincotta et al.,
2003), which converges to a value of two for regular orbits and grows beyond
two at a rate given by the fastest Lyapunov exponent for chaotic orbits. Like
any chaos indicator, the MEGNO can only identify chaotic trajectories with
Lyapunov timescales shorter than the time of integration. Since in all our
plots we seek to study the boundary of trajectories that are chaotic on
secular timescales (which we use to set the timespan of the integrations), we
stretch the color map to highlight small deviations from two. However, the
results are visually indistinguishable even if we extend the color bar by an
order of magnitude. We explore this further and justify our choice of
integration times in Appendix A.

The top panel of Fig. 3 shows the resulting distribution of MEGNO values.
Yellow regions are chaotic, while dark green regions are regular. In purple we
plot contours for the two-planet optical depths from Hadden & Lithwick (2018)
for each adjacent pair of planets, showing the critical eccentricity at which
the MMR covering fraction reaches unity. In black we show the contour where
the sum of the two-planet optical depths crosses one. The agreement on the
left of the plot seems poor due to two wide first-order resonances (the 4:3
between the inner two planets at $P_{2}\approx 1.33$, and the 3:2 between the
outer two planets at $P_{2}\approx 1.36$), but one finds regular (green)
configurations between them for eccentricities approaching the black line
where the two MMRs start crossing, as expected.

In the bottom panel of Fig. 3 we show an identical experiment but with a wider
spacing. The inner planet remains at $P_{1}=1$, but the outer planet has been
moved outward to $P_{3}=3.03$, while the middle planet now approximately
symmetrically samples the range between the 3:2 and 2:1 MMR with both its
neighbors. We see that considering the total MMR optical depth for the middle
planet from its MMRs with both its adjacent neighbors provides a good
approximation to the chaotic boundary.

Figure 3: The contour where the total MMR covering fraction with all pairs of
planets crosses unity (black) provides an approximate boundary for chaotic
motion. These integrations involved three 10 $M_{\oplus}$ planets. The x-axis
shows the middle planet’s orbital period; the outer planets’ periods were held
fixed (top panel $P_{1}=1$, $P_{3}=2.033$, bottom panel $P_{1}=1$,
$P_{3}=3.03$). To suppress any secular evolution and highlight the MMR optical
depth criterion, the eccentricity vectors for the three planets were chosen to
excite only one Laplace-Lagrange eigenmode, and to yield the value of
$\tilde{e}_{12}$ on the y-axis. Purple lines denote the critical eccentricity
contours at which the optical depth of MMRs between the middle planet and one
of its neighbors would reach unity.

Both panels in Fig. 3 use the general expression for the MMR optical depth
(Eq. 7) from Hadden & Lithwick (2018). While the low-eccentricity
approximation (Eq. 9) provides nearly indistinguishable results for the three-
planet criterion in the top panel, it significantly underestimates the optical
depths in the lower panel where the normalized eccentricities in the regular
regions extend $\gtrsim 0.4$, resulting in a boundary that is too high.

###  4.2 Considering Only Adjacent Trios of Planets

When generalizing to even higher multiplicities, we now show that a useful
further simplification is to only consider adjacent trios of planets. Given
that $\tau$ scales only weakly with mass (Eq. 9), we imagine for simplicity a
system of equal-mass planets, and adopt equal spacings as the worst-case
scenario for ignoring the contributions from non-adjacent planets.

Taking the eccentricities as approximately equal in such a five-planet case,
the quadratic scaling of $\tau$ with $S$ (Eq. 9) implies that the middle
planet would have a $\sim 25\%$ contribution to its total MMR optical depth
from non-adjacent planets that are twice as far away as the adjacent planets.
In more general cases, the correction to only considering the most dynamically
packed trio can be significantly smaller. Additionally, once the period ratio
with non-adjacent planets exceeds 2, one is no longer in the compact limit and
Eq. 9 significantly overestimates the optical depth due to the lack of first-
order MMRs. Thus, while for low-planet masses and exotic configurations one
can construct counterexamples (e.g., a closely packed system with two Jupiters
on the ends and two Earths in between), it is generally sufficient to consider
only adjacent planets.

Figure 4: Considering MMR overlap only for the tightest trio of planets
approximates the chaotic boundary. Same numerical experiment as the top panel
of Fig. 3, but with two additional planets, one interior and one exterior to
the previously considered system (see text). The orbital periods are chosen as
a “worst-case scenario” such that the middle trio has approximately the same
spacing as the innermost trio and outermost trio on the left and right of the
plot, respectively. The eccentricity vectors are similarly initialized to only
excite a single Laplace-Lagrange eigenmode and suppress the secular evolution.
The MMR overlap criterion considering only the middle trio of planets (black)
is within $\approx 25\%$ of the optical depth considering all adjacent pairs
of planets (cyan) and the numerical results in this worst-case scenario.

In Fig. 4, we reproduce the numerical experiment in the top panel of Fig. 3,
except we insert two additional $10M_{\oplus}$ planets at orbital periods of
$0.63$ and $3.08$. These values were chosen to ensure that the middle trio of
planets is always the most tightly packed one (the middle trio becomes
comparably spaced to the innermost and outermost trio on the left and right
sides of the plot, respectively). As above, we solve for the Laplace-Lagrange
eigenmodes and initialize the system in a single mode to suppress the secular
evolution.

We see that the result is qualitatively similar, with our MMR overlap
criterion considering only the middle trio of planets (in black) within
$\approx 25\%$ of the optical depth considering all adjacent pairs of planets
(cyan) as motivated above. The strongly chaotic region on the lower left of
the top panel of Fig. 3, caused by the wide first-order MMRs in that period
range as discussed above, is expanded in Fig. 4 by the additional planets.
Apart from this widening, the chaotic structure is qualitatively similar to
Fig. 3. Given that we designed this numerical experiment as a worst-case
scenario for the adjacent-planet approximation, we conclude that only
considering the trio of adjacent planets with the highest $\tau$ provides a
useful reduction in the dimensionality for high multiplicity systems. However,
we mention that even if non-adjacent planets do not contribute meaningfully to
the MMR optical depth directly, they can still significantly alter the secular
evolution, as we now explore.

###  4.3 Secular Evolution

Next, we isolate the effects of secular evolution. We consider the effect of
adding a third body to a system of two tightly packed planets. We choose to
place the third body at a period ratio $>2$ with the remaining planets,
suppressing the optical depth of two-body MMRs it contributes due to the lack
of first-order MMRs beyond the 2:1 (Sec. 4.1). This allows us to focus solely
on the effects of the secular evolution.

We begin by estimating the separation at which one can ignore the secular
effects of unseen planets. This is an important question when trying to draw
conclusions about the close-in systems we predominantly discover, given our
difficulties in observationally detecting distant planets.

####  4.3.1 When Do Distant Unseen Planets Matter?

We found in Sec. 3.3 that in the limit where the third planet can be ignored,
$\boldsymbol{F_{1}}$ is approximately an eigenmode of the system. This secular
problem has been considered in a variety of contexts, most recently in
understanding the excitation of mutual inclinations by an exterior perturber
(Becker & Adams, 2016; Lai & Pu, 2017)777 While we are interested in the
eccentricity rather than the mutual inclination evolution, Pu & Lai (2018)
shows that in both cases one can use the same dimensionless parameter to
describe the coupling strength of the third body.. Rather than rederive the
same results in a different way, we restate them in the context of this
approach.

There are two separate regimes. At wide separations, whenever the ratio
$\epsilon$ of the differential precession induced by the distant perturber on
the inner two orbits is much smaller than the mutual precession induced by the
inner two orbits on one another (Eq. 29), the inner two orbits precess
together rigidly, i.e., at the same average rate (e.g., Lai & Pu, 2017). In
this widely separated limit, $\boldsymbol{e_{12}}$, which tracks the
differential evolution of the inner two planets’ eccentricity vectors, is an
eigenmode of the secular dynamics, and circulates with constant amplitude.
$\boldsymbol{F_{2}}$ tracks the “rigid-body” evolution of the inner system,
and will oscillate in amplitude, but has no effect on the MMR widths (Hadden,
2019). The widths of the MMRs between the compact pair therefore remain
constant and the distant perturber has no effect on the stability boundary.

When the ratio $\epsilon$ exceeds unity, the perturber disrupts the rigid-body
evolution of the inner system. In the limit we are interested in, the inner
two planets are compact ($\alpha_{12}\approx 1$) and the Laplace coefficients
can be approximated through Eq. 33. Inserting Eq. 33 into the expression given
by Lai & Pu (2017) and setting $\alpha_{12}\approx 1$, we have

$\epsilon=\frac{9\pi}{4}\Bigg{(}\frac{m_{p}}{m_{1}+m_{2}}\Bigg{)}\Bigg{(}\frac{P_{2}}{P_{p}}\Bigg{)}^{2}e_{12\rm{cross}}^{3},\>\>\rm{(Outer}\>\rm{Perturber)}$
(34)

where $P_{2}$ and $P_{p}$ are the periods of the middle and perturbing planet,
and $e_{12\rm{cross}}$ is the crossing eccentricity of the inner two planets
(Eq. 3). We note that Denham et al. (2019) derived the a similar condition
under which a third body can strongly influence the secular dynamics of an
inner pair of planets, but they were focused on the regime where all planets
are well separated. In our case where the inner pair of orbits are compact
(period ratio $<2$), Eq. 34 differs from the criterion of Denham et al. (2019)
by over an order of magnitude.

Finally, we mention that at short orbital periods, general relativity and
other short-range forces can also suppress secular oscillations (e.g., Liu et
al., 2015; Wei et al., 2021; Faridani et al., 2021). For closely packed
systems, it is again important to consider whether the differential precession
on the two orbits exceeds their mutual gravitationally induced precession.

####  4.3.2 Secular Modulation of the Stability Boundary

Figure 5: Control experiment of a three-planet system that follows the
2-planet chaotic boundary. Two 10 $M_{\oplus}$ planets with inner period P1
fixed at 1, and P2 varied along the x-axis. There is a 1 $M_{\rm{Jup}}$
perturber always at a period P3=3.2. The eccentricities and pericenters are
chosen to excite only one Laplace-Lagrange mode, suppressing the secular
evolution. The MMRs between the middle and outer planet (shown in dashed
lines) do not meaningfully contribute to the total optical depth of MMRs. This
three-planet system therefore still exhibits the two-planet stability boundary
(purple) from Hadden & Lithwick (2018) corresponding to the Jupiter-mass
planet not being present.

The MMR widths between a compact pair of planets, and thus the two-planet
stability criterion of Hadden & Lithwick (2018), depend only on the
eccentricity combination $\tilde{e}_{12}$ (Eq. 4). As the MMR widths expand
and contract through the secular eccentricity evolution induced by a distant
third planet, the task is therefore to determine the maximum value
$e_{12}^{\rm{max}}$ in the secular cycle. it is thus useful to analyze the
problem in the same rotated basis as in Sec. 3.3,

$\begin{pmatrix}\boldsymbol{F_{1}}\\\ \boldsymbol{F_{2}}\\\
\boldsymbol{F_{3}}\end{pmatrix}=\begin{pmatrix}\cos{\psi}&-\sin{\psi}&0\\\
\sin{\psi}&\cos{\psi}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\boldsymbol{G_{1}}\\\ \boldsymbol{G_{2}}\\\
\boldsymbol{G_{3}}\end{pmatrix},$ (35)

where $\boldsymbol{F_{1}}\propto\boldsymbol{\tilde{e}_{12}}$,
$\boldsymbol{F_{2}}\propto\boldsymbol{e_{\rm{com}}}$ (Eq. 26), and
$\boldsymbol{F_{3}}=\boldsymbol{G_{3}}$.

Figure 6: Secular evolution can strongly modulate the chaotic boundary. Same
planet masses and separations as in Fig. 5, except now the eccentricities and
pericenters are initialized so that $\boldsymbol{e_{12}}$ is made up of two
eigenmodes with initial conditions depicted in the top left panel. The radius
of the black circle given by the vector sum of the two mode contributions is
$e_{12}$, and the units on the axes are the factors by which $e_{12}$ changes
from its original value. While the initial $e_{12}$ changes on the y-axis of
the right panels, the mode amplitudes are chosen so that $e_{12}$ always
changes by the same factor of three. The top right panel shows that estimating
MMR overlap using the initial eccentricities does not match the chaotic
boundary. The bottom left panel shows the maximum in the eccentricity cycle,
and the corresponding bottom right panel shows that the MMR widths at the
secular maximum do provide an accurate overlap boundary.

We consider a suite of configurations with two tightly packed $10M_{\oplus}$
planets around a solar-mass star. The first planet is always initialized at an
orbital period $P_{1}=1$, while we vary $P_{2}$ from the 4:3 to the 3:2 MMR
with the inner planet. In order to avoid considering MMRs with the perturbing
third planet, in all configurations we place the perturber beyond the 2:1
resonance, at a period $P_{3}=3.2$. In order for the perturber to drive
differential precession in the inner pair of planets, we assign it a Jupiter
mass. This corresponds to $\epsilon\approx 0.1-0.3$ (Eq. 34) in the period
range we consider. In other words, the differential precession induced by the
perturber mixes $\approx 10-30\%$ of the center-of-mass degree of freedom
$\boldsymbol{F_{2}}$ into the eigenmode that is almost purely made up of
$\boldsymbol{e_{12}}$ in the absence of the perturber (Sec. 3). This is enough
to drive substantial modulation of the MMR widths depending on initial
conditions, as we show below.

As a first control experiment, we suppress the secular evolution by
initializing the eccentricities in a single eigenmode as in Sec. 4.1,
specifically in the eigenmode associated with $e_{12\rm{com}}$. We vary
$e_{12}$ along the y-axis of Fig. 5, while the eigenmode components set the
ratio of $e_{12\rm{com}}$ and $e_{3}$ to $e_{12}$. Over our range of periods
for the middle planet, $e_{12\rm{com}}/e_{12}$ varies from $\approx 0.8-2.25$,
while planet three remains nearly circular, with $e_{3}/e_{12}$ varying from
$\approx 0.02-0.07$. Thus, the MMR widths should stay fixed in time, and the
prediction is that the stability boundary should follow the two-planet
stability boundary of Hadden & Lithwick (2018).

We plot the results in Fig. 5. In solid black lines, we plot the widths of all
MMRs up to tenth order between the compact pair of planets, with widths given
by $\frac{\delta P}{P}\approx\frac{3}{2}\frac{\delta a}{a}$, and the latter
given by Eq. 5. In dashed black lines, we plot the widths of MMRs between the
middle and outer planet, which we see have negligible widths.

There are few regular trajectories beyond the two-planet criterion in purple
from Hadden & Lithwick (2018). The most obvious exception is inside strong
MMRs, like the 4:3 on the leftmost side of the plot. We reiterate that our
approximations do not cover such cases—certainly the dynamics in those regions
will not be well described by the Laplace-Lagrange secular dynamics we
described above. Many authors have considered the secular dynamics inside MMRs
(Wisdom, 1985; Morbidelli & Moons, 1993; Batygin & Morbidelli, 2013; Pichierri
et al., 2017; Hadden, 2019), but we do not pursue this further.

By engineering a configuration in which MMRs with the third planet can be
ignored, and by initializing the eccentricity vectors so as to suppress the
secular evolution, we set up a three-planet system that follows the two-planet
stability boundary despite the large mass and relative proximity of the
perturbing planet.

Figure 7: Mutual inclinations within a few degrees do not modify the chaotic
boundary. Same numerical experiment as in Fig. 5, but with all planets
assigned the same orbital inclination (given above each subplot), and a random
longitude of ascending node. The panels thus correspond to maximum mutual
inclinations of $1.1^{\circ},3.4^{\circ}$ and $11^{\circ}$, respectively.

In our second experiment we adopt the same masses, separations and initial
$e_{12}$ as above, but this time we initialize the $e_{12}$ as a mix of the
two eigenmodes associated with the compact pair of planets. In particular, for
each configuration, we first solve for the rotated Laplace-Lagrange eigenmodes
$\boldsymbol{A_{m}}$. For a given value of $e_{12}$ we plot on the y axis of
Fig. 6, we then choose the mode amplitudes such that $\boldsymbol{A_{1}}$
contributes $2e_{12}$ to the inner pair’s relative eccentricity, and
$\boldsymbol{A_{2}}$ contributes $e_{12}$. We then start the modes anti-
aligned so the total eccentricity is $2e_{12}-e_{12}=e_{12}$ as desired (top
left panel of Fig. 6). We see on the top right panel that the MMR widths are
too low at the initial value of $e_{12}$ to explain the chaotic regions
through their overlap.

However, once the system starts evolving, and the eigenmodes rotate at their
own eigenfrequencies, $e_{12}$ and thus the resonance widths will change with
time. Given our setup, all configurations will reach a value of $e_{12}$ that
is three times higher than the initial value at the maximum in the secular
cycle when the modes align (bottom left panel). In the bottom right panel we
plot the MMR widths using this inflation factor of three for the
eccentricities. We also plot the boundary from Hadden & Lithwick (2018) for
the critical eccentricity depressed by the same factor of three.

In this picture, the MMR widths are adiabatically modulated by the secular
evolution, and the MMR criterion of Hadden & Lithwick (2018) needs to be
modified so that it does not refer to the critical initial eccentricity, but
rather the maximum eccentricity attained during a secular cycle.

###  4.4 Ignoring Inclinations

Throughout the paper we have worked assuming co-planar orbits. We now examine
the role of small inclinations in modifying the chaotic boundary. We reproduce
the numerical experiment in Fig. 5, with all planets starting with the same
orbital inclination and randomly assigned longitudes of ascending node. The
results are shown in Fig. 7, with inclinations of $0.01,0.03$ and $0.1$
radians as labeled above each subplot. Given the randomly assigned node
longitudes, these correspond to maximum mutual inclinations between pairs of
planets of $1.1^{\circ},3.4^{\circ}$ and $11^{\circ}$, respectively. We see
that inclinations within a few degrees, typical of transiting multiplanet
systems (e.g., Fabrycky et al., 2014) do not affect the chaotic boundary,
although mutual inclinations $\gtrsim 5^{\circ}$ would need to be taken into
account.

##  5 Discussion

We now put everything together. We run a numerical experiment corresponding to
the right panel in Fig. 1. We consider three $10M_{\oplus}$ planets on equally
spaced orbits and equal initial eccentricities (varied along the y-axis). The
initial orbital pericenters are selected at $120^{\circ}$ intervals so as to
drive secular evolution by exciting all the secular eigenmodes.

Figure 8: The chaotic boundary in multiplanet systems. Shown are results from
simulations of three equally spaced $10M_{\oplus}$ planets, as a function of
the dynamical spacing (Eq. 6) and the normalized eccentricities (Eq. 4); see
the text for details. At the tightest separations (left of solid red line),
chaos is driven by the overlap of 3-body MMRs (Quillen, 2011; Petit et al.,
2020). At wider separations, the chaotic boundary is set by the overlap of
2-body MMRs (solid black line). We compare this boundary to the limit for two-
body MMR overlap in 2-planet systems from Hadden & Lithwick (2018), shown in
dashed purple. The difference between the black and purple lines is due to the
secular evolution of the eccentricities driven by the third planet. We also
show in dashed blue the critical eccentricity at which the Laplace-Lagrange
secular evolution would achieve the orbit-crossing condition of Laskar & Petit
(2017).

The result is shown in Fig. 8. First, we see that in this equally spaced,
equal-mass case, dynamical spacings $\lesssim 1.68$ are chaotic at all
eccentricities due to the overlap of 3-body MMRs as predicted by Quillen
(2011). As shown by Petit et al. (2020) and various numerical experiments
(e.g., Faber & Quillen, 2007; Smith & Lissauer, 2009; Obertas et al., 2017),
the resulting instabilities occur within $\sim 10^{10}$ orbits, which is
significantly shorter than most exoplanet hosts’ main-sequence lifetimes. In
cases with unequal masses and/or period ratios, where adjacent planets have
different dynamical spacings, this boundary for 3-body MMR overlap would shift
toward smaller dynamical spacing of the tightest planet pair $<1.68$, with
general expressions given by Petit et al. (2020).

We have argued that beyond the 3-body MMR overlap limit, chaos on secular
timescales is driven by the overlap of 2-body MMRs, and that this boundary can
be significantly modulated by the system’s regular (non-chaotic) secular
evolution. The purple dashed line shows the two-planet limit for the overlap
of 2-body MMRs from Hadden & Lithwick (2018). Our boundary for this case is in
black, with approximately half the deviation from Hadden & Lithwick (2018)
arising from considering MMRs with both of the middle planet’s adjacent
neighbors (Sec. 4.1), and the other half coming from the secular evolution
(Sec. 4.3). In this case, our eccentricity boundary differs from the two-
planet criterion of Hadden & Lithwick (2018) by a factor of $\approx 2$,
though this difference can be larger or smaller depending on the system’s
secular mode amplitudes (Sec. 4.3.2).

We can see that while our boundary matches the numerical results well at tight
separations, it is $\approx 10\%$ too low at period ratios between the 3:2 and
2:1 resonances. This arises at least in part from the breakdown of the compact
approximation at these wide separations in the 2-planet MMR overlap criterion
(Appendix B). We also comment that it is not clear that our approximation of
simply summing MMR widths from both adjacent pairs of planets (Sec. 4.1)
should be valid in this equally spaced case. In this configuration, the MMRs
of the middle planet with its inner neighbor lie on top of those with its
outer neighbor, so one might expect summing their widths to be double-counting
their contribution to the MMR covering fraction. While this special case might
benefit from a more careful analysis, it is reassuring from Fig. 8 that our
simple approach still gives a close approximation to the chaotic boundary. We
have also checked a nearly equally spaced case where we space the outer pair
5% wider than the inner pair to offset the MMRs between adjacent planets, and
find comparable agreement to numerical results as in Fig. 8.

Finally, while Laskar & Petit (2017) and Petit et al. (2017) considered very
long timescales on which the secular eigenmodes might evolve through secular
chaos, it is instructive to compare what they would have predicted on shorter
timescales under the regular Laplace-Lagrange evolution we consider. In
particular, Laskar & Petit (2017) study conditions that would lead to orbit-
crossing, which in the compact limit reduces to $\tilde{e}=1$. In dashed blue,
we therefore plot the critical initial eccentricity that reaches a maximum
$\tilde{e}$ of unity under our assumption of Laplace-Lagrange evolution. We
see that in this case, this estimate is too high by a factor of $\approx 2-5$
compared to the numerical results.

Petit et al. (2017) also consider the possibility that the secular evolution
raises the eccentricities enough for first-order MMRs to overlap. In this
case, at separations beyond the 3-body MMR overlap region, we find that the
orbits always become orbit-crossing before the eccentricities get high enough
for first-order MMRs to overlap. This is not true in general (Petit et al.,
2017), but the condition for overlap of the set of all two-body MMRs (dashed
purple line in Fig. 8) always occurs at lower eccentricities than the
criterion for the overlap of only first-order MMRs.

##  6 Conclusion

We have investigated the boundary in parameter space for chaotic motion in
compact multiplanet systems on secular timescales (longer than the orbital
period by roughly the star-planet mass ratio). At the tightest orbital
separations (left of the red line in Fig. 8), the overlap of three-body MMRs
causes chaos even for initially circular orbits (Quillen, 2011; Petit et al.,
2020). At wider separations, we find that the chaotic boundary is instead set
by the overlap of 2-body MMRs at finite eccentricities. In this picture, a
better measure of how tightly the system is packed than the typically used
Hill separation is the dynamical spacing $S$ (Eq. 6) (Quillen, 2011; Hadden &
Lithwick, 2018; Petit et al., 2020).

We argued that a central reason why compact systems with three or more planets
are chaotic over a significantly broader range of parameter space than two-
planet configurations is a near-symmetry that suppresses changes in the widths
of MMRs in two-planet systems (Sec. 3). By contrast, configurations of three
or more planets undergo long-term secular eccentricity oscillations that cause
MMR widths to expand and contract, adiabatically modulating the chaotic
boundary (Sec. 4.3.2). We find that modeling this secular evolution to lowest
order in the eccentricities through classical Laplace-Lagrange theory is
adequate, and derive semi-analytic expressions for MMR overlap at the maximum
in the secular cycle that closely match the envelopes for chaotic motion
measured through N-body integrations.

A challenge for understanding the long-term stability of multiplanet systems
is the high dimensionality of the problem. We find that two approximations
significantly simplify the problem. First, we argue in Sec. 4.2 that the
widths of two-body MMRs decrease fast enough with orbital separation that it
is sufficient to consider adjacent trios of planets separately, and adopt the
most stringent condition among all trios. Second, we find that the coplanar
problem is a good approximation for mutual inclinations below a few degrees,
typical of transiting multiplanet systems (Sec. 4.4).

For a given trio of planets, this leaves four dimensionless parameters
determining the stability of the system: the two dynamical spacings $S_{12}$
and $S_{23}$ (Eq. 6), and the maximum normalized eccentricities (Eq. 4)
$\tilde{e}_{12}^{max}$ and $\tilde{e}_{23}^{max}$ attained during a secular
cycle. To visualize the results in the figures in the paper we reduce this
dimensionality by variously restricting the masses, spacings and orbital
periods, as well as by fixing the initialization of the secular eigenmodes.
For general cases, we provide routines for calculating the chaotic boundary as
part of the open-source SPOCK package (Tamayo et al., 2020) for community use.

We have checked that, with the exception of a minority of configurations
inside strong MMRs, systems beyond our chaotic boundary are short-lived. On
timescales much longer than the secular timescale, however, it is possible for
overlap of secular resonances (e.g., Lithwick & Wu, 2011) to drive secular
chaos that mixes the secular mode amplitudes and pushes the instability
boundary to lower eccentricities than our short-term chaotic boundary (black
line in Fig. 8). This question was pursued by Laskar & Petit (2017) and Petit
et al. (2017), who analyzed the worst-case scenario where secular chaos drives
all the system’s angular momentum deficit to a single pair of planets.
However, we have seen that the instability boundary they considered for the
short-term dynamics (blue line in Fig. 8) was not accurate; their
corresponding AMD criterion also requires modification. We derive an updated
AMD criterion in Appendix C, which shows that the AMD threshold below which
stability is guaranteed by Laskar & Petit (2017) is lowered by a
multiplicative factor that is exponential in the dynamical spacing. It
represents a significant correction for compact systems with $S\lesssim 3$.

An issue we have not resolved is the fate of planetary configurations on Gyr
timescales. Is secular chaos important in driving long-term instabilities in
observed compact multiplanet systems? Or does MMR overlap typically limit the
eccentricities to such small values on short timescales that secular
resonances are too narrow to overlap, and secular chaos is rare? The answer to
this question has important implications for the timing at which the final
architectures of compact exoplanet systems are set. We plan to address these
long-term dynamics in future work.

We make the scripts used to run the above numerical experiments and generate
their corresponding figures available at
https://github.com/dtamayo/chaospaper. We also provide routines for evaluating
our expressions and estimating whether input planetary configurations are
dynamically chaotic at https://github.com/dtamayo/spock as part of the open-
source SPOCK package (Tamayo et al., 2020).

##  7 Acknowledgements

We thank the referee, Antoine Petit, for a careful and insightful review. We
are also grateful to Sam Hadden and Alysa Obertas for illuminating
discussions. Support for DT was provided by NASA through the NASA Hubble
Fellowship grant HST-HF2-51423.001-A awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. This work was
supported in part by the Natural Sciences and Engineering Research Council of
Canada (NSERC), funding reference number RGPIN-2020-03885. This research was
made possible by the open-source projects REBOUND (Rein & Liu, 2012),
celmech888https://github.com/shadden/celmech., Jupyter (Kluyver et al., 2016),
iPython (Pérez & Granger, 2007), (Harris et al., 2020), and matplotlib
(Hunter, 2007; Droettboom et al., 2016).

##  Appendix A Chaos vs Instability Times

The MEGNO classifies whether orbits are chaotic over the timespan of an N-body
integration. In particular, it does not preclude chaotic behavior on longer
timescales. When it identifies a configuration as chaotic, there is also no
guarantee that the planets will undergo an instability within a reasonable
amount of time. Nevertheless, the MEGNO is a useful tool for rapidly surveying
parameter space, and for highlighting boundaries where the dynamical behavior
changes.

In Fig. 5, we chose to run for three secular cycles of the second-slowest
eigenmode, corresponding approximately to the precession period of the compact
planet pair under the perturbations from the exterior Jupiter, since it is on
this timescale that the MMRs evolve to their maximal widths. If instead we run
the for three secular cycles of the slowest eigenmode, corresponding
approximately to the precession period of the Jupiter perturber, we find that
the vast majority of the parameter space (even at low eccentricities) is
chaotic. This is true even when removing the Jupiter perturber, and
integrating for three secular cycles of the slowest eigenmode for the compact
pair of planets (which becomes long due to the near-symmetry in the two-planet
case discussed in Sec. 3). While we have not identified the source of this
chaos, it is not due to overlap of separate 2-body MMRs, and it is not clear
whether it can drive instabilities over long timescales, which is ultimately
what we are interested in.

In Fig. 9 we measure the instability times in integrations up to $10^{7}$
innermost planet orbits (defined as a close encounter within the sum of
planets’ individual Hill radii), showing that the regions below which 2-body
MMRs overlap are stable on these timescales, despite being chaotic on much
shorter timescales. Conversely, the vast majority of configurations above this
boundary go unstable. While it is unsatisfying to not have a more rigorous
delineation, we conclude that measuring chaos on timescales comparable to the
secular timescale on which the MMRs expand on contract secularly appropriately
captures the boundary for short-term instabilities, and defer the behavior
over much longer timescales to future work. We point out that this issue also
exists for previous criteria for the onset of chaos (e.g., Wisdom, 1980;
Hadden & Lithwick, 2018).

Figure 9:  Comparison of the MEGNO map in Fig. 5 to instability times measured
through integrations lasting up to $10^{7}$ orbits. Measuring the MEGNO over
integrations spanning the secular timescales on which the MMRs expand and
contract capture the appropriate chaotic boundary corresponding to short-term
instabilities.

##  Appendix B Eccentric Two-Planet Chaos Criterion

The value of the expressions for the MMR optical depth (Eq. 7) and critical
eccentricity for chaos (Eq. 8) derived by Hadden & Lithwick (2018) is that
they work over the full range of eccentricities. However, given that
eccentricities in compact multiplanet systems measured through TTVs are
typically small (e.g., Lithwick et al., 2012; Jontof-Hutter et al., 2014;
Hadden & Lithwick, 2017), it can also be useful to find simpler approximations
valid at low eccentricities.

In Fig. 10 we compare our low-eccentricity approximations to the general
expressions of Hadden & Lithwick (2018) for the case of two $1M_{\oplus}$
planets (left panel), two $10M_{\oplus}$ planets (middle panel) and two
$1M_{\rm{Jup}}$ planets (right panel). The period ratio ranges from the
critical separation (Eq. 1) at which first-order MMRs overlap at zero
eccentricity (Wisdom, 1980), out to approximately the period ratio where the
critical eccentricity for MMR overlap (Eq. 9) reaches $\tilde{e}_{12}=0.4$.
Integrations were run and the MEGNO was measured as described in the main
text, and analogous to our other numerical experiments, we integrate for three
periods of the fast Laplace-Lagrange eigenfrequency (since the nearly
conserved eigenmode has a much slower eigenfrequency, Sec. 3). We initialize
the eccentricities purely in the fast eigenmode as in Sec. 4.1. This step is
not necessary at close orbital separations due to the near-symmetry discussed
in Sec. 3, but makes a noticeable difference at wider separations between the
3:2 and 2:1 MMRs.

We find that the approximations in Eq. 9 fall within $10\%$ of those from
Hadden & Lithwick (2018) in the range $0.1<\tilde{e}_{12}<0.4$. We note that
the regime at lower eccentricities is not physically relevant, since first-
order MMRs typically overlap anyway below $\tilde{e}_{12}=0.1$ (see Fig. 10,
where the left-most period ratio has been chosen to match this criterion from
Wisdom 1980).

Figure 10:  Comparison of the critical eccentricity at which 2-body MMRs
overlap for two-planet systems from Hadden & Lithwick (2018) (Eq. 8, purple)
to our low-eccentricity approximation (Eq. 9, black). Both planets have the
same mass (given above each subplot). The low approximation agrees with Hadden
& Lithwick (2018) within 10% for $0.1<\tilde{e}_{12}\lesssim 0.4$.

Figure 11 shows an extension of the left two panels in Fig. 10 to higher
eccentricities and wider period ratios. The criterion of Hadden & Lithwick
(2018) matches the numerical results at all eccentricities, while the simpler
low-eccentricity approximation fails beyond $\tilde{e}_{12}\gtrsim 0.4$.

Figure 11:  Extension of the left two panels in Fig. 10 to higher
eccentricities and wider period ratios, showing the failure of the low-
eccentricity approximation beyond $\tilde{e}_{12}\approx 0.4$.

##  Appendix C An Updated AMD criterion

It is straightforward to generalize the AMD criterion of Petit et al. (2017)
to calculate the minimum AMD that allows for the overlap of higher order MMRs.
This was also explored by Petit (2019), but we derive an explicit expression
for the critical AMD below which the system is AMD stable. Indeed, our
construction of $\boldsymbol{e}_{12}$ from canonical variables through a
rotation makes it particularly simple. The total AMD $\equiv\mathcal{A}$ in
the coplanar case is given by (Eqs. 11, 13 and 14)

$\mathcal{A}=\sum_{i=1}^{N}\Gamma_{i}=\sum_{i=1}^{N}\frac{G_{i}^{2}}{2}$ (C1)

where $N$ is the number of planets. For any given pair of planets, which we
label with indices 1 and 2, we can apply the canonical rotation
$\mathcal{R}_{\psi}$ (Eq. 35), with additional unity entries along the
diagonals for the remaining planets, which remain unchanged. The rotation
preserves lengths, so we similarly have

$\sum_{i=1}^{N}\frac{G_{i}^{2}}{2}=\sum_{i=1}^{N}\frac{F_{i}^{2}}{2}=\mathcal{A}.$
(C2)

Finally, if we assume that all the AMD is transferred to $F_{1}$, which is
simply a scaled version of $\tilde{e}_{12}$ , we obtain the maximum estimate
of $\tilde{e}_{12}$ from Eq. 26,

$\tilde{e}_{12}^{max}=\sqrt{\frac{\Lambda_{1}+\Lambda_{2}}{\Lambda_{1}\Lambda_{2}}}\frac{\sqrt{2\mathcal{A}}}{e_{12}^{\rm{cross}}}$
(C3)

Combining this with Eq. 8, we obtain the critical AMD below which the pair
with indices 1 and 2 are AMD stable,

$\mathcal{A}_{\rm{crit}}=\frac{\Lambda_{1}\Lambda_{2}}{\Lambda_{1}+\Lambda_{2}}\frac{{(e_{12}^{\rm{cross}}})^{2}}{2}\exp\Bigg{[}-\Bigg{(}\frac{3.0}{S_{12}}\Bigg{)}^{4/3}\Bigg{]}$
(C4)

This procedure can be repeated for every adjacent pair of planets, and the
minimum $\mathcal{A}_{\rm{crit}}$ then yields the critical AMD for the whole
system.

Our expression is equivalent to that derived by Laskar & Petit (2017) in the
compact limit ($\alpha\approx 1$), except for the additional exponential term
in Eq. C4. Our corrected threshold can thus differ significantly from that of
Laskar & Petit (2017) for systems with $S\lesssim 3$ (see Fig. 2).

## References

  * Batygin & Morbidelli (2013) Batygin, K., & Morbidelli, A. 2013, Astronomical Journal, 145, 1, doi: 10.1088/0004-6256/145/1/1
  * Becker & Adams (2016) Becker, J. C., & Adams, F. C. 2016, Monthly Notices of the Royal Astronomical Society, 455, 2980
  * Chambers et al. (1996) Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261, doi: 10.1006/icar.1996.0019
  * Cincotta et al. (2003) Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D: Nonlinear Phenomena, 182, 151
  * Cranmer et al. (2021) Cranmer, M., Tamayo, D., Rein, H., et al. 2021, arXiv preprint arXiv:2101.04117
  * Dawson et al. (2016) Dawson, R. I., Lee, E. J., & Chiang, E. 2016, The Astrophysical Journal, 822, 54
  * Deck et al. (2013) Deck, K. M., Payne, M., & Holman, M. J. 2013, Astrophysical Journal, 774, 129, doi: 10.1088/0004-637X/774/2/129
  * Denham et al. (2019) Denham, P., Naoz, S., Hoang, B.-M., Stephan, A. P., & Farr, W. M. 2019, Monthly Notices of the Royal Astronomical Society, 482, 4146
  * Droettboom et al. (2016) Droettboom, M., Hunter, J., Caswell, T. A., et al. 2016, matplotlib: matplotlib v1.5.1, doi: 10.5281/zenodo.44579
  * Duncan et al. (1998) Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, Astronomical Journal, 116, 2067, doi: 10.1086/300541
  * Faber & Quillen (2007) Faber, P., & Quillen, A. C. 2007, Monthly Notices of the Royal Astronomical Society, 382, 1823
  * Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, The Astrophysical Journal, 790, 146
  * Faridani et al. (2021) Faridani, T., Naoz, S., Wei, L., & Farr, W. M. 2021, arXiv e-prints, arXiv:2107.07529.  https://arxiv.org/abs/2107.07529
  * Ferraz-Mello (2007) Ferraz-Mello, S. 2007, Canonical perturbation theories: degenerate systems and resonance, Vol. 345 (Springer Science & Business Media)
  * Funk et al. (2010) Funk, B., Wuchterl, G., Schwarz, R., Pilat-Lohinger, E., & Eggl, S. 2010, Astronomy & Astrophysics, 516, A82
  * Gladman (1993) Gladman, B. 1993, Icarus, 106, 247, doi: 10.1006/icar.1993.1169
  * Goldreich & Tremaine (1981) Goldreich, P., & Tremaine, S. 1981, Astrophysical Journal, 243, 1062
  * Gott (1993) Gott, J. R. 1993, Nature, 363, 315
  * Gratia & Lissauer (2019) Gratia, P., & Lissauer, J. J. 2019, arXiv preprint arXiv:1908.01117
  * Hadden (2019) Hadden, S. 2019, Astronomical Journal, 158, 238, doi: 10.3847/1538-3881/ab5287
  * Hadden & Lithwick (2016) Hadden, S., & Lithwick, Y. 2016, The Astrophysical Journal, 828, 44
  * Hadden & Lithwick (2017) —. 2017, The Astronomical Journal, 154, 5
  * Hadden & Lithwick (2018) —. 2018, The Astronomical Journal, 156, 95
  * Hansen & Murray (2012) Hansen, B. M., & Murray, N. 2012, The Astrophysical Journal, 751, 158
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hénon & Petit (1986) Hénon, M., & Petit, J.-M. 1986, Celestial Mechanics, 38, 67
  * Hernandez & Dehnen (2017) Hernandez, D. M., & Dehnen, W. 2017, Monthly Notices of the Royal Astronomical Society, 468, 2614
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jontof-Hutter et al. (2014) Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, The Astrophysical Journal, 785, 15
  * Kholshevnikov & Vassiliev (1999) Kholshevnikov, K. V., & Vassiliev, N. N. 1999, Celestial Mechanics and Dynamical Astronomy, 75, 75
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, Positioning and Power in Academic Publishing: Players, Agents and Agendas, 87
  * Lai & Pu (2017) Lai, D., & Pu, B. 2017, The Astronomical Journal, 153, 42
  * Laskar (2000) Laskar, J. 2000, Physical Review Letters, 84, 3240
  * Laskar & Gastineau (2009) Laskar, J., & Gastineau, M. 2009, Nature, 459, 817
  * Laskar & Petit (2017) Laskar, J., & Petit, A. 2017, Astronomy & Astrophysics, 605, A72
  * Lissauer & Gavino (2021) Lissauer, J. J., & Gavino, S. 2021, Icarus, 364, 114470
  * Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, The Astrophysical Journal Supplement Series, 197, 8
  * Lithwick & Wu (2011) Lithwick, Y., & Wu, Y. 2011, Astrophysical Journal, 739, 31, doi: 10.1088/0004-637X/739/1/31
  * Lithwick et al. (2012) Lithwick, Y., Xie, J., & Wu, Y. 2012, The Astrophysical Journal, 761, 122
  * Liu et al. (2015) Liu, B., Munoz, D. J., & Lai, D. 2015, Monthly Notices of the Royal Astronomical Society, 447, 747
  * Marchal & Bozis (1982) Marchal, C., & Bozis, G. 1982, Celestial Mechanics, 26, 311, doi: 10.1007/BF01230725
  * Marzari & Weidenschilling (2002) Marzari, F., & Weidenschilling, S. J. 2002, Icarus, 156, 570, doi: 10.1006/icar.2001.6786
  * Morbidelli & Moons (1993) Morbidelli, A., & Moons, M. 1993, Icarus, 102, 316
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge U. Press)
  * Obertas et al. (2017) Obertas, A., Van Laerhoven, C., & Tamayo, D. 2017, Icarus, 293, 52, doi: 10.1016/j.icarus.2017.04.010
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Petigura et al. (2017) Petigura, E. A., Howard, A. W., Marcy, G. W., et al. 2017, The Astronomical Journal, 154, 107
  * Petit (2019) Petit, A. 2019, PhD thesis, PSL Research University
  * Petit et al. (2017) Petit, A. C., Laskar, J., & Boué, G. 2017, Astronomy & Astrophysics, 607, A35
  * Petit et al. (2018) Petit, A. C., Laskar, J., & Boué, G. 2018, Astronomy & Astrophysics, 617, A93, doi: 10.1051/0004-6361/201833088
  * Petit et al. (2020) Petit, A. C., Pichierri, G., Davies, M. B., & Johansen, A. 2020, arXiv preprint arXiv:2006.14903
  * Pichierri et al. (2017) Pichierri, G., Morbidelli, A., & Lai, D. 2017, Astronomy & Astrophysics, 605, A23
  * Pu & Lai (2018) Pu, B., & Lai, D. 2018, Monthly Notices of the Royal Astronomical Society, 478, 197
  * Pu & Wu (2015) Pu, B., & Wu, Y. 2015, The Astrophysical Journal, 807, 44
  * Quarles et al. (2017) Quarles, B., Quintana, E. V., Lopez, E. D., Schlieder, J. E., & Barclay, T. 2017, arXiv preprint arXiv:1704.02261
  * Quillen (2011) Quillen, A. C. 2011, Monthly Notices of the Royal Astronomical Society, 418, 1043
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, Astronomy & Astrophysics, 537, A128
  * Rosenthal et al. (2019) Rosenthal, M., Jacobson-Galan, W., Nelson, B., et al. 2019, The Astronomical Journal, 158, 136
  * Sansottera & Libert (2019) Sansottera, M., & Libert, A.-S. 2019, Celestial Mechanics and Dynamical Astronomy, 131, 1
  * Sessin & Ferraz-Mello (1984) Sessin, W., & Ferraz-Mello, S. 1984, Celestial mechanics, 32, 307
  * Smith & Lissauer (2009) Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381
  * Steffen et al. (2013) Steffen, J. H., Fabrycky, D. C., Agol, E., et al. 2013, Monthly Notices of the Royal Astronomical Society, 428, 1077
  * Tamayo et al. (2021) Tamayo, D., Gilbertson, C., & Foreman-Mackey, D. 2021, Monthly Notices of the Royal Astronomical Society, 501, 4798
  * Tamayo et al. (2017) Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, The Astrophysical Journal Letters, 840, L19
  * Tamayo et al. (2015) Tamayo, D., Triaud, A. H., Menou, K., & Rein, H. 2015, The Astrophysical Journal, 805, 100
  * Tamayo et al. (2016) Tamayo, D., Silburt, A., Valencia, D., et al. 2016, The Astrophysical Journal Letters, 832, L22
  * Tamayo et al. (2020) Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proceedings of the National Academy of Sciences, 117, 18194
  * Tremaine (2015) Tremaine, S. 2015, Astrophysical Journal, 807, 157, doi: 10.1088/0004-637X/807/2/157
  * Udry et al. (2019) Udry, S., Dumusque, X., Lovis, C., et al. 2019, Astronomy & Astrophysics, 622, A37
  * Van Eylen & Albrecht (2015) Van Eylen, V., & Albrecht, S. 2015, The Astrophysical Journal, 808, 126
  * Van Eylen et al. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, The Astronomical Journal, 157, 61
  * Volk & Gladman (2015) Volk, K., & Gladman, B. 2015, The Astrophysical Journal Letters, 806, L26
  * Volk & Malhotra (2020) Volk, K., & Malhotra, R. 2020, arXiv preprint arXiv:2003.05040
  * Wang et al. (2018) Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, The Astronomical Journal, 156, 192
  * Wei et al. (2021) Wei, L., Naoz, S., Faridani, T., & Farr, W. M. 2021, arXiv preprint arXiv:2106.02276
  * Weiss & Marcy (2014) Weiss, L. M., & Marcy, G. W. 2014, The Astrophysical Journal Letters, 783, L6
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, The Astronomical Journal, 155, 48
  * Wisdom (1980) Wisdom, J. 1980, The Astronomical Journal, 85, 1122
  * Wisdom (1985) Wisdom, J. 1985, Icarus, 63, 272, doi: 10.1016/0019-1035(85)90011-9
  * Wisdom (1986) —. 1986, Celestial Mechanics, 38, 175, doi: 10.1007/BF01230429
  * Wisdom (2015) Wisdom, J. 2015, The Astronomical Journal, 150, 127
  * Wu & Lithwick (2011) Wu, Y., & Lithwick, Y. 2011, Astrophysical Journal, 735, 109, doi: 10.1088/0004-637X/735/2/109
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Sciences, 113, 11431
  * Yoshinaga et al. (1999) Yoshinaga, K., Kokubo, E., & Makino, J. 1999, Icarus, 139, 328
  * Zhou et al. (2007) Zhou, J.-L., Lin, D. N., & Sun, Y.-S. 2007, The Astrophysical Journal, 666, 423