The Dynamics of Tightly-packed Planetary Systems in the Presence of an Outer Planet: Case Studies Using Kepler-11 and Kepler-90

The Kepler mission has demonstrated that planetary systems with multiple planets on short orbital periods are common. Roughly 90% of the Kepler candidates are found on orbits less than about 90 days (Burke et al. 2014), with 23% of Kepler stellar hosts containing more than one candidate. In terms of absolute distances and periods, these planetary configurations are much more compact than the solar system's architecture, in which Mercury orbits the Sun in 88 days. Approximately 46% of all candidates reside in known multi-planet systems (Burke et al. 2014). Based on survey results, the frequency of planets at short orbital periods around solar-type stars is between about 30% and 50% (Mayor et al. 2011; Howard et al. 2012), suggesting that these systems may be present around at least 5% of stars. While this may not represent the majority of planetary systems after a billion or more years of evolution, these Systems with Tightly-packed Inner Planets (STIPs) may represent a very common mode of planet formation, with instability and high collisional speeds potentially leading to system decay over long timescales (Volk & Gladman 2015). During planet building, migration of outer planets could frustrate the formation of planets at short orbital periods (Batygin & Laughlin 2015), even if their formation would otherwise be common. Moreover, the presence of outer planets can affect the evolution of inner planetary systems, particularly if the outer system becomes unstable (Mustill et al. 2017). Such a situation has been extensively studied in the context of the solar system's evolution (Brasser et al. 2009, 2013; Agnor & Lin 2012; Kaib & Chambers 2016).

STIPs offer a range of planetary configurations that can be used to test planet formation hypotheses. For example, at face value their low-eccentricity, closely-spaced configurations are consistent with disk migration. However, there is an under abundance of planets in or near commensurabilities, which is a fundamental prediction of convergent disk migration (e.g., Hands et al. 2014) without additional processes. Dynamical migration, such as planet–planet scattering followed by tidal evolution (e.g., Rasio & Ford 1996; Chatterjee et al. 2008; Ford & Rasio 2008), cannot obviously explain the known orbital architectures because planets in STIPs have low eccentricities and mutual inclinations (Fabrycky et al. 2014). The situation becomes further complicated when considering star–planet spin–orbit alignments. Most STIPs with spin–orbit measurements show that the stellar spin axes are aligned with the normals to the planetary orbital planes (Fabrycky & Winn 2009; Southworth 2011). Yet for several systems (e.g., HAT-P-07, KELT-17, Kepler-63, WASP-08, WASP-33), the stellar spin obliquity is large. Such misalignments could be caused by processes intrinsic to the star (Rogers et al. 2012, 2013), but it is not yet clear whether this can explain the observed population. Non-trivial disk-star or disk-environment interactions could instead force the planet-forming disk to be misaligned during planet building (e.g., Bate et al. 2010; Batygin 2012; Crida & Batygin 2014; Lai 2014). The misalignment may also reflect pure dynamical interactions, in which a distant and highly misaligned companion causes an entire STIP to oscillate in inclination while keeping the mutual inclinations of the planets in the STIP small (Batygin et al. 2011; Kaib et al. 2011; Boué & Fabrycky 2014a, 2014b). In this sense, the STIP oscillates as a rigid system. The frequency of oscillation is a function of the total mass of the inner planets, the mass of the outermost body, and the separation of the perturber and the STIP. The Lidov–Kozai mechanism, which would normally cause coupled changes in a planet's orbital eccentricity and inclination, is quenched by planet–planet interactions (Innanen et al. 1997; Takeda et al. 2008).

Finally, secular dynamics may also play a fundamental role in shaping the observed population of planetary systems, which is the subject of this paper. Secular eccentricity resonances can give rise to seemingly spontaneous system instability (Hansen & Murray 2015; Volk & Gladman 2015), causing planet mergers and potentially erosion. Secular inclination resonances can place some planets on large inclinations, without the presence of a highly inclined perturber. This could cause the star to have a large obliquity relative to some planets in the system and could ensure that only a fraction of the system would be observed by transits directly. In this paper, we seek to address the following questions: Can planets in a STIP exhibit large mutual inclination variations while remaining otherwise stable? Can we find plausible examples of this mechanism in action? We also explore the overall stability of STIPs in the presence of a single outer planet for a range of orbital distances.

This paper is organized as follows: In Section 2, we summarize general considerations for stability of planetary systems based on their multiplicity and planet separations. In Section 3, we describe the experimental overview and methodology used throughout this work. Section 4 provides the results from our simulations and secular analyses. We discuss our results in Section 5 and, finally, summarize our findings in Section 6.

A pair of planets in a two-planet system will always be stable if they are separated by ${\rm{\Delta }}\gt 2.4\,{({\mu }_{0}+{\mu }_{1})}^{1/3}$ (Gladman 1993), where μ₀ and μ₁ are the masses of the inner and outer planet relative to the star, respectively. The second planet's semimajor axis is thus ${a}_{1}={a}_{0}(1+{\rm{\Delta }})$, where a₀ is the inner planet's semimajor axis. Planet separations can also be described in units of mutual Hill radii, where the mutual Hill radius ${R}_{{mH}}\,=0.5\,({a}_{0}+{a}_{1}){[({\mu }_{0}+{\mu }_{1})/3]}^{1/3}$. This gives $\eta \,\equiv ({a}_{1}-{a}_{0})/{R}_{{mH}}$, for which two-planet stability requires η ≳ 3.46/(1 + Δ/2). For example, a two Jupiter-mass planet system around a solar-mass star must have planet separations of $3\ {R}_{{mH}}$ or greater.

Unfortunately, there is no planet separation limit to ensure stability at all times when the number of planets ${N}_{P}\gt 2$. In these cases, we can only give a typical timescale for the planets to become orbit crossing (e.g., Obertas et al. 2017). This timescale depends on the initial η and the number of adjacent planets with that η (e.g., Chatterjee et al. 2008). For Earth-mass planets with N_P ≳ 3, a mutual Hill radius separation of 10 can allow long-term stability (Smith & Lissauer 2009). Such spacings and planet multiplicities are relevant to STIPs: based on the confirmed systems in the NASA Exoplanet Archive,¹ there are 74 systems with a multiplicity N_P > 3, 26 with N_P > 4, 6 with N_P > 5, and 2 with N_P > 6. Few STIPs have planets with η < 10 , suggesting that the known planets should be stable over the lifetime of most stars (Lissauer et al. 2014a; Obertas et al. 2017). However, secular interactions can lead to the disruption of a system, even if the planetary spacing alone suggests long-term metastability. A classic example is the evolution of Mercury in the solar system, which has a small but non-negligible probability of being driven to orbit crossing with Venus (Laskar 1994).

Because many STIPs have high multiplicity (which we take to be N_P > 3 in this paper), we must ask whether a large fraction of the systems with only two or three planets are decay products themselves.² There are also observability considerations. In particular, the presence of outer perturbers can affect the observed planet multiplicity of transiting systems, which can have bearing on how we interpret planet-star misalignment (e.g., Winn et al. 2005; Kaib et al. 2011; Boué & Fabrycky 2014a, 2014b), at least in part.

We investigate the observability and stability of STIPs by using direct numerical integration along with secular theory. First, we examine the behavior of Kepler-11 (K11) and Kepler-90 (K90, also known as KIC 11442793 and KOI-351) with and without an additional, perturbing Jupiter-like planet. When included, this perturber is placed exterior to the last known planet in each system, initially at 5.2 au. Because the location of the Jupiter analogue affects the secular frequencies, we "migrate" the planet inwards to sample a broader range of frequency space.

K11 (Lissauer et al. 2011a, 2013) is composed of six planets that have orbital semimajor axes less than 0.5 au. The mass and orbital parameters of K11 planets are well defined, with the exception of the outermost planet, K11g, for which the mass and eccentricity only have upper limits, constrained by stability analyses (Lissauer et al. 2013). The estimated age of K11 is 9.7 ± 1.5 Gyr, and the star has a mass of ${M}_{* }=0.975\,\pm 0.031\,{M}_{\odot }$ and a radius of ${R}_{* }=1.193\pm 0.115\,{R}_{\odot }$ (Lissauer et al. 2013). The other system, K90 (Cabrera et al. 2014; Lissauer et al. 2014b; Schmitt et al. 2014), has seven planets orbiting an ${M}_{* }=1.2\pm 0.1\,{M}_{\odot }$ star.³ The K90 system (along with TRAPPIST-1) has the largest number of confirmed planets thus far, making it one of the closest in planet multiplicity to the solar system. However, unlike the solar system, all of K90's known planets are confined inside 1 au. The masses of K90's planets have not been measured directly, and as a result, the masses used here are estimated from the mass–radius relation described by Wright et al. (2011), which is based on the size distribution reported by Lissauer et al. (2011b). Table 1 summarizes the measured or estimated properties for K11 and K90 and Figure 1 emphasizes their relative spacing, mass, and inclination. The values in Table 1 are adopted for the present study, unless otherwise noted.

Zoom In Zoom Out Reset image size

Table 1.  Nominal Orbital Elements of the Known Planets of K11 and K90

Planet	Mass	a	e	i	ω
	(${M}_{\oplus }$)	(au)		(°)	(°)
Kepler-11	${M}_{* }=0.961\,{M}_{\odot }$		${R}_{* }=1.065\,{R}_{\odot }$
b	1.9	0.09	0.045	0.12	45.0
c	2.9	0.11	0.026	0.07	51.3
d	7.3	0.15	0.004	0.15	146.3
e	8.0	0.19	0.012	0.63	90.0
f	2.0	0.25	0.013	0.05	90.0
g	20.0	0.47	0.100	0.35	90.0
Kepler-90	${M}_{* }=1.2\,{M}_{\odot }$		${R}_{* }=1.2\,{R}_{\odot }$
b	2.4	0.076	⋯	0.28	⋯
c	1.7	0.088	⋯	0.00	⋯
d	7.9	0.307	⋯	0.03	⋯
e	6.9	0.424	⋯	0.11	⋯
f	8.1	0.520	⋯	0.09	⋯
g	69.1	0.736	⋯	0.12	⋯
h	297.9	0.996	⋯	0.08	⋯

Note. The stability of Kepler-11 is very sensitive to the values of the argument of pericenter for planets e, f, and g.

Download table as:  ASCIITypeset image

Published eccentricities and inclinations are used when available; otherwise, the orbital eccentricities are set to zero (e = 0), particularly for K90. The inclination values given in Table 1 are relative to a reference plane. The plane was determined by averaging the published orbital inclinations relative to a perpendicular to the line of sight to the system, which is $\langle {i}_{{\rm{K}}11}\rangle =89\buildrel{\circ}\over{.} 52$ and $\langle {i}_{{\rm{K}}90}\rangle =89\buildrel{\circ}\over{.} 68$. We create 100 realizations of each system, in which the longitudes and anomalies (Ω, ϖ and ${ \mathcal M }$) for each planet were drawn from a uniform random distribution between 0 and 2π. These realizations are first run without the external perturber to establish the stability of the systems over 10 Myr. As will be described in the results (Section 4), the stability of K11 is very sensitive to changes in the argument of pericenter. To further test the stability of K11 and K90 in the presence of a Jupiter analogue, we run the same initial conditions with an additional planet, also for 10 Myr. Unless otherwise noted, the perturber has a mass ${M}_{P}=1\,{M}_{J}$ and initial orbital elements ${a}_{P}=5.2\,\mathrm{au}$, e_P = 0.05, i_P = 13, ω_P = 2738, Ω_P = 1005 and ${{ \mathcal M }}_{P}=93\buildrel{\circ}\over{.} 8$, where the subindex P denotes parameters for the perturber. The realizations that include the perturber will be denoted as K11+ and K90+ throughout this work. The stability of the STIP is expected to be sensitive to the semimajor axis of the perturber (due to the secular frequencies), and as such, we systematically explore the semimajor axis parameter space of the perturber by forcing it to migrate inwards. The simulations that include the inward migration of the perturber are performed using the initial conditions that resulted in stable configurations for K11+ and K90+ after 10 Myr. For K11+ we used the full set of stable realizations, but chose only one of the K90+ stable systems. This choice was based on the analysis of the K11+ simulations, which all showed consistent behavior.

N-body integrations were run using a modified version of the Mercury6 code (Chambers 2012), which includes a general relativity correction due to the proximity of the planets to their host star. The correction used here takes the form

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{\mathrm{GR}}=-\displaystyle \frac{6\,{G}^{2}{M}^{2}}{{r}^{4}{c}^{2}}\ {\boldsymbol{r}},\end{eqnarray} \tag{ 1 }$$

which is appropriate for low-eccentricity orbits (Nobili & Roxburgh 1986). Because we want to resolve close encounters should they occur, we use the hybrid integrator of Mercury6. When in MVS mode, the time step is set to 10⁻² of the initial period of the innermost planet. Again, the total integration time, if not otherwise stated, is 10 Myr.

The inward migration of the Jupiter analogue is achieved by applying an acceleration term to the planet of the form

$$\begin{eqnarray}{{\boldsymbol{a}}}_{\mathrm{mig}} & = & -\displaystyle \frac{2\pi }{{\tau }_{\mathrm{mig}}}\left(\displaystyle \frac{1\,\mathrm{au}}{a}\right)[3(\hat{{\boldsymbol{r}}}\cdot \dot{{\boldsymbol{r}}})\hat{{\boldsymbol{r}}}+(\hat{{\boldsymbol{r}}}\times \dot{{\boldsymbol{r}}})\times \hat{{\boldsymbol{r}}}]\\ & & =-\displaystyle \frac{2\pi }{{\tau }_{\mathrm{mig}}}\left(\displaystyle \frac{1\,\mathrm{au}}{a}\right)[2(\hat{{\boldsymbol{r}}}\cdot \dot{{\boldsymbol{r}}})\hat{{\boldsymbol{r}}}+\dot{{\boldsymbol{r}}}],\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{r}}$ is the perturber's position vector from the star, $\dot{{\boldsymbol{r}}}$ its velocity vector, a its semimajor axis, and ${\tau }_{\mathrm{mig}}$ is the timescale for migration at 1 au. Over the region of interest, migration is nearly constant. For these simulations, we set ${\tau }_{\mathrm{mig}}=30\,\mathrm{Myr}$ which is equivalent to $\dot{a}\approx 0.42\,\mathrm{au}$ Myr⁻¹. During the 10 Myr migration, the orbital precession of the perturber increases by 1–2 orders of magnitude. The perturber's eccentricity also changes, decreasing among the stable systems from 0.05 to ∼0.035 at a_J = 3 au and to 0.02 at a_J = 1 au, corresponding to t = 6 and 10 Myr, respectively.

We complement the numerical simulations with secular theory as laid out by Murray & Dermott (1999). The theory considers only the secular contributions of the disturbing function to the equations of motion, expanded to second order in eccentricity, e, and inclination, i, as well as to first order in planetary mass. For systems in which the planets have small e and i, the theory can identify locations at which a test particle would be in an eccentricity or inclination secular resonance, although higher-order methods (e.g., Laskar 1985, 1986) or direct N-body calculations are needed to determine the outcome of a given resonance. If planets have strong interactions, such as what might be expected near mean motion resonances (MMR), then additional frequencies can become important and/or the predicted frequencies can become shifted with respect to the second-order theory. Because K11 and K90 exhibit TTVs, we do expect deviations from second-order theory, as is the case in the solar system with Jupiter and Saturn (Brouwer & van Woerkom 1950). Nevertheless, such effects are excluded in our calculations. We also exclude the effects of GR on the secular frequencies. Regardless, as will be shown in the results and discussion sections, secular theory seems to identify the locations of resonances with reasonable accuracy.

After building the second-order secular theory for any given planetary system configuration (e.g., the current STIP and perturber's orbital elements), we use the theory to highlight the resonant structure in the system. First, we introduce a test particle over a range of semimajor axes to sample the forced eccentricities and inclinations. Resonances will correspond to distances at which the forced eccentricity or inclination show very large increases, resulting from a singularity introduced when a precession frequency matches a system eigenfrequency. Using this approach, we can compare the resonant structure from secular theory with the behavior of the N-body simulations. For example, if a given planet becomes orbit crossing or develops large inclination variations, we can build the secular theory for the system, but exclude the highly perturbed planet(s). In K90(+), the two innermost planets are the most easily excited, so we examine the resonance structure by removing these planets from the secular theory, allowing us to see if the forced eccentricity and inclination at their semimajor axes suggest a resonance, at least for a test particle. This is only reasonable whenever the removed planets have low masses compared with the other planets. K90b and K90c are tightly coupled, so this must be done with some caution. Nevertheless, as we will show, this approach is reasonable enough to be useful for the present situation. As with K90(+), we examine the secular resonant structure in K11(+) by removing K11b, which is often the first planet to become orbit crossing, or by removing K11b and K11c as the planets are coupled. The eigenfrequencies for the second-order theory are obtained using a version of the secular theory python script available at https://github.com/norabolig/resmap.

We also analyze the resonant structure of K11(+) and K90(+) and examine the effects of strong planet interactions on the secular theory by calculating the synthetic frequency spectra for select systems. This is done by rerunning the given simulation for 0.3 Myr with output every 60 years. The synthetic spectra are found by taking the discrete Fourier transform of

$$\begin{eqnarray}&&e\,\exp (j\varpi )\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\sin (i)\,\exp (j{\rm{\Omega }}),\end{eqnarray} \tag{ 4 }$$

for each planet, where e is the eccentricity, i the inclination, while ϖ and Ω are the longitudes of pericenter and ascending node, respectively, and $j=\sqrt{-1}$. The Fourier transform of Equations (3) and (4) return the principal apsidal and nodal precession frequencies and amplitudes that compose the synthetic secular theory.

The N-body simulations of K11 and K90 (without the perturber) allow us to examine the dynamical stability of our realizations of these systems. We find that the stability of K11 is very sensitive to perturbations to the nominal argument of pericenter whenever we consider non-zero eccentricity. In particular, 44% of the simulations become unstable when using the nominal eccentricities but randomizing the longitudes. This is qualitatively consistent with Mahajan & Wu (2014), who found that K11 is preferentially unstable if any of the eccentricities e > 0.04, using the same planetary masses used here. In contrast, 100% of the realizations in this study are stable if the initial eccentricity is zero, which agrees with the zero-eccentricity stability tests of Lissauer et al. (2011a, 2013). For K90, 20% of the our realizations present instability despite having initial orbital eccentricities set to zero for all planets. This result is unexpected at face value, as we might anticipate the system to be stable (as with K11) under such conditions. This agrees with the stability test of K90's discovery paper (Cabrera et al. 2014). The system's secular perturbations could be the origin of the noted instability and/or MMR, as will be discussed later.

The simulations of K11 and K90 can now be used as a reference for exploring additional dynamical perturbations. Using the same realizations for K11 and K90, we now include a Jupiter analogue in each system as described in Section 3. We find that the stability of the STIPs is not affected by the presence of the single perturber when it is placed at 5.2 au. The same 44% and 20% analogues become unstable in less than 10 Myr for K11+ and K90+, respectively, which suggest that the instability is driven solely by the inner planets, at least for this short time period and for the initial placement of the Jupiter analogue.

While the stability of the STIPs is unaffected, the systems do respond to the presence of the outer perturber. This is highlighted in Figures 2 and 3, which show the inclination evolution for the planets in stable realizations of K11 and K90, respectively, with and without the presence of the Jupiter analogue. In each case, the planets exhibit variation among their inclinations with time. When the perturber is included, however, there is an additional large-scale variation in all of the orbital inclinations, i.e., each planet's inclination oscillates with respect to a common orbital plane while the orientation of the shared orbital plane changes. The presence of the perturber precesses the longitudes of ascending node of all the inner planets at the same rate, which is also shown in Figures 2 and 3 (with ΔΩ). Such "dynamical rigidity" has been seen in simulations of other systems when perturbers are introduced (Kaib et al. 2011) at high inclination, and has been further explored analytically (Boué & Fabrycky 2014a, 2014b). Several test simulations were run with the Jupiter analogue at high inclination (50°), which showed that the longitudes of ascending node remained locked, with the STIP undergoing large and coherent inclination variations (Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have thus far only considered one location for a perturber, arbitrarily introduced at 5.2 au. We now explore a wider range of semimajor axes by forcing the perturber to migrate inwards from 5.2 au to ∼1 au. First, the overall behavior of the STIPs is unchanged by simply moving the Jupiter analogue inwards. Figures 5 and 6, for example, show that the systems continue to exhibit coherent changes in the orbital plane, but with an increasing precession rate of that plane as the Jupiter analogue approaches the inner system. However, for certain system configurations, the outcomes can be very different. In most cases, this appears to be due to shifts in the secular frequencies. For example, eccentricity resonances can drive the STIP toward instability, causing planetary collisions (as occurs for K11). An inclination resonance can force a planet (or planets) out of the original common orbital plane, effectively creating multiple orbital planes in a stable system (as occurs in K90). Specifically, runs with an initially stable K11+ develop an instability if the perturber reaches a semimajor axis of ∼1.0 au, which places the perturber in a 3:1 near-MMR with K11g. This commensurability could be causing the instability in K11+, but the secular resonances could also have a strong contribution to the destabilization of K11+, as will be discussed shortly. At this time, 59% of the unstable K11+ realizations result in either K11b crossing the orbit of K11c followed by K11f crossing K11e's orbit or vice versa.⁴ In K90+, when a_J ≲ 3.2 au, K90b and K90c evolve together away from the rest of the planets by about 16° in inclination, with the system remaining stable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The behavior of both K11+ (Figures 7 and 8) and K90+ (Figure 9) can be understood by looking at their forced eccentricities and inclinations. As the Jupiter analogue moves inwards, the locations of the innermost inclination and eccentricity resonances move outwards in the STIP, eventually crossing the innermost planets. For K11+, an inclination and eccentricity resonance overlap the semimajor axis of K11b when a_J = 0.98 au, which is when instability occurs in the corresponding simulations (Figure 7). For this analysis, K11b and K11c were removed from the secular theory calculation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of K90+, when the perturber is at a_J ≈ 3.3 au the two innermost planets appear to be trapped in an inclination secular resonance as determined by secular theory (with planets K90b and K90c removed, Figure 9). Together, Figures 6 and 9 suggest that the separation of the K90+ system into two distinct orbital planes (for a_J ≲ 3.2 au) is due to an inclination resonance overlapping the K90b and K90c positions. This secular resonance excites the inclination of both planets, which are strongly coupled, and its strength increases as the perturbing planet's distance to K90 decreases. This second orbital plane for K90b and K90c can acquire a maximum inclination i_bc ∼ 18° if a_J ∼ 3.0 au. While there is also an inclination overlapping the location of K11b in K11+ (when a_j ≈ 0.98 au), there is also an eccentricity resonance present (Figure 7) possibly leading to instability before any large inclination changes could occur.

As already discussed, removing planets b and c from both K11+ and K90+ allows us to explore whether there are forced inclination and eccentricities at their locations. This assumes that the planets are massless, which is not true. As such, the method is not guaranteed to highlight resonances, although the correspondence with the N-body simulations suggest the approximation is valid in this case. The calculation nonetheless reduces the complexity of the secular structure by removing frequencies. To highlight this, we show in Figure 8 a secular map in which we only remove K11b, allowing K11c to contribute to the secular model. An inclination and eccentricity resonance is present at r ≈ 0.09 au, just inside K11b. The position of this resonance is fixed regardless of the perturber's proximity to the inner system, suggesting that the location of this resonance is due to the K11+ inner planets. There are additional resonances that are at even smaller semimajor axes, which are affected by the perturber's location. As the perturber moves inwards, this inner resonance structure moves outwards and the resonance wings overlap. At face value, Figure 8 suggests that K11b might not overlap a secular resonance. We will later show that K11b and K11c are strongly coupled and Figure 7 might better reflect the outer system's influence on K11b and K11c. Nevertheless, Figure 8 shows the proximity of these resonances to K11b, which could explain the observed instability of the system and its apparent fine tuning even when the perturber is absent. A small perturbation in the configuration of planets could cause K11b to enter a secular eccentricity resonance and collide with K11c, the most common outcome in our simulations that become unstable.

We investigate the secular structure of the systems further by calculating synthetic secular frequencies (precession spectra) with and without the perturber. The precession spectra are determined directly from the N-body calculations, using 0.3 Myr of output with a sample time of 60 years. In K11+ and K90+, the perturber was set to a_J = 0.98 au and a_J = 3.0 au, respectively. The comparison of the precession spectra with apsidal and nodal secular eigenfrequencies is shown in Figures 10 and 11. As expected, the value and number of eigenfrequencies change from K11 to K11+ and K90 to K90+, respectively, due to the additional planet in the perturbed systems, with a greater displacement for the apsidal eigenfrequencies. The amplitude of the global Fourier spectra, in both eccentricity and inclination, for K11+ and K90+ is higher than in the unperturbed analogues; there is about one order of magnitude difference in K11+ and two to three in K90+. The nodal eigenfrequencies show very good agreement with the synthetic spectra for both K11(+) and K90(+), although there are some deviations. The apsidal precession frequencies of both systems exhibit stronger disagreement between the synthetic spectra and the secular theory, but the features in the Fourier spectra remain recognizable, particularly for K11(+). The misalignment of the principal apsidal frequencies relative to the estimated secular eigenfrequencies might indicate that there are strong interactions among the planets, which are neglected in the secular calculation. This could include nonlinear secular resonances and near MMR. Malhotra et al. (1989) showed that the e − ϖ eigenfrequencies are significantly shifted if a near first order MMR is present in the system because the averaged effect of the term $e\,\cos (j\lambda -(j+1)\lambda ^{\prime} +\varpi )$ in the secular expansion is non-zero. There is a weaker effect of such a near-resonant term on i − Ω because the nodal term is associated with i². The MMR can also add eigenfrequencies to the system, as occurs in the solar system due to the near-MMR between Jupiter and Saturn (Brouwer & van Woerkom 1950).

Zoom In Zoom Out Reset image size

The main precession period of each planet, in K11(+) and K90(+), and the corresponding amplitude are displayed in Table 2. In the four systems, the two innermost planets have the same main precession period ($\dot{\varpi }$ and $\dot{{\rm{\Omega }}}$) and their precession spectra are almost indistinguishable, confirming the strong coupling of these planetary pairs. Furthermore, in K11, the nodes of the six planets precesses with the same period, suggesting that K11g is dominating the nodal precession, which causes a dynamical rigidity in the orbital plane. On the other hand, there are two distinct pericenter precession periods, in which K11b and K11c have an apsidal period of ${P}_{{\rm{K}}11{\rm{b}} \mbox{-} {\rm{c}}}\simeq 1346\,\mathrm{years}$ and K11d through K11g have periods of ${P}_{{\rm{K}}11{\rm{d}} \mbox{-} {\rm{g}}}\simeq {\rm{33,340}}\,\mathrm{years}$. The perturber in K11+ changes the frequencies and decouples K11g; K11g's node precesses at the same rate as the perturber, which is ∼4 times faster than the rest of planets in K11+. In contrast to K11, the main nodal precession periods of K90 do not suggest an intrinsic dynamical rigidity, although there is strong coupling between K90b and K90c, as well as K90g and K90h. However in K90+, there are two distinct nodal precession periods that are different by only 7%. These nearly overlapping eigenfrequencies (Figure 11) suggest a secular resonance is causing the large change in inclinations for K90b and K90c (Figure 6). The amplitude of the broad frequency peak seen in Figure 11 is dominant by one order of magnitude or more. When the K90+ simulation is allowed to evolve to 0.4 Myr, it is possible to identify the two close frequency components in the dominant broad peak of Figure 11, as expected from secular theory for K90+. These are shown in Figure 12, with the peaks corresponding to ${P}_{\mathrm{node}}\simeq {\rm{23,100}}$ years and ${P}_{\mathrm{node}}\simeq {\rm{21,400}}\,\mathrm{years}$, which straddle the eigenfrequency associated with K90 in that range. The amplitude of these two frequencies is interchanged when the spectra is calculated from a subsample. The whole precession spectra of the longer N-body integration is not shown in Figure 12 due to severe aliasing at large frequencies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The identification of the frequency components in the apsidal precession spectra of K90 and K90+ (top panels in Figure 11) is also complicated due to aliasing. Even so, we can discern a substantial misalignment of the main frequency components compared with the nearest eigenfrequency, as mentioned above. In K90, there is a double-peaked feature corresponding to K90b and K90c at $\dot{\varpi }\sim 0\buildrel{\circ}\over{.} 05$ yr⁻¹. This is absent in the spectra of the other planets, which could originate from a nonlinear secular resonance or from a MMR.

The discrepancies between the frequency components of the precession spectra of K11 and K90 and their secular eigenfrequencies demonstrate that either additional physics should be included in the secular theory or that we need to consider a higher order expansion in e, i, and mass. Some of the physical effects that we know we are neglecting in the secular code include GR and MMRs. The GR contribution to the precession rates is small; it only affects the apsidal precession frequencies of the synthetic spectra by ≲2 %, with K90 being the most affected. GR does not shift the nodal spectra in either system. The MMR, depending on the order, can affect linear terms in eccentricity and inclination (first order).

K11 and K90 exhibit TTVs (Lissauer et al. 2011a; Cabrera et al. 2014), which occur if there are strong gravitational interactions between the planets. The discovery papers of both K11 and K90 discussed MMRs for different planet combinations. The period ratios of K11b–K11c, and K90b–K90c suggest a 5:4 near commensurability in both systems. Additionally, Cabrera et al. (2014) discusses a possible 2:3:4 near resonance between K90d, K90e and K90f. These first claims of MMR in the system were based on the period ratios of neighboring planets, but dynamically a MMR might not occur, due to the dependence on the eccentricity and inclination of the planets. In order to confirm the existence of a MMR in a pair of planets it is necessary to estimate their resonant angle. A MMR is present if the resonant angle librates near a constant value.

We investigate first-order MMRs between neighboring planets by calculating their resonant angle from the N-body simulations directly using

$$\begin{eqnarray}&&\varphi ^{\prime} =j\lambda ^{\prime} +(1-j)\lambda -\varpi ^{\prime} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\varphi =j\lambda ^{\prime} +(1-j)\lambda -\varpi ,\end{eqnarray} \tag{ 6 }$$

where λ is the mean longitude and ϖ is the longitude of pericenter. The apostrophe denotes the external planet. Using the given masses and orbital elements for K11, all of the resonant angle combinations circulate at all times. Thus, we cannot confirm a 5:4 near-MMR between K11b and K11c. In K90's case, two pairs of neighboring planets display a librating resonant angle (Figure 13). These plots suggest that K90b is in 5:4 near-MMR with K90c and that K90g is in a 3:2 near-MMR with K90h. In both cases, the resonant angle librates with large-amplitude around 0° and 180°. The existence of the 5:4 between K90b-c suggests that some migration could have taken place. The possible 3-body resonance between d, e, and f noted by Cabrera et al. (2014) would correspond to a resonant angle $\varphi =2{\lambda }_{d}-6{\lambda }_{e}+4{\lambda }_{f}$. After plotting this angle from our simulations, we find that the resonant angle circulates through the length of the realization. The 3:2 near-MMR between K90g and K90h could be the origin of the instability in the K90 realizations. The instability in K90 depends on the exact masses of the planets, particularly those for K90g and K90h (Cabrera et al. 2014), which are the most massive. Although, the nominal masses of K90 planets are unknown, the impact of these massive planets on the system might be strengthened by their apparent MMR.

Zoom In Zoom Out Reset image size

The behavior of the K11 and K90 analogues explored here show a range of dynamical outcomes that STIPs could have in the presence of an outer (undetected) perturber. To date, it is uncertain whether the detected planets among the known exoplanet systems represent a complete set or whether there are undetected planets in those systems. Follow up measurements of orbital inclinations of STIPs could reveal dynamic rigidity in the orbital plane and as a consequence, could be used as an indirect detection method for an outer planet. The observational determination of the longitudes of the nodes of planets in a given STIP could tell us whether there is an external perturber should the Ωs for all the planets be ≲20°. In a non-perturbed STIP, the nodes would not preferentially align to a value, although it might be the case that a subset of the planets present nodal alignment due to coupling. Moreover, an external perturber will cause a coherent change in the inclination of the planets, causing the duration of their transits to change coherently. An outer planet could also cause large breaks in the nominal orbital plane, which could explain some spin–orbit misalignment of low mass planets or reduce the number detected by transit.

In this paper, we have only considered a single Jovian perturber as a first step. Future work will need to include the effects of multiplicity among outer planetary systems. Having two or more Jovian planets would modify the orbital precession frequencies and increase the complexity of the system' secular frequency spectrum, particularly if the giant planets are mutually interacting. This would have consequences for the dynamics and stability of STIPs and could increase the fraction of inner systems that become unstable. As an example, consider studies based, at least in part, on the solar system, which show that the stability of the inner solar system or inner system analogues can depend sensitively on the orbital configuration of the solar system's outer planets (e.g., Lithwick & Wu 2011; Agnor & Lin 2012; Clement & Kaib 2017).

In our simulations, the Jovian perturber was placed on a low eccentricity orbit (e ≤ 0.05). If the perturber's orbital eccentricity were to be increased, which affects the width and power of secular resonances, we would expect additional instability. For example, Clement & Kaib (2017) showed that increasing the nominal eccentricity of Jupiter and Saturn by a factor of two enhances chaos within the inner solar system and reduces the system's stability.

Finally, we integrated our simulations for 10 Myr, and used a corresponding "migration" timescale for the perturber. This timescale is a compromise between integrating these systems with sufficiently small time-steps and evolving the systems for a dynamically meaningful duration. In either set of simulations (with or without a perturber), we expect the fraction of unstable systems to increase with time (Volk & Gladman 2015). While we find that the presence of an outer planet does not typically affect the stability of the planetary systems studied here, the migration timescale used to sample different system configurations could have caused the perturber to move through an unstable configuration too quickly, not allowing the instability to develop. As such, the fractional difference between the number of systems that do become unstable with and without an outer planet could increase for longer simulations and should be explored in future work.

5.1. Comparison with Previous and Ongoing Studies

In this paper, we studied the stability and observability of two known high-multiplicity STIPs, K11 and K90, in the presence of an outer Jupiter-like planet through N-body simulations and secular theory. The presence of the perturber causes dynamical rigidity about a common orbital plane for the inner planets, while the stability of the system remains unaltered when compared with unperturbed realizations for most perturber locations. The observed instability seems to be inherent to STIPs, suggesting secular resonances among the planets. The rigid behavior of the orbital plane occurred for most of the parameter space that we explored as long as no instability developed. The presence of the perturber also caused two possible effects on systems that are otherwise stable: (1) the orbital plane of the planets could be separated into two distinct planes, as in K90+, and (2) the system could become unstable for particular perturber locations. The N-body simulations and secular analysis demonstrate that the instability and multiple orbital planes are consequences of the eccentricity and inclination secular resonances, respectively. For K11, we suggest that the eccentricity resonance close to K11b is the source of the system's inherent instability.

Comparing STIP secular eigenfrequencies to the synthetic counterparts provides a deeper insight into the coupling and possible presence of MMR between planets. K11's nodal precession frequencies indicate dynamical rigidity, seemingly due to K11g, although this is dependent on the actual mass of K11g. In K90, our simulations show a 5:4 (near) MMR between K90b and K90c, as well as a 3:2 (near) MMR between K90g and K90h.

Observations of the rigid behavior of a STIP would indicate the existence of an outer planetary system. It is possible that some of the detected planetary systems with low multiplicities are part of higher multiplicity STIPs that are affected by secular resonances.

We thank B. Gladman, C. van Laerhoven, and the anonymous referee for comments that improved this manuscript. This work was supported by an NSERC Discovery grant, The University of British Columbia, the Canadian Foundation for Innovation, and the BC Knowledge Development Fund. The simulations presented here were enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). This research 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.