Dynamical instabilities in systems of multiple short-period planets are likely driven by secular chaos: a case study of Kepler-102

The Kepler space telescope carried out a large systematic exoplanet survey during the Kepler (Borucki et al., 2010) and K2 missions (Howell et al., 2014), and has provided a wealth of data on planets and planetary systems in the Galaxy. A large subset, about 40%, of the Kepler planet candidates are in systems with two or more planet candidates (e.g., Coughlin et al., 2016). Studies of the Kepler data on multiple-planet systems have concluded that planetary systems are typically coplanar to within a few degrees (e.g., Lissauer et al., 2011; Tremaine & Dong, 2012; Fang & Margot, 2012; Johansen et al., 2012; Figueira et al., 2012; Fabrycky et al., 2014) and that they generally have low eccentricities as well (e.g., Xie et al., 2016; Van Eylen & Albrecht, 2015). In contrast with the Solar system, most of the Kepler systems are tightly packed, within less than 0.5 au of their stellar host (e.g., Lissauer et al., 2011; Fressin et al., 2013). Theoretical studies of dynamical stability of systems with similar masses and orbital spacings to those of the observed Kepler systems conclude that such systems are close to the threshold of instability (e.g., Fang & Margot, 2013). Pu & Wu (2015) and Volk & Gladman (2015) have suggested this could be because inner planet systems form with even tighter spacings and higher multiplicity, then undergo a sequence of dynamical instabilities that pare down the numbers of planets and widen their orbital spacings, a process involving planet-planet collisions and mergers.

The orbital spacings in planetary systems are usefully described in units of the mutual Hill radius ($R_{mH}$) of adjacent planets. For two planets of masses $m_{1}$ and $m_{2}$ orbiting a star of mass $M_{*}$ in orbits of semimajor axes $a_{1}$ and $a_{2}$, the parameter $K$ describes the number of mutual Hill radii separating the planetary orbits:

The dynamical spacing requirement for stable co-planar, low-eccentricity two-planet systems is relatively simple and closely follows that for the coplanar restricted three body problem (see, e.g., Gladman, 1993, and references therein),

In contrast with two planet systems, the stability in higher-multiplicity systems is poorly understood and is generally expected to be of a statistical nature, although we can expect that the minimum dynamical separation of Eq. 2 must be satisfied (Malhotra, 2015). Empirical scaling laws for the relationship between stability times and dynamical separations have been derived with numerical integrations of evenly-spaced, equal-mass planets (e.g., Chambers et al., 1996; Smith & Lissauer, 2009; Funk et al., 2010; Obertas et al., 2017) as well as for inhomogeneous systems (that is, those with a dispersion in planet masses, orbital spacings and orbital eccentricities and inclinations, Pu & Wu (2015)). However, these empirical scaling relationships have limited applicability to the observed multi-planet distribution because they mostly apply to instabilities at dynamical separations of $K\lesssim$ 10–15, whereas the observed multi-planet systems have a much broader distribution of $K$ (see Figure 1). Additionally, although empirical rules for dynamical (in)stability have been found in numerical simulations, there is relatively little detailed understanding of the underlying mechanisms that generate the instabilities in simulated planetary systems.

In the present work, we expand upon these previous studies by using the observed Kepler and K2 planetary architectures as a starting point for studying the underlying source of dynamical stability and instability in the diverse set of observed compact systems of small planets. Figure 1 shows the estimated $K$ distribution of all the currently catalogued adjacent pairs of planets observed during the Kepler and K2 missions in systems with 4 or more planets (population statistics taken from the NASA Exoplanet Archive¹¹1exoplanetarchive.ipac.caltech.edu). The observational data can be used to measure the fractional orbital separations from the orbital periods, but do not directly provide planet masses. Therefore, to compute the $K$ values, we use the Wolfgang et al. (2016) statistical mass-radius relationship to assign a range of possible masses to each observed planet based on the planet radius reported in the NASA Exoplanet Archive’s Composite Planet Data table. For each observed planetary system, we generate 100 versions of that system from the mass-radius relationship and calculate $K$ for each adjacent pair in the system. The resulting histogram of $K$ values (in which each observed planet pair has a weight of one) is shown in Figure 1.

Considering that most of the Kepler multi-planet host stars are several gigayears old (e.g., Silva Aguirre et al., 2015), it is notable that so many planet pairs appear to have quite small values of $K$, $K\lesssim 15$, where the above-mentioned scaling laws would suggest dynamical lifetimes shorter than gigayears. Perhaps some of the planets at small spacings are spurious detections or have masses smaller than we have assumed. However, we find that even large dynamical spacings do not guarantee stability. In our numerical simulations of these systems (described in Section 2), $\sim 20\%$ of simulated systems became unstable on timescales of a few billion orbits ($\lesssim 100$ Myr). Approximately half of the instabilities occur between planet pairs that begin with $K>15$ (see the green histogram in Figure 1), separations predicted to be stable based on the simple scaling laws mentioned above.

This motivates us to investigate the dynamics of these planetary systems in order to identify the source of dynamical instabilities as well as markers for dynamically stable architectures. The rest of this paper is organized as follows. In Section 2, we describe our simulation results for the distribution of dynamical lifetimes of multi-planet systems spanning the range of orbital architectures similar to those of the observed sample of Kepler multis. In Section 3, we describe our case study of the Kepler-102 system: we show with analytical and numerical estimates that the overlap of mean motion resonances is unlikely to be the direct cause of dynamical instabilities in this system, we investigate the secular dynamics in this system, and we use an analysis of the system’s angular momentum deficit (AMD) to show that secular chaos is the likely cause of dynamical instabilities in simulations of this system. In Section 4 we describe the usefulness of spectral analysis of relatively short simulations to diagnose long-term stable and unstable planetary architectures; based on this analysis, we tentatively conclude that secular chaos likely drives instability in simulations of many of the observed Kepler and K2 multiplanet systems. We summarize in Section 5.

To investigate the potential drivers of instability in realistic multiplanet systems, we performed a large suite of simulations based on the observed properties of the Kepler and K2 systems with four or more confirmed planets and an estimated stellar mass. A complete analysis of these simulations is deferred to a future paper, but here we use some of their basic results to motivate a case study (Section 3.1) of the detailed dynamics that drive instabilities in multi-planet systems.

For each observed planetary system in our sample, we use the properties of the system listed in the NASA Exoplanet Archive’s Composite Planet Data table combined with a statistical mass-radius relationship as the basis for our simulation initial conditions. From this database, we take the stellar mass, planetary radius, and orbital periods (along with all associated uncertainties) for each system. The planetary radii and uncertainties are fed into the Wolfgang et al. (2016) statistical mass-radius relationship²²2Their code is available at github.com/dawolfgang/MRrelation to generate 100 masses for each planet; this includes 50 masses calculated using Wolfgang et al. (2016)’s transit timing variation priors and 50 calculated using the radial velocity priors. For each generated set of planetary masses, a stellar mass is selected uniformly from within the uncertainties (to ensure we fully sample the most likely stellar mass range with our relatively small set of simulations), and the observed planetary orbital periods are then used to calculate a semimajor axis for each planet. The eccentricities and inclinations of each planet are then selected from rayleigh distributions of widths $\sigma_{e}=0.02$ and $\sigma_{i}=1.4^{\circ}$, respectively; these widths are consistent with the estimated intrinsic eccentricity and inclination widths for Kepler multiplanet systems (e.g., Hadden & Lithwick, 2014; Xie et al., 2016; Van Eylen & Albrecht, 2015). The orbital angles (argument of perihelion, longitude of ascending node, and mean anomaly) for each planet are randomized in the range 0 to $2\pi$. These initial conditions do not represent the exact Kepler systems, but represent plausible variations of each observed planetary system.

The resulting $\sim 6000$ sets of planetary system initial conditions were then integrated for $5\times 10^{9}$ orbits of the innermost planet (corresponding to $\sim 10-100$ Myr timescales, depending on the system). The integrations were performed using the mercurius integrator (based on Chambers 1999) within rebound (Rein & Liu, 2012). A post-Newtonian general relativity correction was added to the code (using the same prescription as in Lissauer et al., 2011). The initial timestep was set to 1/40th of the innermost planet’s orbital period, and the integrator switches to an adaptive timestep routine to handle close encounters between planets. The integrations were stopped if two planets collided (the collision radius was set to the physical radius of the planets).

Overall, $\sim 20\%$ of these simulations resulted in a collision between planets in the system. As noted by Volk & Gladman (2015), the distribution of instability timescales is linear in log time (Figure 2). We find that systems with dynamical spacings $K\lesssim 12-15$ have the strongest correlation between $K$ and the probability of dynamical instability; this is consistent with the stability scaling laws discussed in Section 1. However, at larger values of $K$ (where the vast majority of the observed multiplanet systems appear to lie) the relationship between stability and dynamical spacing is not clear. We find that, for $K\gtrsim 15$, neither the fraction of surviving systems nor the instability timescale appears to be strongly correlated with dynamical spacing. Analysis of our suite of simulations shows that the most consistent correlation between any simple dynamical parameter and the probability that the system experiences instability is with the period ratios of adjacent planet pairs. Systems that have planet pairs with period ratios below 2 are the ones that most often display instabilities, as seen in Figure 3. In the next section, we investigate a single planetary system in detail to explore what these trends from the larger simulation set mean for the source of instabilities.

To identify the triggers of long-term instability in multi-planet systems, we investigate in detail the five-planet Kepler-102 system. We chose Kepler-102 for three reasons: it has five confirmed planets (and thus could display rich dynamical behavior over the sampled parameters); our simulations of this system had a high incidence of instabilities; and the observed period ratios of the confirmed planets suggest mean motion resonances could be important. Given this last point and that the probability of our simulated planetary systems experiencing an instability is not strongly correlated with dynamical spacings and only appears to be related to the period ratios of adjacent planet pairs, we first investigate whether proximity to or overlap between mean motion resonances is a viable explanation for instabilities in this system (Section 3.2); however, for Kepler-102 (and later for a wider range of systems), we find that this is not a likely explanation for reasonable planet mass ranges. Instead, the secular structure of the simulated planetary systems appears to be much more predictive of instability. In Section 3.3, we describe evidence that secular chaos causes the transfer of angular momentum deficit amongst planets within the planetary system such that even initially low-eccentricity planets can be driven to orbit-crossing and collisions.

The analysis of the Kepler-102 system in the previous section has provided useful insights into possible drivers of dynamical instability, but does not provide a way to predict long term stability. Here we use the results of our large suite of long integrations supplemented with a suite of short integrations of the same systems, to identify a tool to predict stability/instability from the short simulations. This tool is based on the concept of spectral entropy of conservative dynamical systems introduced in Noid et al. (1977) and Powell & Percival (1979). Briefly, the fast Fourier transform of the time series from a numerical integration of a dynamical system is dominated by a small number of frequency components if the system is regular (stable), but it has many weaker frequency components if the system is chaotic (unstable); thus, the number of significant frequencies (“spectral number") in the power spectrum is diagnostic of (in)stability. This method has been recently used to investigate orbital stability and chaos of stellar orbits in the Galaxy (Lépine et al., 2017, and references therein) and to investigate regions of regular and chaotic orbital motion in a few individual exoplanet systems (Tadeu dos Santos et al., 2012; Alves et al., 2016). In the context of the solar system, Michtchenko et al. (2002) employed this method to diagnose chaotic diffusion of the asteroid 1459 Magnya and its collisional family, and Kotoulas & Voyatzis (2004) used this method to map chaotic regions in the vicinity of Neptune’s mean motion resonances in the Kuiper belt. Below, we employ this method to demonstrate that the power spectrum of the time series of the planets’ AMD, eccentricity, and inclination from short integrations can reveal whether a system is long term stable.

For each of the multi-planet Kepler and K2 systems that we simulated in our long-term integration set (Section 2), we re-integrated each set of initial conditions over $5\times 10^{6}$ initial orbital periods of the innermost planet, recording the orbital elements and AMD of each planet at 3000 evenly spaced points. This time span is sufficiently long to capture multiple cycles of the longest-period linear secular modes in the vast majority of these systems. We then performed a fast Fourier transform, FFT (using the numpy fft package in python), of the time series of the semimajor axis, inclination, eccentricity, and AMD of each planet in each simulation. For an input time series with 3000 points evenly spaced in time, this yields an estimate of the power associated with 1500 evenly linearly spaced frequencies where the highest frequency is twice the sampling frequency and the total length of the time series determines the lowest frequency. Following Michtchenko et al. (2002), we characterize the results of the FFT by considering how many frequencies in the FFT have power $\geq 5\%$ of the peak frequency in the power spectrum. While Michtchenko et al. (2002) considered the absolute number of frequencies above this threshold (spectral number) as their parameter of interest, here we make a slight modification to instead define a “spectral fraction", i.e., the fraction of the frequencies in the FFT that exceed the $5\%$ threshold; this ensures our parameterization does not depend on the length or cadence of the input time series (because the number of entries determines the number of frequencies in the FFT).

Figure 8 shows the FFT of the AMD of planet b in two versions of the Kepler-102 system. In the longer integrations, the case shown in the left panel survived the full $5\times 10^{9}$ orbits, while the case on the right did not. In the short integration, the FFT of the stable case has only a few frequencies that rise above the 5% of the peak power threshold and its spectral fraction is $2\times 10^{-3}$; these few frequencies occur in a well defined single peak in the FFT close to one of the linear secular modes of this system. However, in the short integration of the unstable case there are a significantly larger number of spectral frequencies above the 5% threshold, and its spectral fraction is $4.5\times 10^{-2}$. We see that these two cases have qualitatively different power spectra in the short integrations, and the spectral fraction is small for the stable case and larger for the unstable case.

We note that small spectral fractions for the AMD evolution of a planet is what we would expect in cases where linear secular theory accurately describes the evolution of a planet’s orbit; i.e., when the eccentricity and inclination evolution are decoupled and can be described as a sum of a limited number of secular frequencies. When this is the case, the AMD evolution (which reflects a combination of the eccentricity and inclination evolution, see equation 4) will have approximately the same number of dominant frequencies that are similar in magnitude to the secular frequencies. For reference, we show the predicted linear secular frequencies for our simulated Kepler-102 systems in Figure 8. In the stable system, we can see that there are relatively few local peaks in planet b’s AMD power spectrum and that they occur at frequencies similar to the range of the secular frequencies; the same cannot be said of the unstable system’s AMD power spectrum for planet b.

We calculated the spectral fractions based on short integrations for all of the planets in all of the simulated Kepler and K2 systems to assess how well they predict stability or instability in the longer simulations. We exclude from our analysis cases that went unstable on timescales shorter than $3\times 10^{6}$ orbits as the power spectra from such short simulations will not include some of the longer period secular modes (such short-lived systems also reveal themselves with very little cpu time and thus do not need predictors of instability). How system stability in our long simulations depends on the spectral fractions from the short simulations is summarized in Figure 9. To construct each panel of Figure 9, we binned our simulations according to the largest spectral fraction from the AMD time series of any individual planet in the system (x axes) and the largest spectral fraction for the semimajor axis time series (y-axes; left panels), the eccentricity time series (y-axes; middle panels), and the inclination time series (y-axes; right panels). The color assigned to each spectral fraction bin indicates the fraction of systems in the bin that survived for $5\times 10^{9}$ orbits of the innermost planet in our long simulations. (Using the average instead of the largest spectral fraction yields very similar plots.)

For reference, we also ran a numerical integration of the solar system’s eight major planets and computed the spectral fractions for the giant planets (calculated over a timespan equal to $5\times 10^{6}$ Jupiter orbits) and terrestrial planets (calculated over a timespan equal to $5\times 10^{6}$ Mercury orbits); the maximum spectral fractions for the giant planets and for the terrestrial planets are indicated by ‘g’ and ‘t’, respectively, in Figure 9. This is a useful comparison because long-term simulations of the giant planets have shown their orbits to be stable (even in cases where their maximum Lyapunov exponents indicate chaotic behavior, see Hayes 2008), whereas the current orbits of the terrestrial planets do allow for potential future collisional trajectories (Laskar & Gastineau, 2009) that might be attributable to secular chaos (Lithwick & Wu, 2014).

The top panels of Figure 9 display the results for the entire set of planetary systems described in Section 2, and it is immediately apparent that the stable planetary systems are not randomly distributed in spectral fraction. The stable systems are concentrated at small AMD spectral fractions; they also tend toward low spectral fractions in eccentricity and inclination (which are both also represented in the AMD analysis). There is not a strong trend with the spectral fraction of the semimajor axis; this is unsurprising as non-resonant planets should not have strong features in their semimajor axis power spectra.

Nevertheless, there are a few stable systems in the top panels of Figure 9 that do not have small spectral fractions in AMD, $e$, and/or $i$. Conjecturing that these outliers are systems that are affected by mean motion resonances, we took our set of simulation initial conditions and determined which cases had pairs of planets closer than 1.5 resonance widths to mutual first-order resonances (using equation 3 to calculate widths for low-eccentricity orbits consistent with our simulation initial conditions). We excluded these cases (which represent instances only $\sim 5$ of the $\sim 60$ observed planetary system architectures on which the simulations are based). The stability of the remaining systems as a function of spectral fraction are shown in the bottom panels of Figure 9. This removes the majority of the clear outliers; it is then clear that in systems not strongly influenced by MMRs, a spectral fraction above $\sim 0.01-0.02$ in short simulations is correlated with a much lower chance of stability in long simulations. This supports the hypothesis that secular chaos is a significant driver of long term instability. The location of the solar system’s terrestrial planets near the secular stability boundary in the spectral fraction stability maps also supports this conclusion.

As a final step, we examine whether the stability timescales of the non-resonant planetary systems (from the bottom panels of Figure 9) are correlated with the AMD spectral fractions. The median survival time for a planetary system (in the long simulations) as a function of the AMD spectral fraction (from the short simulations) is shown in Figure 10. We find a rather sharp transition at a spectral fraction of 0.01: below this value all of the cases are stable for the full $5\times 10^{9}$ orbital periods of the long simulations; above this value, survival times generally decrease.

The strong correlations between long term dynamical stability and AMD spectral fraction support the hypothesis that secular chaos is an important driver of the evolution of planetary systems. The “spectral fraction” approach is a very promising tool for quickly diagnosing the likelihood of dynamical instabilities in a planetary system. The short simulations on which the spectral fraction calculations are based typically take only $\sim 5-10$ minutes per system to complete on a desktop computer, compared to the hundreds of cpu hours required to integrate an individual system for billions of orbits. A spectral fraction parameter might be particularly valuable in attempts to predict planetary system stability using machine learning (e.g., Tamayo et al., 2016).

We have performed a large suite of numerical simulations of planetary systems based on the observational sample of Kepler and K2 transiting systems with four or more confirmed planets. We find that $\sim 20\%$ of the cases show dynamical instabilities (leading to planet-planet collisions) within $5\times 10^{9}$ orbital periods of the innermost planet. This result indicates that dynamical instabilities are not uncommon for plausible variations on the observed planetary architectures. We find that these instabilities occur at a wide range of the planets’ dynamical spacings, indicating that the source of instabilities in many of these simulations is not simply a result of planets being closely packed. We find a wide range of instability timescales in our simulations. Correlation analysis of the stable versus unstable systems reveals only that the unstable systems typically have smaller period ratios of adjacent planet pairs (Section  2).

The bulk simulations motivated us to take a close look at what drives the instabilities in these systems. In our case-study of the Kepler-102 system (Section 3), we find that mean motion resonances are unlikely to be the initial trigger of instability. Despite period ratios that appear close to low-integer ratios, for plausible values of planet masses and orbital eccentricities, the widths of mean motion resonances are too narrow for resonance overlap to occur and to directly influence dynamical stability (Section 3.2). However, we do note that only modest eccentricity increases (up to $e\sim 0.1$) are required for the 4:3 resonance between the innermost pair of the Kepler-102 system to affect their dynamics. In Section 3.3, we report evidence that the inward transfer of AMD via secular chaos is the most likely driver of eccentricity growth in the orbits of Kepler-102’s inner planets in our simulations. An estimate of the total system AMD required to raise the eccentricity of Kepler-102 b enough to induce instability provides an approximate lower limit of $\sim 0.1M_{\earth}$ on the mass of planet b (Fig. 6). This conclusion is supported by the numerical evidence of inward transfer of AMD in long term unstable versions of the Kepler-102 system (Fig. 7). Even for very low initial eccentricities ($e\sim 0.02$), both the AMD-stability estimate and the long-term simulations indicate that cases where Kepler-102 b’s mass exceeds $\sim 0.1M_{\earth}$ are more favorable for long-term stability. This lower limit on planet b’s mass corresponds to Earth-like bulk densities, and excludes the much lower densities allowed by statistical mass-radius relationships.

In Section 4, we performed a frequency analysis of the AMD and the orbital elements for the planets in our simulated systems in short integrations (spanning a few secular cycles) and classified them according to a "spectral fraction" that quantifies whether a planet’s secular evolution is dominated by a few well-defined frequencies (small spectral fraction) or has a noisy power spectrum (large spectral fraction). We find that for planetary systems lacking first order mean motion resonances (which is the majority of the mutiplanet systems considered here), small spectral fractions (below $\sim 0.01$) of the AMD calculated from integrations of a few million orbital periods are strongly associated with long term stability (found in much more cpu-intensive simulations spanning billions of orbits). We also find that the median stability timescale generally decreases with increasing AMD spectral number. This supports the hypothesis that secular chaos (which allows the transfer of AMD between planets via near-degenerate or resonant secular frequencies) is the dominant driver of instabilities in many multiplanet systems. We conclude that the spectral fraction approach also provides a promising tool to predict the longer-term stability of planetary systems by means of computationally inexpensive short term numerical simulations.

Acknowledgements: This work was supported by NASA (grant 80NSSC18K0397). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. An allocation of computer time from the UA Research Computing High Performance Computing (HPC) at the University of Arizona is gratefully acknowledged. This paper includes data collected by the Kepler and K2 missions. Funding for the Kepler and K2 missions is provided by the NASA Science Mission directorate.