## Eccentricities and the Stability of Closely-Spaced Five-Planet Systems

Pierre Gratia a, ∗ , Jack J. Lissauer 1

a CIERA and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208

b Space Science &amp; Astrobiology Division, MS 245-3, NASA Ames Research Center, Moffett Field, CA 94035

## Abstract

Observations of exoplanets have revealed that systems with planets on closely-spaced orbits are common, which motivates the question 'How closely can planets orbit to one another and still be dynamically-stable for very long times?'. To address this question, we investigate the stability of idealized planetary systems consisting of five planets, each equal in mass to the Earth, orbiting a one solar mass star. All planets orbit in the same plane and in the same direction, and the planets are uniformly spaced in units of mutual Hill Sphere radii. Most of the systems that we integrate begin with one or more planets on eccentric orbits, with eccentricities e as large as e = 0 05 being considered. . For a given initial orbital separation, larger initial eccentricity of a single planet generally leads to shorter system lifetime, with little systematic dependence of which planet is initially on an eccentric orbit. The approximate trend of instability times increasing exponentially with initial orbital separation of the planets found previously for planets with initially circular orbits is also present for systems with initially eccentric orbits. Mean motion resonances also tend to destabilize these systems, although the reductions in system lifetimes are not as large as for initially circular orbits. Systems with all planets having initial e = 0 05 and aligned periapse angles typically survive far longer than systems with . the same spacing in initial semi-major axis and one planet with e = 0 05, but they . have slightly shorter lifetimes than those with planets initially on circular orbits.

Keywords: Extrasolar planets, Celestial mechanics, Dynamical simulations

∗ Corresponding author

Email addresses: pierregratia@gmail.edu (Pierre Gratia), jack.lissauer@nasa.gov (Jack J. Lissauer)

## 1. Introduction

NASA's Kepler mission has discovered hundreds of multiple planet systems (Fabrycky et al. (2014); Rowe et al. (2014)). Many of these multi-planet systems are quite tightly packed, where adjacent planets orbit much closer to one another than do those in the Solar System (Lissauer et al. (2011b)). For instance, the Kepler-11 system harbors six known planets, five of which have orbital periods between 10 and 47 days (Lissauer et al. (2011a)). Studying the orbital stability of such tightly-packed planetary systems provides clues to how they form and how long they survive. The present study can be seen as an extension of both Obertas et al. (2017) (henceforth OVT17), who investigated the stability of closely-spaced five-planet systems with the same masses and relative spacings as we use herein and all planets on initially circular orbits; and of Smith and Lissauer (2009), who simulated far fewer simulations, but include a few systems with both initial eccentricity and inclination. Here, we expand on those findings and look at systems for which one planet starts out with a nonzero initial eccentricity, at a 10 times higher in initial orbital resolution than Smith and Lissauer (2009). We also consider various systems in which all planets begin on eccentric orbits with e = 0 05. .

The stability of systems of two planets has been characterized analytically; see Hill (1878) and Gladman (1993). In contrast, few constraints on the stability of more than two planet systems can be derived analytically. Such an investigation therefore requires extensive numerical integrations, for example Chambers et al. (1996); Marzari and Weidenschilling (2002); Smith and Lissauer (2009); Pu and Wu (2015); Morrison and Kratter (2016); Obertas et al. (2017).

The structure of the paper is as follows. We start by describing the setup of our simulations in Section 2. In Section 3, the focus is on circular systems and the well-known log-linear relationship between initial planetary separation and time to close encounter, which we reproduce with almost identical regression parameters. We then move on to the core part of this work. In Section 4, we study systems that start out with the innermost planet eccentric. We consider several different eccentricities and compare the systems' lifetimes. In Section 5, we perform similar integrations, giving the outermost planet nonzero eccentricity, and in Section 6, we conduct a more limited study with one of the intermediate planets given nonzero eccentricity. Section 7 considers systems with all planets initially having eccentricity of 0.05; we study systems in which all planets initially have the same longitude of periapse, as well as systems in which the longitudes of periapse are chosen randomly. We compare the coefficients of the log-linear lifetime relationships for many batches of runs in Section 8. Our conclusions are summarized in Section 9.

## 2. Methodology

Our systems consist of five identical, 1 M ⊕ (Earth-mass) planets, orbiting a 1 M glyph[circledot] (solar mass) star in the same direction, and initially separated by the same multiple of their mutual Hill radii (equation 2). We restrict ourselves to a two dimensional study by setting all inclinations to zero. We investigate the effects of nonzero initial eccentricities on these five-planet systems.

In what follows, time is measured in units of the initial orbital period of the innermost planet, i.e., Earth-years if the innermost planet begins at 1 AU. Our simulations were performed with REBOUND, the N-body simulation code by Rein and Liu (2012). In particular, we use the symplectic Wisdom-Holman integrator (WHFast) provided by REBOUND (Rein and Tamayo (2015)). WHFast is fast and unbiased, but ill-suited for resolving collisions and close encounters. Thus, we stop a simulation as soon as a close encounter occurs, and compare lifetimes computed in this manner. Rice et al. (2018) investigated similar but shorter-lived systems, and followed their evolution beyond their close encounter time (for systems up to and including initial separations of eight mutual Hill radii), while we focus more closely on long-lived systems and stopped our simulations at the time of close encounter, including systems separated by initial distances well beyond eight mutual Hill radii.

## 2.1. Mutual Hill radius

Before describing the initial conditions, we recall the definition of the Hill radius, which is our unit of distance between initial planetary orbits. The Hill radius of an astronomical object, here a planet, orbiting a much more massive object, here a star, is the radius of the zone in which the planet's gravitational pull is the primary force affecting motion relative to the planet. For small eccentricities, it can be approximated as

$$R _ { H } \approx a \left ( \frac { m } { 3 M _ { * } } \right ) ^ { 1 / 3 },$$

where M glyph[star] and m are the masses of the star and planet, respectively, and a is the semi-major axis. The Hill radius is then proportional to the semi-major axis. In this paper, following convention for this type of simulations, we use a slightly modified version: the mutual Hill radius between planet i + 1 and its immediate neighbor orbiting closer to the host star, planet i :

$$R _ { H _ { i, i + 1 } } \approx \left ( \frac { m _ { i + 1 } + m _ { i } } { 3 M _ { * } } \right ) ^ { 1 / 3 } \left ( \frac { a _ { i } + a _ { i + 1 } } { 2 } \right ).$$

Here the denominator of the first ratio continues to be the host star's mass, as in OVT17, but Quarles and Lissauer (2018) use the more accurate formula M glyph[star] +( i -1) M ⊕ . The initial separations between the planets are expressed in terms of their mutual Hill radii, with the innermost planet placed at 1 AU.

## 2.2. Initial positions and stopping conditions

We start by stating the initial conditions and describing the parameters used for our simulations.

## 2.2.1. Initial positions

The initial semi-major axes of the five planets are always chosen to be separated by the same number (not necessarily integer), β , times their mutual Hill radii (we follow the notation of Smith and Lissauer (2009), while OVT17 use ∆ instead of β ):

$$a _ { i + 1 } - a _ { i } = \beta R _ { H _ { i, i + 1 } }.$$

From Equations (2) and (3), we have:

$$a _ { i + 1 } = a _ { i } \left [ 1 + \frac { \beta } { 2 } \left ( \frac { m _ { i } + m _ { i + 1 } } { 3 M _ { * } } \right ) ^ { 1 / 3 } \right ] \left [ 1 - \frac { \beta } { 2 } \left ( \frac { m _ { i } + m _ { i + 1 } } { 3 M _ { * } } \right ) ^ { 1 / 3 } \right ] ^ { - 1 }.$$

We perform batches of integrations in which we increment the initial β by 0.01 for each new simulation, while keeping everything else fixed at some specified nonzero initial value. Systems of two planets on initially circular orbits are stable if β &gt; 2 √ 3 ≈ 3 4641, so we begin our numerical studies at . β = 3 47. Although this critical . two-planet separation criterion is strictly valid only for circular systems (Petit et al., 2018), our initial eccentricities are small enough that this starting value for β is justified - not least because we are mainly interested in the long-lived systems.

The true longitude is defined as

$$\theta \equiv \Omega + \omega + \nu,$$

where Ω is the longitude of ascending node, ω the argument of pericenter, and ν the true anomaly. Since we are working in a co-planar setting, the longitude of ascending node Ω is zero, and for convenience we also set the argument of peripasis ω = 0. This leaves us with the initial true anomaly ν , which is the angle from the pericenter to the planet's position, and the only nonzero parameter, so the true longitude is θ = ν . In this work, we choose it as θ i = 2 πiλ radians, where i is the i -th planet starting from the innermost one with i = 1, and λ = 1+ √ 5 ≈ 1 618 is the golden ratio. . This choice avoids any pair of planets starting near conjunction.

Finally, our time step has been chosen to be 0.049281314 years ≈ 18 days.

## 2.2.2. Stopping conditions

We stopped any given simulation if a close encounter occurred, which we define as the distance between any two objects at the end of any time step being below 0 01 AU. We did not encounter a case where a planet has been ejected. . We also set two more stopping conditions: one, due to limited computational resources, for

any given simulation, namely a maximum run time of 10 10 orbits of the innermost planet's initial orbit (i.e., Earth-years, years hereafter); and one for a given batch of simulations running over a β -range, namely five of them exceeding 2 × 10 9 ≈ 10 9 3 . years. For systems that did not have a close encounter within ten billion years, we used ten billion years as the lifetime in the statistics, which produces a smaller bias in our fits than ignoring those runs would.

## 3. Initially Circular Orbits

For the circular case, we ran our simulations in the region β ∈ [3 47 . , 8 73], . the upper bound representing the initial separation of the fifth system that exceeded two billion years before a close encounter occurred. System lifetimes are shown as a function of β in Figure 1. The general patterns that we find, with lifetimes increasing roughly exponentially with initial separation apart from dips near mean motion resonances, are similar to those noted by SL09 and OVT17. Our results show somewhat less scatter for similar values of β than OVT17, presumably because they performed far more simulations and randomly selected initial planetary longitudes (with the latter choice being especially important for small values of β ).

In the two subsequent figures, we also added the exponential fits for a selection of β ranges. The blue line in Figure 1 corresponds to the best-fitting exponential dependence of system lifetime as a function of initial orbital separation over the entire interval in β that we integrated; the orange line shows the fit for 5 32 . ≤ β ≤ 8 73, . the interval for which we have good stability lifetime values for all sets of runs with the outer four planets beginning on circular orbits; the green line shows the best fit including only systems with β ≤ 8 4; its slope can be compared to previous work by . Smith and Lissauer (2009) and OVT17 (see Table 1). The lines differ only slightly over these β -ranges, with the largest slope (in orange) belonging to the highest considered lower bound, where β ∈ [5 32 . , 8 73]. .

Figure 1: Close encounter times for systems whose planets all start out with zero eccentricity as a function of the separation in initial semi-major axis measured in mutual Hill radii β . On the vertical axis, t e stands for 'time to encounter'.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-b20a53b0-4426-474d-bc71-540adfe79d6a.jpg)

The exponential increase in system lifetimes with separation is generally quantified using a log-linear relationship between the close encounter time, t e , and initial separation of the planetary orbital semi-major axes in units of the mutual Hill radius β :

$$\log \left ( \frac { t _ { e } } { t _ { 0 } } \right ) = b \beta + c,$$

where t 0 is the initial orbital period of the innermost planet (one year), and b and c are constants. As noted in OVT17, however, this exponential fit does not capture all of the variations in close encounter times, which are up to two orders of magnitude off from the regression line. Note also that attempts have been made to fit a power law between initial separation and close encounter time (Zhou et al. (2007)).

We compare our coefficients with those found for the nonzero eccentricity cases in Section 4. A full tabulation of all computed coefficients and their standard errors σ can be found in Table 3. There, however, we considered new parameters b ′ and c ′ , computed from a new variable β ′ defined by

$$\beta ^ { \prime } \equiv \beta - 2 \sqrt { 3 },$$

following Quarles and Lissauer (2018):

$$\log \left ( \frac { t _ { e } } { t _ { 0 } } \right ) = b ^ { \prime } \beta ^ { \prime } + c ^ { \prime }.$$

Table 1: Comparison of slopes and intercepts of exponential fits to lifetimes of initially circular five planet systems in our simulations with those of previous studies.

| Reference                 | b                 | c                   | Range               |
|---------------------------|-------------------|---------------------|---------------------|
| Smith and Lissauer (2009) | 1.012             | - 1 . 686           | 3 . 4 ≤ β ≤ 8 . 4 √ |
| Obertas et al. (2017)     | 0.951             | - 1 . 202           | 2 3 < β < 8 . 4 √   |
| Obertas et al. (2017)     | 1.086             | - 1 . 881           | 2 3 < β < 10 √      |
| This work                 | 0 . 972 ± 0 . 018 | - 1 . 333 ± 0 . 107 | 2 3 < β ≤ 8 . 4 √   |
| This work                 | 1 . 080 ± 0 . 019 | - 1 . 895 ± 0 . 119 | 2 3 < β ≤ 8 . 73    |

## 4. Innermost planet eccentric

We now turn to the investigation of systems with initial conditions analogous to those considered in Section 3, with the exception that the innermost planet starts out with a nonzero eccentricity. In our first batch of simulations, we set e 1 = 0 01. . The resulting lifetimes are plotted in Figure 2. Predictably, compared to the circular case, the close encounter times of those systems tend to be shorter for given β ; hence we needed to go to higher values of β to find five systems that survived for &gt; 2 × 10 9 years.

Figure 2: Close encounter times for systems whose innermost planet starts out with eccentricity e 1 = 0 01; the upward-pointing triangle in the upper right represents a system that survived for . the entire 10 10 year interval simulated. In blue the corresponding regression line; in orange, the regression line for the β interval for which we have values for all batches with the inner planet eccentric and the outer four planets beginning on circular orbits, which is [5 32 . , 10 05]. . All slopes and intercept values can be found in Table 3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-b4ce405a-1179-4706-a24a-090bba8125e4.jpg)

As the eccentricity increases, we expect the lifetimes to decrease for fixed β . This is indeed what happens, and Figure 3 shows the behavior of the highest eccentricity inner planet systems we simulated, e 1 = 0 05. . The persistence of a log-linear relationship between initial separation and close encounter times for eccentricities at least up to 0 05 can be seen in Figure 4, with various slopes. . Several systems (one for e 1 = 0 02, eight for . e 1 = 0 03, and twenty-two for . e 1 = 0 05), suffered a close . approach extremely rapidly, with lifetimes less than ten years. Those encounters are between planets 1 and 2 and occur during their first synodic period: given our choices of initial longitudes and periapse angle of the inner planet, the first passage occurs when the inner planet is near apoapse. Those excessively short lifetimes are thus direct consequences of our initial longitudes.

Figure 3: Close encounter times for systems whose innermost planet starts out with eccentricity e 1 = 0 05. . In blue the corresponding regression line, starting just after the widest separation system that has a collision during the first synodic period; in orange, the regression line for the same interval in β as the orange line in Figure 2

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-7eeb3b27-a8f5-4a4e-90c8-39ac63f0b591.jpg)

Figure 4: Close encounter times for initial e 1 = 0 0 01 0 02 0 03, and 0 05, respectively. , . , . , . .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-5ac12592-6994-473d-b10a-515ebe008260.jpg)

We conclude that for any given value of e 1 simulated, system lifetime tends to roughly follow the log-linear relationship with β of the form given in Eq. 6, although the coefficients differ from those of initially circular orbits.

We now quantify the differences in close encounter times by computing their average over a range of three Hill radii, ending with the last point available for the circular case (i.e., the fifth initially circular system exceeding two billion years of stability), which amounts to the range β ∈ [5 74 . , 8 73]. .

Thus in Figure 5, each point represents the log-difference in close encounter time between one of the four eccentricity cases studied and the initially circular system. Not surprisingly, the higher the eccentricity, the larger the close encounter time differences. The four averages over the whole range are recorded in Table 2.

We omit similar plots for other planets eccentric - the ratios follow the same pattern for comparisons of different eccentricities of the same planet. However, an interesting pattern is found for the comparison of the average close encounter time between different planets starting out with the same eccentricity (Figure 10).

Figure 5: Close encounter time differences between simulations with all planets starting on circular orbits from those with e 1 = 0 01 . , e 1 = 0 02, and . e 1 = 0 05. . In the legend, [ e 1 = 0] is a time unit and means log t e [ e 1 = 0]; and so on. The average time difference between the circular system and the corresponding (same β ) eccentric system is given in Table 2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-29129ee7-5c23-4df2-a375-ba9cb83fca74.jpg)

Table 2: The 'average multiplicative factor' in close encounter time difference for each batch of systems compared to the circular systems. We compute this factor by exponentiating the average of the differences in the logarithms of the encounter time of simulations with the same value of β , within a specified range in β . The range consists of the last 300 overlapping points β ∈ [5 74 . , 8 73]. . The batch of simulations labeled e align = 0 05 has all planets on eccentric and aligned orbits; see . Section 7 for details.

|   circular | e 1 = 0 . 01   |   e 1 = 0 . 02 |   e 1 = 0 . 03 |   e 1 = 0 . 05 |   e align = 0 . 05 |
|------------|----------------|----------------|----------------|----------------|--------------------|
|          1 | 0 . 543        |          0.237 |          0.091 |          0.013 |              0.331 |

## 5. Outermost planet eccentric

In this section, we present simulations analogous to those in Section 4 except that the initial eccentricity is given to the outer planet. We expect results qualitatively

similar to those in which the inner planet is eccentric, but quantitative differences may exist.

We introduce here the Angular Momentum Deficit (AMD). For the planar problem that we are studying herein, the AMD of a planet is defined as the difference between the angular momentum it would have on a circular orbit of radius equal to its semi-major axis and its actual orbital angular momentum:

$$\text{AMD} = m \sqrt { G M a } ( 1 - \sqrt { 1 - e ^ { 2 } } ).$$

The total AMD of a system is not changed by secular interactions among the planets, and thus the close encounter times of our systems may be strongly correlated with their initial AMD. (See also Laskar (1997), and Laskar and Petit (2017) for an analytically derived AMD-stability criterion.) Comparing the AMD of a system with inner planet eccentric with one having planets at the same semi-major axes but the outer planet having the same eccentricity, we see that because the outer planet has larger a , the system with the outer planet eccentric has larger AMD.

What role does angular momentum deficit play in the destabilization of the closely-packed planetary systems that we are investigating? We look to answer this question in the following way. By giving nonzero eccentricity not to the innermost planet, but to a different one instead, we can calculate the expected outcome if AMD is at play. In particular, if the initial value of the AMD is key, then there should be no systematic difference in lifetimes if the nonzero initial eccentricity of the new planet is adjusted relative to its respective initial semi-major axis such that it matches the initial AMD of the innermost planet, whereas systems at the same β with the outermost planet eccentric would on average be shorter-lived than those with the innermost planet having the same value of eccentricity.

glyph[negationslash]

We performed six batches of simulations with only the outer planet starting on an eccentric orbit. Three of these used fixed values of e 5 . The other three adjusted e 5 to match the AMD of the corresponding system where the innermost planet starts out with nonzero initial eccentricity e 1 . For each matched AMD simulation (i.e., each value of β ), the initial eccentricity, e 5 , is adjusted such that the initial angular momentum deficit of that system is the same as that for the e 1 = 0 system:

$$m _ { 1 } \sqrt { G M _ { \odot } a _ { 1 } } \Big ( 1 - \sqrt { 1 - e _ { 1 } ^ { 2 } } \Big ) = m _ { 5 } \sqrt { G M _ { \odot } a _ { 5 } } \Big ( 1 - \sqrt { 1 - e _ { 5 } ^ { 2 } } \Big ). \quad \quad ( 1 0 )$$

This equation can be solved for e 5 . Taking into account equal planetary masses m 1 = m 5 , we have

$$e _ { 5 } & = \left [ 1 - \left ( 1 - \sqrt { \frac { a _ { 1 } } { a _ { 5 } } } \left ( 1 - \sqrt { 1 - e _ { 1 } ^ { 2 } } \right ) \right ) ^ { 2 } \right ] ^ { 1 / 2 }.$$

We use Equation (4) to get the semi-major axes ratio, and for simplicity define

$$X \equiv \frac { 1 } { 2 } \left ( \frac { 2 M _ { \oplus } } { 3 M _ { \odot } } \right ) ^ { 1 / 3 } \approx 0. 0 0 6 3.$$

We recursively apply Equation (4) giving us a 5 in terms of a 1 :

$$a _ { 5 } = \left ( \frac { 1 + \beta X } { 1 - \beta X } \right ) ^ { 4 } a _ { 1 }.$$

Plugging this back into Equation (11), we finally have

$$e _ { 5 } = \left [ 1 - \left ( 1 - \left ( \frac { 1 + \beta X } { 1 - \beta X } \right ) ^ { - 2 } \left ( 1 - \sqrt { 1 - e _ { 1 } ^ { 2 } } \right ) \right ) ^ { 2 } \right ] ^ { 1 / 2 }.$$

The choice of initial e 5 thus depends on both the initial separation β and the corresponding system's initial eccentricity e 1 of the innermost planet.

Figure 6: Close encounter times for initial e 5 = 0 0 01 0 02 0 05, and the AMD-adjusted , . , . , . e 5 . The AMD-adjusted systems tend to be slightly shorter-lived than those starting out at e 5 = 0 05 for all . β , which is to be expected, since the adjustment puts the initial e 5 eccentricity on an interval going from 0 048 for . β = 3 47 down to 0 043 for . . β = 12 32, thus always remaining smaller than but close . to 0 05 eccentricity. . However, the interesting part is to compare those systems with those whose innermost planet begins at e 1 = 0 05, since their initial AMD is identical by construction, while . the initial eccentricity is in different planets. This figure is the counterpart to Figure 4, where we plotted the lifetimes of runs with different eccentricities for the innermost planet.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-b5ff92c8-8720-4a49-908a-4e3cda5545a0.jpg)

Table 3 presents coefficients of the exponential fits to all six batches of systems with just the outer planet eccentric. Figure 6 shows the close encounter times, again as a function of β , of all three batches with constant values of the outermost planet's initial eccentricity, e 5 , as well as the batch with e 5 chosen to match the AMD of the run with e 1 = 0 05. . General trends are, as expected, similar to those seen for analogous batches with the inner planet eccentric.

Figure 7: Close encounter times of systems with initial e 1 = 0 05 . , e 5 = 0 05, and . e 5 = f amd ( e 1 = 0 05 . , β ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-7c6d42ba-f731-4259-82f4-9627ade72820.jpg)

Figure 7 compares lifetimes for batches of systems with either the inner or outer planet beginning with eccentricity e ≈ 0 05, and Figure 10 shows difference in close . encounter time between innermost and outermost planet for both the standard e 5 batch and the AMD adjusted batch, complementing the comparison of different eccentricities for the same planet shown in Figure 5. There, we focused on the change in close encounter times for different initial eccentricities of the same (innermost) planet, whereas here our focus is on the change in close encounter times for different planets with the same initial eccentricity (we choose the largest eccentricities, e ≈ 0 05). .

Figure 10 shows that the systems with e 5 ≈ 0 05 . and those with e 1 = 0 05 . have similar lifetimes. Only the 300 most widely-spaced systems for which we have data, β ∈ [9 33 . , 12 32] . are shown here, because the long-lived systems with these separations are more relevant to observed exoplanet systems than the more closelyspaced systems that we modeled. The e 5 = 0 05 systems on average have a somewhat .

shorter close encounter time than those with e 1 = 0 05, but the difference is small, . with the average multiplicative factor being 0.881. On the other hand, the red points for the AMD-adjusted batch lie significantly below the line of zero difference: The AMD-adjusted systems live much longer than the e 1 = 0 05 systems, with an average . ratio of 2.393. This is mostly due to the smaller initial eccentricities of those systems: they range from e 5 = 0 0444 at . β = 9 33 to . e 5 = 0 0428 at . β = 12 32. .

The slopes, b ′ , and intercepts, c ′ , of exponential fits to all batches of runs are given in Table 3, using a separation shifted by the critical separation for Hill stability of two-planet systems (eq. 7. Table 3 also gives the values of σ b ′ and σ c ′ , the standard errors of our estimates of the slopes and intercepts). For most batches, we considered multiple intervals in β . We begin by comparing our numbers with circular initial orbits (lines one and two). We then record the data for all batches, starting from the first point beyond which there are no more invalid points, until the fifth system within the batch reaching two billion years of lifetime. Next, we compare innermost planet eccentricities including e 1 = 0 05, before comparing all batches . except any single e i = 0 05 (but including . e align = 0 05). . Next are the batches with one planet having eccentricity ≈ 0 05. . Finally, fit coefficients to batches of systems with all planets eccentric and randomized periapses within a given range are given.

glyph[negationslash]

Table 3: Slopes and intercepts and their standard deviations for log-linear fits to lifetimes of simulated batches of systems.

| Eccentricity                          | Range                           | b ′     | σ b ′   | c ′    | σ c ′   |
|---------------------------------------|---------------------------------|---------|---------|--------|---------|
| 0                                     | 3 . 47 ≤ β ≤ 8 . 73             | 1.080   | 0.019   | 1.85   | 0.06    |
| e 1 = 0 . 01                          | 3 . 47 ≤ β ≤ 8 . 73             | 0.999   | 0.016   | 1.87   | 0.05    |
| last                                  | invalid → fifth beyond 2 × 10 9 |         |         |        |         |
| e 1 = 0 . 01                          | 3 . 47 ≤ β ≤ 10 . 05            | 1.030   | 0.014   | 1.82   | 0.05    |
| e 1 = 0 . 02                          | 3 . 48 ≤ β ≤ 10 . 34            | 0.984   | 0.010   | 1.61   | 0.04    |
| e 1 = 0 . 03                          | 4 . 19 ≤ β ≤ 11 . 21            | 0.867   | 0.011   | 1.58   | 0.05    |
| e 1 = 0 . 05                          | 5 . 32 ≤ β ≤ 12 . 61            | 0 . 785 | 0.010   | 1 . 02 | 0.06    |
| e 2 = 0 . 05                          | 5 . 49 ≤ β ≤ 12 . 99            | 0 . 827 | 0.010   | 0 . 52 | 0.06    |
| e 3 = 0 . 01                          | 3 . 47 ≤ β ≤ 10 . 09            | 1.017   | 0.012   | 1 . 82 | 0.05    |
| e 3 = 0 . 02                          | 4 . 00 ≤ β ≤ 11 . 09            | 0.921   | 0.012   | 1 . 64 | 0.06    |
| e 3 = 0 . 05                          | 5 . 07 ≤ β ≤ 12 . 99            | 0.810   | 0.009   | 0 . 56 | 0.06    |
| e 5 = 0 . 01                          | 3 . 47 ≤ β ≤ 9 . 98             | 1.016   | 0.014   | 1 . 84 | 0.05    |
| e 5 = 0 . 02                          | 3 . 47 ≤ β ≤ 10 . 80            | 0.911   | 0.012   | 1 . 79 | 0.05    |
| e 5 = 0 . 05                          | 4 . 86 ≤ β ≤ 12 . 38            | 0.759   | 0.010   | 1 . 13 | 0.06    |
| e 5 = e amd 1 | 0 . 01                | 3 . 47 ≤ β ≤ 10 . 05            | 0.981   | 0.013   | 1 . 93 | 0.05    |
| e 5 = e amd 1 | 0 . 02                | 3 . 47 ≤ β ≤ 10 . 30            | 0.983   | 0.011   | 1 . 68 | 0.05    |
| e 5 = e amd 1 | 0 . 05                | 4 . 85 ≤ β ≤ 12 . 32            | 0.826   | 0.011   | 1 . 02 | 0.06    |
| e align = 0 . 05                      | 3 . 47 ≤ β ≤ 10 . 08            | 0.837   | 0.016   | 2 . 18 | 0.06    |
| e 1 = 0:                              | overlapping range:              |         |         |        |         |
| e 1 = 0 . 01                          | 5 . 32 ≤ β ≤ 10 . 05            | 1.056   | 0.025   | 1.69   | 0.11    |
| e 1 = 0 . 02                          | 5 . 32 ≤ β ≤ 10 . 05            | 0.973   | 0.020   | 1.60   | 0.09    |
| e 1 = 0 . 03                          | 5 . 32 ≤ β ≤ 10 . 05            | 0.907   | 0.017   | 1.44   | 0.08    |
| e 1 = 0 . 05                          | 5 . 32 ≤ β ≤ 10 . 05            | 0 . 767 | 0.016   | 1 . 11 | 0.07    |
| all except e ≈ 0 . 05:                | overlapping range:              |         |         |        |         |
| e 1 = 0 . 01                          | 4 . 21 ≤ β ≤ 10 . 05            | 1.033   | 0.017   | 1.80   | 0.07    |
| e 1 = 0 . 02                          | 4 . 21 ≤ β ≤ 10 . 05            | 0.969   | 0.014   | 1.63   | 0.05    |
| e 1 = 0 . 03                          | 4 . 21 ≤ β ≤ 10 . 05            | 0.893   | 0.012   | 1.50   | 0.05    |
| e 3 = 0 . 01                          | 4 . 21 ≤ β ≤ 10 . 05            | 1.005   | 0.016   | 1 . 87 | 0.06    |
| e 3 = 0 . 02                          | 4 . 21 ≤ β ≤ 10 . 05            | 0.954   | 0.013   | 1 . 56 | 0.05    |
| e 5 = 0 . 01                          | 4 . 21 ≤ β ≤ 10 . 05            | 1.009   | 0.017   | 1 . 88 | 0.07    |
| e 5 = 0 . 02                          | 4 . 21 ≤ β ≤ 10 . 05            | 0.951   | 0.013   | 1 . 69 | 0.05    |
| e 5 = e amd 1 | 0 . 01                | 4 . 21 ≤ β ≤ 10 . 05            | 0.963   | 0.017   | 2 . 02 | 0.07    |
| e 5 = e amd 1 | 0 . 02                | 4 . 21 ≤ β ≤ 10 . 05            | 0.961   | 0.015   | 1 . 72 | 0.06    |
| e align = 0 . 05                      | 4 . 21 ≤ β ≤ 10 . 05            | 0.790   | 0.020   | 2 . 37 | 0.08    |
| e ≈ 0 . 05:                           | overlapping range:              |         |         |        |         |
| e 1 = 0 . 05                          | 5 . 49 ≤ β ≤ 12 . 32            | 0 . 755 | 0.011   | 1 . 16 | 0.06    |
| e 2 = 0 . 05                          | 5 . 49 ≤ β ≤ 12 . 32            | 0 . 778 | 0.011   | 0 . 73 | 0.06    |
| e 3 = 0 . 05                          | 5 . 49 ≤ β ≤ 12 . 32            | 0.765   | 0.011   | 0 . 76 | 0.07    |
| e 5 = 0 . 05                          | 5 . 49 ≤ β ≤ 12 . 32            | 0.742   | 0.011   | 1 . 22 | 0.07    |
| e 5 = e amd 1 | 0 . 05                | 5 . 49 ≤ β ≤ 12 . 32            | 0.829   | 0.012   | 1 . 00 | 0.07    |
| all e ≈ 0 . 05, randomized periapses: | range:                          |         |         |        |         |
| e all | 45 ◦ = 0 . 05                 | 3 . 47 ≤ β ≤ 11 . 24            | 0.749   | 0.012   | 1 . 81 | 0.06    |
| e all | 90 ◦ = 0 . 05                 | 3 . 47 ≤ β ≤ 12 . 60            | 0.628   | 0.009   | 1 . 48 | 0.05    |
| e all | 135 ◦ = 0 . 05                | 3 . 47 ≤ β ≤ 13 . 10            | 0.593   | 0.009   | 1 . 04 | 0.05    |
| e all | 180 ◦ = 0 . 05                | 3 . 47 ≤ β ≤ 14 . 19            | 0.532   | 0.008   | 0 . 95 | 0.05    |
| e all | 360 ◦ = 0 . 05                | 3 . 47 ≤ β ≤ 15 . 18            | 0.479   | 0.006   | 0 . 69 | 0.04    |

## 6. One intermediate planet eccentric

We close the main part of this investigation by considering systems in which one of the intermediate planets is initially eccentric. For most of these runs we choose the middle (third) planet initially eccentric, and also included a run with e 2 = 0 05. . Figure 8 confirms the log-linear relationship for those planetary systems.

We compare the 0 05 systems in Figure 9. We find slightly lower close encounter . times for eccentric intermediate planets, as shown in Figure 10. We thus confirm that a log-linear relationship continues to hold for intermediate planets eccentric, and find the close encounter times decrease on average by roughly a factor of two when the initially eccentric planet has both interior and exterior neighbors. The slopes ( b ′ ) and intercepts ( c ′ ) fit for the log-linear relationship in the overlap range β ∈ [5 49 . , 12 32] . are similar for all four runs, with the systems with the initially eccentric planet in the middle having slightly steeper slopes but somewhat smaller intercept values than those with the inner or outer planet eccentric. The e 5 = e amd 1 | 0 05 . batch has a steeper slope because the eccentricity drops (albeit slowly) as β is increased.

Figure 8: Close encounter times for e 3 = 0 01 . , e 3 = 0 02 . , and e 3 = 0 05, and, in purple, . e 2 = 0 05. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-dd9e1f18-43a7-4db7-a66a-9ed5fbe1b988.jpg)

Figure 9: Close encounter times for e 1 = 0 05 . , e 2 = 0 05 . , e 3 = 0 05, and . e 5 = 0 05. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-7e4f7d1a-0273-4c1e-a2be-b08a7ec2ec9d.jpg)

Figure 10: Close encounter differences between simulations with initial e 1 = 0 05 and those with . e 2 = 0 05, . e 3 = 0 05, . e 5 = 0 05, or . e 5 = f amd ( e 1 = 0 05 . , β ). This plot represents the quantitative comparison between the systems shown in Figures 7 and 9. Here, the average multiplicative factors in close encounter time compared to the e 1 = 0 05 system are, in the plotting order, 0.525, 0.436, . 0.881, and 2.393, respectively. Thus, the systems with intermediate planet eccentric (both e 2 and e 3 ) are about half as long-lived on average as those without two planetary neighbors.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-dbd171b6-00d7-41fc-8040-b4071112284e.jpg)

## 7. All planets on eccentric orbits

In the previous sections, we found that for eccentricities up to 0 05, systems with . one planet eccentric follow a log-linear behavior as a function of initial separation, and this is true no matter which planet starts out with the nonzero eccentricity. Here, we consider systems with all planets on initially eccentric orbits with e = 0 05. . In one batch, all periapses angles are initially aligned, e align = 0 05, whereas the others . have periapse angles that are randomly-drawn over specified ranges in azimuth of width 45 ◦ -360 . ◦

Figure 11 compares our results with previously-shown results for circular orbits and for e 1 = 0 05. The difference between the systems with just one planet eccentric . and all of them eccentric is stark: when all planets start out aligned and at the same eccentricity of e align = 0 05, the system lifetimes are considerably longer than those . with just one planet having e = 0 05, reaching the stopping point of five runs lasting . longer than 2 × 10 9 years at separation similar to the batches with one planet having initial e = 0 01. . In contrast, randomized initial angular parameters dramatically reduce the close encounter times, and no resonance-related structure is appearing up until at least β = 15.

It is clear that the initial difference between the eccentricity vectors of the planets (a quantity sometimes referred to as relative eccentricity of the planets) is more important than individual planetary eccentricity values. This result is in qualitative agreement with Pu and Wu (2015)'s finding that the minimum initial approach distance of orbits in a system of mildly eccentric planets is a key factor in determining the system's stability. However, as the lifetimes of the systems with aligned eccentric orbits are somewhat less than those of initially circular orbits for the largest orbital separations that we investigated, we find that for a given initial minimum separation between orbits, larger eccentricities typically reduce system lifetimes by a smallish factor. Also, the lifetimes of aligned eccentric runs are not as well-approximated by a log-linear fit as other configurations that we studied. Comparison with two other sets that reached five long-lived runs at similar values of β (to minimize effects of including differing resonances at different ranges), ( e 1 = 0 01 and . e 3 = 0 01) shows . a smaller coefficient of determination, R 2 , compared to the latter two batches, 0 803 . vs. 0 898 and 0 910, respectively. . .

Figure 11: Superposition of close encounter times for the circular systems, all planets eccentric and aligned, all planets eccentric with randomized initial periapse angles, and the e 1 = 0 05 systems. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-367dba49-6dae-4b3c-9a6a-047769ad88a7.jpg)

The substantial differences in lifetimes of systems with all planets having initial e = 0 05 for aligned and randomly-oriented initial . longitudes of periapse seen in Fig. 11 motivated us to investigate intermediate cases in which initial periapse angles are spread over a limited range in angle. Results are presented Figure 12. We see a gradual transition between the aligned and fully random results for intermediate amounts of periapse dispersion. Nonetheless, all batches approximately follow an exponential relationship between initial spacing and encounter time.

Figure 12: Superposition of close encounter times for systems with periapses aligned and those with periapses randomized within 45 ◦ , 90 ◦ , 135 ◦ , 180 ◦ and 360 . ◦

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-9cc9ddbc-85a4-48c9-961a-0f17a6d1a780.jpg)

## 8. Exponential Fitting Coefficients for Differing Batches of Run

Figure 13 plots most of the slopes and intercepts from Table 3 in the ( b , c ′ ′ ) plane. We also show lines of constant lifetimes (10 8 years and 10 10 years) for specified values of β to facilitate comparison of the stability of differing batches. Details are described in the caption. The systems with initial eccentricity of 0 01 (circular symbols, apart . from one that is magenta) lie near the line corresponding to lifetimes of 10 8 years for β = 9 5. . Similarly, the most eccentric systems in our simulations (one planet having initial e ≈ 0 05, pentagons) accumulate around the lines of highest drawn separation, . both at 10 8 years ( β = 12 5) and 10 . 10 years ( β = 12 5). . Generally, points in the lower left correspond to batches of runs that are relatively short-lived, whereas those in the upper right correspond to long-lived batches. Results for different ranges of the same set (filled, open and partly filled points of the same size, shape and color) don't differ by much, and the differences that are present don't trend in a systematic

manner with with the size of the fit region. Slopes are less steep for larger e because there is not much difference in t c for closely-spaced orbits but substantial differences in lifetime are present for large values of β . There is only a slight reduction in c ′ (the fit lifetime of the most closely-spaced systems) as initial eccentricity is increased from 0 to 0.03, but it decreases by a more significant amount when e is increased to 0.05.

The all eccentric and aligned fits are in their own region on the plot. For small separations they are almost as long-lived as initially circular runs, but at large separations the differences become greater. Overall, they aren't as well fit by a single exponential relationship as are the sets with only one planet initially eccentric.

Our findings are applicable to near-resonant pairs of exoplanets whose masses and eccentricities have been constrained using transit timing variations (TTVs) observed by the Kepler spacecraft. Hadden and Lithwick (2016) found an analytical formula for the relative eccentricity of a pair of planets orbiting near first-order mean motion resonances. At the same time, their individual eccentricities remain poorly constrained, as similarly described by Jontof-Hutter et al. (2016). Our simulations complement those results by demonstrating that systems of planets on moderately eccentric orbits are almost as stable as those on nearly circular orbits if the relative eccentricities are small, as is the case when each planet starts with e = 0 05 and . periapse angles are aligned within ∼ 45 . ◦ As a consequence, stability arguments typically cannot be used to rule out moderately large absolute eccentricities.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01117v2-figure-fecb6bc7-7e8d-4b11-b210-afa0f0722fcb.jpg)

′

Figure 13: Slopes and intercepts for the exponential fits to the lifetimes of systems in various batches of runs over the ranges in β specified below. The shape of the symbol shows the initial eccentricity of the specified planet(s), with circles representing both the circular runs and those with e = 0 01, upward-pointing triangles representing . e = 0 02, squares . e = 0 03 and pentagons . e = 0 05. We used a star shape for the all planets eccentric ( . e = 0 05) and aligned systems, since the . relevant feature for their lifetimes are not the individual planets' eccentricities, but the alignment of their periapses. Blue symbols denote batches with the innermost planet initially eccentric, grey (green) symbols denote those with the second (third) initially eccentric, red symbols signify that the outermost planet is initially eccentric, and magenta symbols have either all planets circular or all planets eccentric and aligned. Filled symbols correspond to fits over the full range from just after the widest separation system that has a collision during the first synodic period (or β = 3 47 . if no system in the batch suffered this fate) through the fifth system to survive for &gt; 2 × 10 9 years. Open and semi-open symbols represent restricted ranges, with empty symbols fit over the range 5 32 . ≤ β ≤ 10 05, left-filled . symbols relying on data from the range 4 21 . ≤ β ≤ 10 05 and . right-filled symbols for the range is 5 49 . ≤ β ≤ 12 32. .

The curves represent contours of constant β for log-linear fit lifetimes of 10 8 and 10 10 years.

## 9. Conclusions

We investigated the orbital stability of idealized systems of five identical, closelypacked planets that begin uniformly-spaced in terms of mutual Hill radii. In most of the cases considered, one planet was initially on an eccentric orbit and the other four planets began on circular orbits. As expected, system lifetimes usually increased for larger orbital separations and dropped for larger initial eccentricity. For a given initial eccentricity, we found that systems approximately obeyed a log-linear relationship between initial separation and close encounter time, in agreement with previous studies of planets starting with circular orbits, but the slopes dropped as the eccentricity increased. Mean motion resonances reduced system stability, but by smaller amounts as initial eccentricities increased.

Lifetimes of systems with the outer planet initially eccentric were on average similar to those with the inner planet starting with the same eccentricity; giving the specified initial eccentricity to an intermediate planet resulted, on average, in slightly shorter system lifetimes. The dependence of system lifetime on which planet is initially eccentric do not scale with total system angular momentum deficit, suggesting that mean motion resonances dominate over secular resonances in destabilizing the planetary systems studied.

Studies in which all planets began with an eccentricity e = 0 05 produced inter-. esting results. If the planetary periapses were selected randomly, system lifetimes were, as expected, far shorter than comparable systems with only one planet initially eccentric. However, if all of the planets' periapses were initial aligned or nearly so, systems were generally much longer-lived than systems in which only one planet had initial e = 0 05. . For initially fully aligned systems with all planets beginning with e = 0 05, planetary spacing required for survival . &gt; 10 9 years is similar to those for which four of the planets begin on circular orbits and the other planet starts with e = 0 01. .

## 10. Acknowledgements

We thank Nicolas Faber, Eric Feigelson, Sam Hadden, Daniel Jontof-Hutter, Yoram Lithwick, Fred Rasio, David Rice, Jason Steffen, Daniel Tamayo, and Christa Van Laerhoven for useful comments and discussions.

PGwas supported by NASA-Ames / Fonds National de la Recherche Luxembourg Grant SP0044174/60048152. JJL's work on this project was supported by NASA's Astrophysics Data Analysis Program 16-ADAP16-0034.

## 11. References

## References

Chambers, J. E., Wetherill, G. W., and Boss, A. P. (1996). The Stability of MultiPlanet Systems. Icarus , 119:261-268. doi: https://doi.org/10.1006/icar. 1996.0019 .

Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., Rowe, J. F., Steffen, J. H., Agol, E., Barclay, T., Batalha, N., Borucki, W., Ciardi, D. R., Ford, E. B., Gautier, T. N., Geary, J. C., Holman, M. J., Jenkins, J. M., Li, J., Morehead, R. C., Morris, R. L., Shporer, A., Smith, J. C., Still, M., and Van Cleve, J. (2014). Architecture of Kepler's Multi-transiting Systems. II. New Investigations with Twice as Many Candidates. ApJ , 790:146. doi: https://doi.org/10.1093/mnras/sty2418 .

Gladman, B. (1993). Dynamics of systems of two close planets. Icarus , 106:247. doi: https://doi.org/10.1006/icar.1993.1169 .

Hadden, S. and Lithwick, Y. (2016). Numerical and analytical modeling of transit timing variations. The Astrophysical Journal , 828(1):44.

Hill, G. W. (1878). Researches in the lunar theory. American journal of Mathematics , 1(1):5-26. doi: https://doi.org/10.2307/2369430 .

Jontof-Hutter, D., Ford, E. B., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., Van Laerhoven, C., Agol, E., Deck, K. M., Holczer, T., and Mazeh, T. (2016). Secure mass measurements from transit timing: 10 kepler exoplanets between 3 and 8 m ⊕ with diverse densities and incident fluxes. The Astrophysical Journal , 820(1):39.

Laskar, J. (1997). Large scale chaos and the spacing of the inner planets. Astronomy and Astrophysics , 317:L75-L78.

Laskar, J. and Petit, A. (2017). Amd-stability and the classification of planetary systems. Astronomy &amp; Astrophysics , 605:A72.

Lissauer, J. J., Fabrycky, D. C., Ford, E. B., Borucki, W. J., Fressin, F., Marcy, G. W., Orosz, J. A., Rowe, J. F., Torres, G., Welsh, W. F., Batalha, N. M., Bryson, S. T., Buchhave, L. A., Caldwell, D. A., Carter, J. A., Charbonneau, D., Christiansen, J. L., Cochran, W. D., Desert, J.-M., Dunham, E. W., Fanelli, M. N., Fortney, J. J., Gautier III, T. N., Geary, J. C., Gilliland, R. L., Haas, M. R., Hall, J. R., Holman, M. J., Koch, D. G., Latham, D. W., Lopez, E., McCauliff, S., Miller, N., Morehead, R. C., Quintana, E. V., Ragozzine, D., Sasselov, D., Short, D. R., and Steffen, J. H. (2011a). A closely packed system of low-mass,

low-density planets transiting kepler-11. Nature , 470(7332):53. doi: https:// doi.org/10.1038/nature09760 .

Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., Steffen, J. H., Ford, E. B., Jenkins, J. M., Shporer, A., Holman, M. J., Rowe, J. F., Quintana, E. V., Batalha, N. M., Borucki, W. J., Bryson, S. T., Caldwell, D. A., Carter, J. A., Ciardi, D., Dunham, E. W., Fortney, J. J., Thomas N. Gautier, I., Howell, S. B., Koch, D. G., Latham, D. W., Marcy, G. W., Morehead, R. C., and Sasselov, D. (2011b). Architecture and dynamics of kepler's candidate multiple transiting planet systems. The Astrophysical Journal Supplement Series , 197(1):8. doi: https://doi.org/10.1088/0067-0049/197/1/8 .

Marzari, F. and Weidenschilling, S. (2002). Eccentric extrasolar planets: the jumping jupiter model. Icarus , 156(2):570-579. doi: https://doi.org/10.1006/icar. 2001.6786 .

Morrison, S. J. and Kratter, K. M. (2016). Orbital stability of multi-planet systems: behavior at high masses. The Astrophysical Journal , 823(2):118. doi: https: //doi.org/10.3847/0004-637X/823/2/118 .

Obertas, A., Van Laerhoven, C., and Tamayo, D. (2017). The stability of tightlypacked, evenly-spaced systems of earth-mass planets orbiting a sun-like star. Icarus , 293:52 - 58. doi: https://doi.org/10.1016/j.icarus.2017.04.010 .

Petit, A. C., Laskar, J., and Bou´, G. (2018). e Hill stability in the amd framework. arXiv preprint arXiv:1806.08869 .

Pu, B. and Wu, Y. (2015). Spacing of kepler planets: Sculpting by dynamical instability. The Astrophysical Journal , 807(1):44. doi: https://doi.org/10. 1088/0004-637X/807/1/44 .

Quarles, B. and Lissauer, J. J. (2018). Long-term stability of tightly packed multiplanet systems in prograde, coplanar, circumstellar orbits within the α centauri ab system. The Astronomical Journal , 155(3):130. doi: https://doi.org/10.3847/ 1538-3881/aaa966 .

Rein, H. and Liu, S.-F. (2012). Rebound: an open-source multi-purpose n-body code for collisional dynamics. Astronomy &amp; Astrophysics , 537:A128. doi: https: //doi.org/10.1051/0004-6361/201118085 .

Rein, H. and Tamayo, D. (2015). Whfast: a fast and unbiased implementation of a symplectic wisdom-holman integrator for long-term gravitational simulations.

Monthly Notices of the Royal Astronomical Society , 452(1):376-388. doi: https: //doi.org/10.1093/mnras/stv1257 .

Rice, D. R., Rasio, F. A., and Steffen, J. H. (2018). Survival of non-coplanar, closely packed planetary systems after a close encounter. Monthly Notices of the Royal Astronomical Society , 481(2):2205-2212. doi: https://doi.org/10.1093/mnras/ sty2418 .

Rowe, J. F., Bryson, S. T., Marcy, G. W., Lissauer, J. J., Jontof-Hutter, D., Mullally, F., Gilliland, R. L., Issacson, H., Ford, E., Howell, S. B., Borucki, W. J., Haas, M., Huber, D., Steffen, J. H., Thompson, S. E., Quintana, E., Barclay, T., Still, M., Fortney, J., Gautier, III, T. N., Hunter, R., Caldwell, D. A., Ciardi, D. R., Devore, E., Cochran, W., Jenkins, J., Agol, E., Carter, J. A., and Geary, J. (2014). Validation of Kepler's Multiple Planet Candidates. III. Light Curve Analysis and Announcement of Hundreds of New Multi-planet Systems. ApJ , 784:45. doi: https://doi.org/10.1088/0004-637X/784/1/45 .

Smith, A. W. and Lissauer, J. J. (2009). Orbital stability of systems of closely-spaced planets. Icarus , 201:381-394. doi: https://doi.org/10.1016/j.icarus.2008. 12.027 .

Zhou, J.-L., Lin, D. N., and Sun, Y.-S. (2007). Post-oligarchic evolution of protoplanetary embryos and the stability of planetary systems. The Astrophysical Journal , 666(1):423.
Pierre Gratia Jack J. Lissauer
Eccentricities and the Stability of Closely-Spaced Five-Planet Systems
Observations of exoplanets have revealed that systems with planets on
closely-spaced orbits are common, which motivates the question "How closely can
planets orbit to one another and still be dynamically-stable for very long
times?". To address this question, we investigate the stability of idealized
planetary systems consisting of five planets, each equal in mass to the Earth,
orbiting a one solar mass star. All planets orbit in the same plane and in the
same direction, and the planets are uniformly spaced in units of mutual Hill
Sphere radii. Most of the systems that we integrate begin with one or more
planets on eccentric orbits, with eccentricities $e$ as large as $e= 0.05$
being considered. For a given initial orbital separation, larger initial
eccentricity of a single planet generally leads to shorter system lifetime,
with little systematic dependence of which planet is initially on an eccentric
orbit. The approximate trend of instability times increasing exponentially with
initial orbital separation of the planets found previously for planets with
initially circular orbits is also present for systems with initially eccentric
orbits. Mean motion resonances also tend to destabilize these systems, although
the reductions in system lifetimes are not as large as for initially circular
orbits. Systems with all planets having initial $e= 0.05$ and aligned periapse
angles typically survive far longer than systems with the same spacing in
initial semi-major axis and one planet with $e= 0.05$, but they have slightly
shorter lifetimes than those with planets initially on circular orbits.