Prevalance of Chaos in Planetary Systems Formed Through Embryo Accretion
Matthew S. Clement, Nathan A. Kaib
TL;DR
This study uses numerical simulations to analyze chaos in planetary systems formed through embryo accretion, revealing that chaos is common and mainly driven by Jupiter's perturbations, with implications for understanding our solar system's dynamics.
Contribution
It provides a detailed characterization of chaos in simulated terrestrial planet systems, highlighting Jupiter's role and the prevalence of late instabilities.
Findings
Chaos occurs in about half of the simulated systems.
Removing Jupiter reduces chaos significantly.
Chaos is more common in systems with 4-5 planets.
Abstract
The formation of the solar system's terrestrial planets has been numerically modeled in various works, and many other studies have been devoted to characterizing our modern planets' chaotic dynamical state. However, it is still not known whether our planet's fragile chaotic state is an expected outcome of terrestrial planet accretion. We use a suite of numerical simulations to present a detailed analysis and characterization of the dynamical chaos in 145 different systems produced via terrestrial planet formation in Kaib and Cowan (2015). These systems were created in the presence of a fully formed Jupiter and Saturn, using a variety of different initial conditions. They are not meant to provide a detailed replication of the actual present solar system, but rather serve as a sample of similar systems for comparison and analysis. We find that dynamical chaos is prevalent in roughly half…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20| Run | Run Time | Jupiter | Saturn | Outer Planet |
|---|---|---|---|---|
| Name | eccentricity shift | |||
| 1 | 100 Myr | unchanged | ||
| a | 100 Myr | removed | removed | unchanged |
| b | 100 Myr | removed | unchanged | |
| c | 100 Myr | 150% | ||
| d | 100 Myr | 200% | ||
| e | 100 Myr/1 Gyr | + 0.05 |
| system | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| ann | 71% | 12% | 2% | 12% | 0% | 4% |
| cjs | 48% | 18% | 18% | 14% | 2% | 0% |
| body | cjs10 a | cjs13 a | sol a | cjs10 m | cjs13 m | sol m |
|---|---|---|---|---|---|---|
| (AU) | (AU) | (AU) | () | () | () | |
| 1 | 0.51 | 0.54 | 0.39 | 0.61 | 0.68 | 0.055 |
| 2 | 0.76 | 0.78 | 0.72 | 0.63 | 0.18 | 0.82 |
| 3 | 1.1 | 1.1 | 1.0 | 1.0 | 0.77 | 1.0 |
| 4 | 1.9 | 1.6 | 1.5 | 0.46 | 0.59 | 0.11 |
| 5 | 2.5 | 2.7 | 2.8 | 0.32 | 0.80 | 0.00015 |
| System | MEGNO | Inner | Run 1111“Y” indicates chaos was detected, “N” indicates it was not. | Run a111“Y” indicates chaos was detected, “N” indicates it was not. | Run b111“Y” indicates chaos was detected, “N” indicates it was not. | Run 1 | |
|---|---|---|---|---|---|---|---|
| (years) | Planets | MMRs222(J) and (S) indicate the resonance was with Jupiter or Saturn respectively. | |||||
| cjs1 | 100.7 | 5 | Y | Y | Y | ||
| cjs2 | 45.35 | 4 | Y | N | N | 7:4 | |
| cjs3 | 187.2 | 5 | Y | Y | Y | ||
| cjs4 | 160.9 | 5 | Y | Y | Y | 2:1 | |
| cjs5 | 183.5 | 5 | Y | Y | Y | ||
| cjs6 | 160.0 | 5 | Y | Y | Y | ||
| cjs7 | 126.0 | 5 | Y | N | Y | 7:4,5:1, | |
| 7:1 (S), | |||||||
| 5:2 (J,S) | |||||||
| cjs8 | 2.004 | 4 | N | N | N | ||
| cjs9 | 1.997 | 3 | N | N | N | ||
| cjs10 | 1.896 | 5 | N | N | N | ||
| cjs11 | 1.999 | 2 | N | N | N | ||
| cjs12 | 90.97 | 5 | Y | Y | Y | ||
| cjs13 | 2.003 | 5 | N | N | N | ||
| cjs14 | 1.988 | 3 | N | N | N | ||
| cjs15 | 8.794 | 5 | Y | Y | Y | 7:3 (J) | |
| cjs16 | 212.8 | 4 | Y | Y | Y | ||
| cjs17 | 4.467 | 4 | Y | N | N | 3:1 (J) | |
| cjs18 | 21.04 | 3 | Y | N | N | 3:1 (J) | |
| cjs19 | 53.36 | 4 | Y | Y | Y | ||
| cjs20 | 2.000 | 2 | N | N | N | ||
| cjs21 | 1.996 | 2 | N | N | N | ||
| cjs22 | 2.000 | 3 | N | N | N | ||
| cjs23 | 169.0 | 4 | Y | Y | Y | ||
| cjs24 | 20.29 | 4 | Y | N | Y | 5:2 (J,S) | |
| cjs25 | 2.040 | 4 | N | N | N | 2:1 | |
| cjs26 | 183.2 | 5 | Y | Y | Y | 6:1 | |
| cjs27 | 1.999 | 3 | N | N | N | ||
| cjs28 | 1.999 | 2 | N | N | N | ||
| cjs29 | 12.82 | 4 | Y | Y | Y | 5:3:1 | |
| cjs30 | 4.816 | 3 | Y | N | N | ||
| cjs31 | 2.001 | 2 | N | N | N | ||
| cjs32 | 1.966 | 4 | N | N | N | ||
| cjs33 | 1.999 | 2 | N | N | N | ||
| cjs34 | 117.0 | 4 | Y | Y | Y | ||
| cjs35 | 2.000 | 3 | N | N | N | ||
| cjs36 | 186.5 | 4 | Y | N | Y | ||
| cjs37 | 193.5 | 5 | Y | Y | Y | ||
| cjs38 | 5.732 | 3 | Y | N | Y | ||
| cjs39 | 1.999 | 4 | N | N | N | ||
| cjs40 | 1.999 | 4 | N | N | N | ||
| cjs41 | 113.5 | 4 | Y | Y | Y | 5:3 | |
| cjs42 | 1.998 | 4 | N | N | N | 5:1 | |
| cjs43 | 49.51 | 4 | Y | N | Y | 4:1 | |
| cjs44 | 189.4 | 4 | Y | Y | Y | ||
| cjs45 | 194.1 | 4 | Y | N | Y | ||
| cjs46 | 2.011 | 4 | N | N | N | ||
| cjs47 | 129.4 | 4 | Y | N | Y | ||
| cjs48 | 8.483 | 4 | Y | Y | Y | ||
| cjs49 | 20.44 | 4 | Y | N | Y | ||
| cjs50 | 6.271 | 4 | Y | N | Y | ||
| eejs1 | 178.4 | 3 | Y | N | Y | ||
| eejs2 | 87.08 | 4 | Y | Y | Y | 5:3 | |
| eejs3 | 202.9 | 4 | Y | N | Y | 8:5,7:1 | |
| eejs4 | 174.4 | 1 | Y | N | N | ||
| eejs5 | 2.000 | 2 | N | N | N | ||
| eejs6 | 29.38 | 3 | Y | Y | Y | 3:1,7:1 (S) | |
| eejs7 | 34.65 | 4 | Y | N | Y | ||
| eejs8 | 1.988 | 3 | N | N | N | ||
| eejs9 | 9.676 | 5 | Y | Y | Y | 8:3,8:5 | |
| eejs10 | 151.8 | 1 | Y | N | N | ||
| eejs11 | 2.000 | 2 | N | N | N | ||
| eejs13 | 7.045 | 4 | Y | Y | N | ||
| eejs15 | 1.999 | 1 | N | N | N | ||
| eejs16 | 2.002 | 3 | N | N | N | 5:2 (J,S) | |
| eejs18 | 2.017 | 4 | N | N | N | ||
| eejs19 | 2.000 | 1 | N | N | N | ||
| eejs20 | 1.992 | 3 | N | N | N | ||
| eejs21 | 1.919 | 3 | N | N | N | ||
| eejs22 | 1.998 | 2 | N | N | N | ||
| eejs23 | 183.4 | 4 | Y | Y | Y | ||
| eejs24 | 57.05 | 5 | Y | Y | Y | 8:1 (S) | |
| eejs25 | 2.041 | 2 | N | N | N | ||
| eejs26 | 1.912 | 3 | N | N | N | ||
| eejs27 | 1.813 | 3 | N | N | N | ||
| eejs28 | 4.370 | 3 | Y | N | N | ||
| eejs29 | 173.0 | 1 | Y | N | N | 5:2 (J,S) | |
| eejs30 | 158.7 | 2 | Y | N | Y | ||
| eejs31 | 184.6 | 4 | Y | Y | Y | ||
| eejs32 | 1.969 | 3 | N | N | N | 7:3 | |
| eejs33 | 2.004 | 3 | N | N | N | 7:2,8:5 | |
| eejs34 | 196.5 | 1 | Y | Y | N | ||
| eejs35 | 154.0 | 2 | Y | N | Y | 7:3 | |
| eejs36 | 180.2 | 1 | Y | N | N | ||
| eejs37 | 183.2 | 4 | Y | Y | Y | ||
| eejs38 | 171.8 | 2 | Y | N | N | ||
| eejs39 | 2.005 | 2 | N | N | N | ||
| eejs40 | 189.1 | 4 | Y | Y | Y | ||
| eejs41 | 33.10 | 4 | Y | Y | Y | 5:3 | |
| eejs42 | 2.511 | 4 | N | N | N | ||
| eejs43 | 1.998 | 1 | N | N | N | ||
| eejs44 | 193.8 | 3 | Y | N | Y | ||
| eejs45 | 178.1 | 4 | Y | N | Y | ||
| eejs46 | 2.009 | 4 | N | N | N | ||
| eejs47 | 8.082 | 3 | Y | N | N | 8:5 | |
| eejs49 | 1.996 | 3 | N | N | N | 5:2 | |
| eejs50 | 1.929 | 2 | N | N | N | 7:3 | |
| ann1 | 2.014 | 3 | N | N | Y | ||
| ann2 | 1.998 | 3 | N | N | N | ||
| ann3 | 1.988 | 4 | N | N | N | ||
| ann4 | 2.000 | 3 | N | N | N | ||
| ann5 | 29.32 | 3 | Y | N | Y | 8:5,7:1 (J) | |
| ann6 | 7.441 | 4 | Y | N | N | 7:4 | |
| ann7 | 3.132 | 3 | Y | N | N | ||
| ann8 | 185.6 | 3 | Y | N | Y | ||
| ann9 | 2.010 | 3 | N | N | N | 2:1 | |
| ann10 | 1.998 | 2 | N | N | N | ||
| ann11 | 2.294 | 3 | N | N | N | ||
| ann12 | 195.1 | 4 | Y | Y | Y | ||
| ann13 | 125.7 | 2 | Y | Y | Y | ||
| ann14 | 60.09 | 4 | Y | N | N | 4:1 | |
| ann15 | 2.103 | 3 | N | N | N | ||
| ann16 | 1.99 | 2 | N | N | N | ||
| ann17 | 4.372 | 4 | Y | N | N | 5:2 (J,S) | |
| ann18 | 2.159 | 3 | N | N | N | ||
| ann19 | 183.8 | 3 | Y | Y | Y | ||
| ann20 | 175.4 | 4 | Y | Y | Y | ||
| ann21 | 21.76 | 3 | Y | N | Y | 3:1 | |
| ann22 | 44.64 | 3 | Y | N | Y | ||
| ann23 | 53.55 | 3 | Y | N | Y | 7:4 | |
| ann24 | 1.999 | 3 | N | N | N | ||
| ann25 | 20.22 | 4 | Y | Y | Y | 8:3 | |
| ann26 | 9.212 | 4 | Y | Y | Y | 5:3 | |
| ann27 | 13.65 | 4 | Y | Y | N | 5:2,8:5 | |
| ann28 | 15.54 | 5 | Y | Y | Y | ||
| ann29 | 115.7 | 4 | Y | N | Y | ||
| ann30 | 193.1 | 3 | Y | Y | Y | ||
| ann31 | 8.47 | 3 | Y | Y | Y | ||
| ann32 | 182.3 | 4 | Y | Y | Y | ||
| ann33 | 191.6 | 4 | Y | Y | Y | 2:1 | |
| ann35 | 23.00 | 3 | Y | N | N | ||
| ann36 | 1.998 | 3 | N | Y | N | 7:3 | |
| ann37 | 192.5 | 4 | Y | Y | Y | ||
| ann38 | 2.001 | 3 | N | N | N | ||
| ann39 | 189.6 | 3 | Y | Y | Y | ||
| ann40 | 192.9 | 4 | Y | Y | Y | ||
| ann41 | 1.998 | 3 | N | N | N | ||
| ann42 | 2.000 | 2 | N | N | N | ||
| ann43 | 2.000 | 3 | N | N | N | ||
| ann44 | 11.41 | 3 | Y | N | N | ||
| ann45 | 2.000 | 2 | N | N | N | ||
| ann46 | 152.4 | 3 | Y | Y | Y | ||
| ann47 | 2.004 | 3 | N | N | N | ||
| ann48 | 54.18 | 4 | Y | Y | Y | 5:2 (J,S) | |
| ann49 | 185.7 | 3 | Y | N | Y | 7:4 | |
| ann50 | 2.008 | 3 | N | N | N |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Prevalence of Chaos in Planetary Systems Formed Through Embryo Accretion
Matthew S. Clement11affiliation: HL Dodge Department of Physics Astronomy, University of Oklahoma, Norman, OK 73019, USA & corresponding author email: [email protected] & Nathan A. Kaib11affiliation: HL Dodge Department of Physics Astronomy, University of Oklahoma, Norman, OK 73019, USA & corresponding author email: [email protected]
Abstract
The formation of the solar system’s terrestrial planets has been numerically modeled in various works, and many other studies have been devoted to characterizing our modern planets’ chaotic dynamical state. However, it is still not known whether our planets’ fragile chaotic state is an expected outcome of terrestrial planet accretion. We use a suite of numerical simulations to present a detailed analysis and characterization of the dynamical chaos in 145 different systems produced via terrestrial planet formation in Kaib & Cowan (2015). These systems were created in the presence of a fully formed Jupiter and Saturn, using a variety of different initial conditions. They are not meant to provide a detailed replication of the actual present solar system, but rather serve as a sample of similar systems for comparison and analysis. We find that dynamical chaos is prevalent in roughly half of the systems we form. We show that this chaos disappears in the majority of such systems when Jupiter is removed, implying that the largest source of chaos is perturbations from Jupiter. Chaos is most prevalent in systems that form 4 or 5 terrestrial planets. Additionally, an eccentric Jupiter and Saturn is shown to enhance the prevalence of chaos in systems. Furthermore, systems in our sample with a center of mass highly concentrated between 0.8–1.2 AU generally prove to be less chaotic than systems with more exotic mass distributions. Through the process of evolving systems to the current epoch, we show that late instabilities are quite common in our systems. Of greatest interest, many of the sources of chaos observed in our own solar system (such as the secularly driven chaos between Mercury and Jupiter) are shown to be common outcomes of terrestrial planetary formation. Thus, consistent with previous studies such as Laskar (1996), the solar system’s marginally stable, chaotic state may naturally arise from the process of terrestrial planet formation. Keywords: Chaos, Planetary Formation, Terrestrial Planets Received; Accepted
1. Introduction
Our four terrestrial planets are in a curious state where they are evolving chaotically, and are only marginally stable over time (Laskar, 1996, 2008; Laskar & Gastineau, 2009). This chaos is largely driven by interactions with the 4 giant planets. However our understanding of the dynamical evolution of the gas giants, particularly Jupiter and Saturn, has changed drastically since the introduction of the Nice Model (Gomes et al., 2005; Morbidelli et al., 2005; Tsiganis et al., 2005).
The classical model of terrestrial planetary formation, where planets form from a large number of small embryos and planetesimals that interact and slowly accrete, is the basis for numerous studies of planetary evolution (e.g. Chambers, 2001a; O’Brien et al., 2006; Chambers, 2007; Raymond et al., 2009b; Kaib & Cowan, 2015). Using direct observations of proto-stellar disks (Currie et al., 2009), it is clear that free gas disappears long before the epoch when Earth’s isotope record indicates the conclusion of terrestrial planetary formation (Halliday, 2008). For these reasons, a common initial condition taken when numerically forming the inner planets is a fully formed system of gas giants at their current orbital locations. Many numerical models have produced planets using this method. However, none to date have analyzed the chaotic nature of fully evolved accreted terrestrial planets up to the solar system’s current epoch. It should be noted that other works have modeled the outcome of terrestrial planetary formation up to 4.5 Gyr. Laskar (2000) evolved 5000 such systems from 10000 planetesimals and showed correlations between the resulting power-law orbital spacing and the initial mass distribution. Furthermore, many works have performed integrations of the current solar system, finding solutions that showed both chaos and a very real possibility of future instabilities (Laskar, 2008; Laskar & Gastineau, 2009). Our work is unique in that we take systems formed via direct numerical integration of planetary accretion, evolve them to the solar system’s age, probe for chaos and its source, and draw parallels to the actual solar system.
Although the classical terrestrial planet formation model has succeeded in replicating many of the inner solar system’s features, the mass of Mars remains largely unexplained (Chambers, 2001a; O’Brien et al., 2006; Chambers, 2007; Raymond et al., 2009b; Kaib & Cowan, 2015). Known as the Mars mass deficit problem, most simulations routinely produce Mars analogues which are too massive by about an order of magnitude. Walsh et al. (2011a) argue for an early inward, and subsequent outward migration of a fully formed Jupiter, which results in a truncation of the proto-planetary disc at 1 AU prior to terrestrial planetary formation. If correct, this “Grand Tack Model” would explain the peculiar mass distribution observed in our inner solar system. Another interesting solution involves local depletion of the disc in the vicinity of Mars’s orbit (Izidoro et al., 2014). A detailed investigation of the Mars mass deficit problem is beyond the scope of this paper. It is important, however, to note that accurately reproducing the mass ratios of the terrestrial planets is a significant constraint for any successful numerical model of planetary formation.
Through dynamical modeling, we know chaos is prevalent in our solar system (Sussman & Wisdom, 1988; Laskar, 1989; Sussman & Wisdom, 1992; Laskar, 2008). It is important to note the difference between “stability” and “chaos.” While a system without “chaos” can generally be considered stable, a system with “chaos” is not necessarily unstable (Milani & Nobili, 1992; Deck et al., 2013). As is convention in other works, in this paper “chaos” implies both a strong sensitivity of outcomes to specific initial conditions, and a high degree of mixing across all energetically accessible points in phase space (Deck et al., 2013) Conversely “Instability” is used to describe systems which experience specific dynamical effects such as ejections, collisions or excited eccentricities.
The chaos in our solar system mostly affects the terrestrial planets, particularly Mercury, and can cause the system to destabilize over long periods of time. Laskar (2008) even shows a 1–2% probability of Mercury’s eccentricity being excited to a degree which would risk planetary collision in the next 5 Gyr. What we still don’t fully understand is whether these chaotic symptoms (highly excited eccentricities, close encounters and ejection) are an expected outcome of the planetary formation process as we presently understand it, or merely a quality of our particular solar system. The work of Laskar (2000) showed us that the outcomes of semi-analytic planetary formation models of our own solar system show symptoms of chaos, and are connected to the particular initial mass distribution which is chosen. However these systems were formed without the presence of the gas giants, and planetesimal interactions were simplified to minimize computing time. Perhaps our solar system is a rare outlier in the universe, with it’s nearly stable, yet inherently chaotic system of orbits occurring by pure chance. Of even greater interest, if it turns out that systems like our own are unlikely results of planetary formation, we may need to consider other mechanisms that can drive the terrestrial planets into their modern chaotic state.
This work takes 145 systems of terrestrial planets formed in Kaib & Cowan (2015) as a starting point. The systems are broken into three ensembles. The first set of 50 simulations, “Circular Jupiter and Saturn” (cjs), are formed with Jupiter and Saturn on nearly circular (e<0.01) orbits, at their current semi-major axes. The simulations use 100 self-interacting embryos on nearly circular and coplanar orbits between 0.5 and 4.0 AU, and 1000 smaller non-self-interacting planetesimals. The smaller planetesimals interact with the larger bodies, but not with each other. Additionally, the initial embryo spacing is uniform and embryo mass decreases with semi-major axis to yield an r-3/2 surface density profile. The second ensemble (containing 46 integrations), “Extra Eccentric Jupiter and Saturn” (eejs) evolve from the same initial embryo configuration as cjs, with Jupiter and Saturn initially on higher (e=0.1) eccentricity orbits. The final batch of integrations (49 systems), “Annulus” (ann), begin with Jupiter and Saturn in the same configuration as cjs, however no planetesimals are used. 400 Planetary embryos for ann are confined to a thin annulus between 0.7–1.0 AU, roughly representative of the conditions described following Jupiter’s outward migration in the Grand Tack Model (Walsh et al., 2011b).
After advancing each system to t=4.5 Gyr, we perform detailed 100 Myr simulations and probe multiple chaos indicators. By careful analysis we aim to show whether chaotic systems naturally emerge from accretion models, and whether the source of the chaos is the same as has been shown for our own solar system.
2. Methods
2.1. System Formation and Evolution
We use the simulations modeling terrestrial planet formation in Kaib & Cowan (2015) as a starting point for our current numerical work. In Kaib & Cowan (2015), all simulations are stopped after 200 Myrs of evolution, an integration time similar to previous studies of terrestrial planet formation (e.g. Chambers, 2001a; Raymond et al., 2004; O’Brien et al., 2006; Walsh et al., 2011b). Because we ultimately want to compare the dynamical state of our solar system (a 4.5 Gyr old planetary system) with the dynamical states of our simulated systems, we begin by integrating the systems from Kaib & Cowan (2015) from Myr to Gyr. Since bodies can evolve onto crossing orbits and collide before Gyr, accurately handling close encounters between massive objects is essential. Thus, we use the MERCURY hybrid integrator (Chambers, 1999) to integrate our systems up to Gyr. During these integrations, we use a 6-day timestep and remove bodies if their heliocentric distance exceeds 100 AU. Because we are unable to accurately integrate through very low pericenter passages, objects are also merged with the central star if their heliocentric distance falls below 0.1 AU. Though by no means ideal, the process of removing objects at 0.1 AU is commonplace in direct numerical models of planetary formation due to the limitations of the integrators used for such modeling. Chambers (2001a) showed that this does not affect the ability to accurately form planets in the vicinity of the actual inner solar system, since objects crossing 0.1 AU must have very high eccentricities. These excited objects interact weakly when encountering forming embryos due to their high relative velocity, and rarely contribute to embryo accretion. It should be noted that many discovered exoplanetary systems have planets with semi-major axis interior to 0.1 AU. However, we are not interested in studying such systems since we aim to draw parallels to our actual solar system. The WHFAST integrator used in the second phase of this work (Section 2.2), however, can integrate the innermost planet to arbitrarily high eccentricities, so the 0.1 AU filter is no longer used. Finally, to assess the dynamical chaos among planetary-mass bodies, any “planetesimal” particles (low-mass particles that do not gravitationally interact with each other) that still survive after 4.5 Gyrs are manually removed from the final system.
2.2. Numerical Analysis
Numerical simulations for detailed analysis of the fully evolved systems are performed using the WHFAST integrator in the Python module Rebound (Rein & Liu, 2012). WHFAST (Rein & Tamayo, 2015) is a freely available, next generation Wisdom Holman symplectic integrator (Wisdom, 1981; Wisdom & Holman, 1991; Kinoshita et al., 1991) ideal for this project due to its reduction on the CPU hours required to accurately simulate systems of planets over long timescales. WHFAST’s reduction in error arising from Jacobi coordinate transformations, incorporation of the MEGNO (Mean Exponential Growth factor of Nearby Orbits) parameter, improved energy conservation error and tunable symplectic corrector up to order 11 motivate the integrator choice. The accuracy of many mixed variable symplectic integration routines are degraded by integrating orbits through phases of high eccentricity and low pericenter. For this reason, in Figure 1 we plot the variation in energy from our simulation with the lowest pericenter ( AU). In the upper panel, we plot the energy of the innermost planet, since this should be fixed in the secular regime. We see that energy variations stay well below one part in 103. In the lower panel, we plot the fractional change in the total energy of this system. Again, we find that energy variations rarely exceed one part in 104. Finally, in Figure 2 we show a histogram of the fractional energy change between the start and end of the integration for all of our simulations. For the vast majority of systems, the fractional energy change is far less than one part in 104, and no simulations exceed 103. Values in excess of 104 are from simulations where the eccentricities of the giant planets were artificially inflated and the performance of the integrator is degraded by close encounters. These systems with frequent close encounters are obviously chaotic, so our chaos determination is not affected.
MEGNO is the primary tool for identifying chaotic systems. Introduced in Cincotta et al. (2003), MEGNO represents the time averaged ratio of the derivative of the infinitesimal displacement of an arc of orbit in N-dimensional phase space to the infinitesimal displacement. For quasi-periodic (stable) motion, MEGNO will converge to a value of 2 in the infinite limit. For chaotic systems, however, MEGNO will diverge (Cincotta et al., 2003). Maffione et al. (2013) showed that MEGNO is an extremely useful and accurate tool for detecting chaos. Systems which are non-chaotic will maintain stable MEGNO values of 2 for the duration of the simulation, while chaotic systems diverge from 2. For this project, systems which attained a maximum value of MEGNO 3.0 were classified as chaotic.
For use in certain analyses, the Lyapunov Timescale () is also output. WHFAST calculates the inverse of by least squares fitting the time evolution of MEGNO (Rein & Tamayo, 2015). Systems classified as chaotic tended to have a less than 10–100 Myr.
Some simulations which quickly displayed chaos were terminated early to save computing time. Terminating these simulations early did not affect the chaos determination since MEGNO had already clearly diverged, nor did shorter simulations affect follow on data analysis and reduction (such as the detection of resonances described in section 2.4).
Another tool we use to characterize the chaos in our systems is Angular Momentum Deficit (AMD) (Laskar, 1997). AMD (equation 1) measures the difference between the z-component of the angular momentum of a given system to that of a zero eccentricity, zero inclination system with the same masses and semi-major axes. Evaluating the evolution of AMD over the duration of a simulation will probe whether angular momentum is being exchanged between giant planets and terrestrial planets as orbits excite and deexcite due to induced chaos (Chambers, 2001b).
[TABLE]
2.3. Simulation Parameters
Simulations are run for 100 Myr, with an integration timestep of 3.65 days. Orbital data is output every 5000 years. The order of the symplectic corrector is the order to which the symplectic correction term in the interaction Hamiltonian (dt in Rein & Tamayo (2015)) is expanded. Here, we set this to order 1 (WHFAST allows for corrections up to order 11) (Rein & Tamayo, 2015). 10 sample systems were integrated at different corrector values in order to determine the lowest corrector order necessary to accurately detect chaos. Additionally, a total of 16 1 Gyr simulations consisting of both chaotic and non-chaotic systems of 3, 4 and 5 terrestrial planets were performed to evaluate long-term behavior and verify the adequacy of 100 Myr runs. Finally, 6 sets of 145 simulations are performed, results and findings for which are reported in section 3. The 6 runs are summarized in Table 1.
2.4. Detecting Mean Motion Resonances
A Mean Motion Resonance (MMR) occurs when the periods of orbital revolution of 2 bodies are in integer ratio to one another. For a given MMR, the resonant angle will librate between 2 values (Àngel Jorba, 2012). Many possible resonant angles exist. For this paper, however, we only consider 4 of the more common planar resonant angles (Elliot et al., 2005). Only planar resonances are considered because our systems typically have very low inclinations.
To detect MMRs, the average keplerian period is calculated for all bodies in the simulation. 4 resonant angles are calculated for all sets of bodies with period ratios within 5% of a given integer ratio. 21 different MMRs (all possible permutations of integer ratios between 2:1 and 8:7) are checked for. Using a Komolgorov-Smirnov test, each resulting time-resonant angle distribution (e.g. figure 3) is compared to a uniform distribution (e.g. figure 3d), yielding a p-value. All distributions with p-values less than 0.01 are evaluated by eye for libration. Figure 3 shows 4 different example distributions and their classification.
3. Results
3.1. System Evolution Beyond 200 Myrs
In Kaib & Cowan (2015), systems of terrestrial planets were generated via simulations of terrestrial planet accretion. These simulations were terminated after 200 Myrs of evolution, as each simulation had evolved into a system dominated by 1–6 terrestrial planet-mass bodies. Terminating accretion simulations after 200–400 Myrs of system evolution is common practice since the great majority of accretion events occur well before these final times are reached. However, it remains unknown how these newly formed systems evolve over the next several Gyrs. Do planetesimals and embryos naturally accrete into indefinitely stable configurations of terrestrial planets? Or are the systems that arise from terrestrial planet accretion often only marginally stable, with major instabilities occurring hundreds of Myrs or Gyrs after formation?
To begin answering this question, we take the systems from Kaib & Cowan (2015) and integrate them for another 4.3 Gyrs with MERCURY. In Figure 4, we show the cumulative distributions of times at which these systems lose their last terrestrial planet mass body ( 0.055 M*⊕*). These planets can be lost via collision with a larger planet, collision with the Sun, or ejection from the system ( AU). We find that there are many systems that undergo substantial dynamical evolution after their first 200 Myrs. As Figure 4 shows, between 20 and 50% of systems lose at least 1 planet after Myrs. This fraction varies with the simulation batch. Systems in the cjs set are the most likely to lose planets at late times. This is likely due to the fact that these systems often form planets well beyond 2 AU (Raymond et al., 2009c), where dynamical timescales are longer and planets require a longer time period to undergo ejections or final collisions compared to those at 1 AU. In eejs simulations, a smaller fraction of systems (25%) lose a planet after their first 200 Myrs of evolution. In these systems, the effects of an eccentric Jupiter and Saturn greatly deplete the mass orbiting beyond 1.5–2 AU (Raymond et al., 2009c), and this absence of more distant material may explain the decrease in late instabilities. Finally, the conditions are even more extreme in the ann simulations, where the initial planetesimal region is truncated at 1 AU. These simulations have the lowest rate of late ( Myrs) instabilities at 18%.
It should also be noted that some systems lose planets at extremely late times. 8 out of 150 systems (5%) lose planets after Gyr. 5 of these systems are from cjs, while eejs and ann yield 1 and 2 systems, respectively. This small, yet non-negligible fraction of systems undergoing late instabilities may help explain the existence of transient hot dust around older main sequence stars (Wyatt et al., 2007). These very late instabilities in our systems occur even though the orbits of Jupiter and Saturn are effectively fixed for the entire integration. The rate of instabilities would likely be significantly higher if the orbits of the gas giants evolved substantially over time (Brasser et al., 2009; Agnor & Lin, 2012; Brasser et al., 2013; Kaib & Chambers, 2016).
In Figure 5, we look at the mass distributions for the last planet lost from each system with an instability after Myrs. In general, we see that the last planets lost from systems with late instabilities have masses below 0.5 M*⊕. This is not surprising, since during an instability event it is typical for the smallest planets to be driven to the highest eccentricities, resulting in their collision or ejection (Rasio & Ford, 1996; Chatterjee et al., 2008; Jurić & Tremaine, 2008; Raymond et al., 2009a). However, not all systems abide by this. In particular, 3 of the 13 eejs systems that undergo late instabilities lose planets with masses well over 1 M⊕*.
This suggests there may be a different instability mechanism in eejs systems. Indeed, when we look at how the last planets are lost from eejs systems, we find that 10 of the 13 systems with late instabilities lose their planets via collision with the Sun. This contrasts strongly with the cjs and ann systems, where there is only one instance of a planet-Sun collision among the 34 systems that have late instabilities. Moreover, there are 4 eejs systems with only 1 planet at Myrs, which go on to have a planet-Sun collision before Gyrs. In these cases, the gas giants are clearly driving instabilities. This is not surprising, since the heightened eccentricities of Jupiter and Saturn will enhance the secular and resonant perturbations they impart on the terrestrial planets. When interactions between a gas giant and the terrestrial planets are the main driver of an instability, the relative masses of the terrestrial planets lose their significance because they are all so small relative to the gas giants. This allows for more massive planets to be lost from these systems.
Figure 6A shows an example of an instability within a cjs system. In this case, a system of 5 terrestrial planets are orbiting at virtually fixed semi-major axes for 3.8 Gyrs when an instability develops between the inner 3 planets. The second and third planets collide and the resulting 4-planet system finishes the simulation with smaller orbital eccentricities than it began with. On the other hand, the evolution of an eejs system is shown in Figure 6B. Here we see the eccentricities of 3 relatively well separated planets driven up around 400–500 Myrs, leading to a collision between the second and third planets. After the collision, the outermost planet’s eccentricity is again quickly excited and eventually approaches 0.8. Shortly after this point, the planet collides with the Sun (after a scattering event with the inner planet). While the detailed dynamics of this system are undoubtedly complex, the behavior is clearly different from that of the cjs system, and is almost certainly a consequence of the enhanced gas giant perturbations produced from their increased eccentricities.
We also study how the properties of systems with late instabilities differ from systems that do not lose planet-mass bodies after Myrs. In Figure 7 we look at the number of planets that each system has. In panels A–C, we see that after 200 Myrs of evolution, cjs systems typically have 4–6 planet-mass bodies (an average of 4.68 planets per system). This is significantly higher than eejs and ann systems, which have an average of 2.45 and 3.00 planets per system, respectively. The differences can largely be attributed to the lack of distant planets in these systems, owing to their initial conditions. Panels A–C also show which systems go on to lose planets at later times. For cjs and ann simulations, these systems tend to have more planets than the overall distribution. In contrast, no such trend is seen among eejs systems. Regardless of planet number, the eejs systems all seem to have roughly the same probability of losing a planet at late times. This is again a symptom of the gas giants driving instabilities within these systems, unlike the cjs and ann systems, where interactions between terrestrial planets play a larger role in late instabilities. Finally, panels D–F show the distributions of planets per system after 4.5 Gyrs of evolution. At the end of our integrations, the cjs, eejs, and ann systems have an average of 3.76, 2.12, and 2.78 planets per system respectively. For all of our simulation batches, we see that systems with late instabilities tend to have lower numbers of planets than the overall distribution of systems. Thus, in the case of cjs and ann systems, late instabilities tend to transform systems with relatively high numbers of planets into systems with relatively few planets. We also note that 4 eejs systems finish with no terrestrial planets whatsoever.
Finally, we show the AMD of each of our terrestrial planet systems at Myrs and Gyrs in Figure 8. Panels A–C show our systems’ AMD distribution at Myrs. For each of our simulation batches, the median AMD is greater than the solar system’s value. Our cjs, eejs and ann simulations have median AMD values of 2.9, 4.0, and 1.7 times the value of the modern inner solar system. Again, we also show the AMD distributions for systems that go on to have late instabilities. These systems tend to have larger AMD values. For the cjs, eejs and ann systems that undergo late instabilities the median AMD values at Myrs are 4.1, 14, and 4.3 times the solar system’s AMD, respectively. Interestingly, though, panels D–F demonstrate that these systems are not always destined to maintain a relatively large AMD. Systems in the cjs and ann batches that undergo late instabilities have median AMD values of 2.7 and 3.7 times the value of the solar system after 4.5 Gyrs of evolution, respectively. Thus, a late instability does not necessarily increase the AMD of the system, and in some situations can result in moderate decreases. On the other hand, in eejs systems, the excited orbits of Jupiter and Saturn continue to wreak havoc on the terrestrial planets. Systems that experience late instabilities have a median AMD of 132 times that of the solar system!
3.2. Prevalence of Chaos
A selection of results from our simulations are provided in Appendix A (Table LABEL:table:data). and MEGNO are listed for run 1 for all systems. Additionally, we provide our chaos determination (yes or no) for runs 1, a and b, as well as MMRs detected for run 1. Figure 9 compares the fraction of all systems which are chaotic between runs 1, a and b. We find that removing Jupiter and Saturn has the greatest effect on reducing chaos in our systems. In general, % of systems exhibit some form of chaos, when Saturn is removed only % of systems are chaotic and when Jupiter is removed that number is only %. This indicates that the chaos in most of our systems is likely driven by perturbations from Jupiter. In fact, when Jupiter was removed, all systems but 1 had ’s which either increase, or are within 1.5 orders of magnitude of the original value.
Figure 9 also clearly shows the disparity of chaos between systems with different numbers of terrestrial planets. This is most pronounced in 5-planet configurations, where only 2 such systems are free of chaos with the outer planets in place, and only 3 when they are removed. Having 4 terrestrial planets may not to be a significant source of chaos in our own solar system. This can be seen in figure 10, which provides MEGNO plots for the solar system with and without the 4 outer planets. In fact, the solar system may better be described as a 3-planet configuration when compared to our results. Most 4 and 5-planet systems in our study differ greatly from our own since they typically contain only planets with masses comparable to Earth and Venus (see Table 3 for 2 such 5-planet examples). Mars analogues are rare in our systems, and Mercury sized planets are almost non-existent. If we consider our solar system a 3 terrestrial planet arrangement, it’s inherent chaos fits in well with our results; where about half of systems show chaos with the giant planets in place, and only around 1 in 6 when they are removed. A shortcoming of this comparison is that many studies have shown that Mercury is a very important source of the chaos in our own solar system (Laskar, 2008; Laskar & Gastineau, 2009). However, when we integrate the solar system without Mercury, the system is still chaotic. Therefore, though the actual solar system does match our results, this comparison is limited by the fact that present models of the terrestrial planet formation systematically fail to produce Mercury analogs.
MMRs between planets are common features in many of our chaotic systems, implying that they are often important sources of the dynamical chaos. We detect 365 MMRs among all simulations in this phase of the project, 82% of which occur in chaotic systems. Further analysis shows that the MMRs which do occur in non-chaotic systems tend to be of higher order between smaller terrestrial planets. It should be noted that the vast majority of these MMRs are intermittent, and last only a fraction of the entire simulation duration.
Figure 11 shows the fraction of systems which are chaotic in runs 1, c, d and e. It is clear that an eccentric Jupiter and Saturn can quickly introduce chaos to an otherwise non-chaotic system. One interesting result from this batch of simulations is that when the eccentricities of the outer planets are inflated, the likelihood of a 5:2 MMR between Jupiter and Saturn developing increases. In almost all systems, this resonant perturbation introduces chaos, and can possibly destabilize the system. Since systems labeled cjs were formed with the giant planets on near circular orbits, simply multiplying the already low eccentricity by 1.5 or 2 was not enough to produce a noticeable effect. For this reason, run e was performed using a step increase of 0.05. There is a clear parallel between the results of this scenario and a Nice Model instability (Gomes et al., 2005; Morbidelli et al., 2005; Tsiganis et al., 2005), where the outer planets rapidly transition from nearly circular to relatively eccentric orbits.
To further probe this effect, we repeat run e for cjs and ann systems using the MERCURY hybrid integrator in order to accurately detect collisions and ejections. Systems are integrated for 1 Gyr using simulation parameters similar to those discussed in section 2.1. We find instabilities are relatively common in these systems. 29% of ann systems and 52% of cjs systems lose one or more planets over the 1 Gyr integration. Table 2 shows the percentage of systems which lose a given number of terrestrial planets. In fact, the resulting systems are quite similar to those produced after integrating the eejs batch to the current epoch. A small fraction of systems lose all inner planets, and some can have instabilities occur very late in the simulations (Figure 12). Overall we show that an event similar to the Nice Model scenario, where the Giant planets eccentricities quickly inflate, can result in a non-negligible probability of inner planet loss.
If we again classify our solar system as a 3 terrestrial planet system, we can draw further parallels between the results of such cjs run e configurations and the Nice Model instability, since Jupiter and Saturn begin on near circular orbits. In run e, 6 out of 8 such systems were chaotic, as compared to only 2 out of 8 in run 1. Indeed, we see that even a perturbation in the giant planet’s eccentricities of 0.05 is successful in rapidly making a system chaotic. In fact, when our own solar system is integrated with the outer planets on circular (e<0.001, i0) orbits, the chaos disappears.
3.3. Angular Momentum Deficit and Last Loss Analysis
For our systems, we evaluate the difference between the average AMD of the inner planets over the first and last 3 Myr of our 100 Myr simulations. Taking the average removes the contributions from periodic forcing in the AMD from Jupiter and Saturn. We find a weak trend for chaotic 4 and 5 planet systems to have larger changes in AMD over the duration of the simulation than their non-chaotic counterparts. Of the 13 systems which had total changes in AMD greater than the actual solar system’s value, 9 were classified as chaotic. The largest outlier, a non-chaotic eejs 2 planet system (eejs25), is discussed further in section 3.5.1. We also search for any correlation between the time and mass of the last object ( 0.055 M*⊕*) lost, and the chaos of a system. Though we show that some unstable, chaotic systems can stabilize after losing a planetary mass body, we are unable to identify any conclusive trends as to whether a late instability will shape the ultimate chaotic state of a system.
3.4. Mass Concentration Statistic and Center of Mass Analysis
To evaluate the degree to which mass is concentrated at a given distance away from the central star, we utilize a mass concentration statistic () (Chambers, 2001b):
[TABLE]
The expression in parenthesis in (2) is essentially the level of mass concentration at any point as a function of semi-major axis. Chambers (2001b) utilizes a logarithm in the equation since, in our own solar system, the semi-major axes of the planets lie spaced in rough geometric series (the famous Titus-Bode law). is the maximum value of the mass concentration function. A system where most of the mass is concentrated in a single, massive planet would have a very steep mass concentration curve, and a high value of . A system of multiple planets with the same mass would have a smoother curve and yield a lower . As a point of reference, the value of the solar system’s 4 inner planets, where most of the mass is concentrated in Venus and Earth, is 90. values are provided in the same format as AMD values in figure 13. A general, weak correlation can be seen between chaotic systems and slightly higher values of , however this trend is not very conclusive. We also provide the center of mass for each system of terrestrial planets in figure 13. A clear trend is visible where non-chaotic systems tend to have a center of mass between 0.8–1.2 AU (the value being slightly greater as number of planets increases). In general, the more mass is concentrated closer to the central star, or closer to Jupiter, the greater the likelihood of chaos developing. This is likely related to Jupiter’s role in introducing chaos to systems. This trend is true for all three simulation subsets, but particularly strong in the cjs and eejs batches.
3.5. Systems of Particular Interest
3.5.1 eejs25
The system with the largest change in AMD over the duration of the simulation is surprisingly non-chaotic. This outlier (eejs25), is a system of just 2 inner planets. The innermost planet is 117% the mass of Earth, residing at a semi-major axis of 0.62 AU, and the second planet is 96% the mass of Venus at a semi-major axis of 1.35 AU. The innermost planet is locked in a strong, secularly driven resonance with Jupiter (figure 14). This causes the eccentricity of the innermost planet to periodically oscillate between 0.15 and 0.7 over a period of 8 Myr. These oscillations are remarkably stable. In fact, due to the fortuitous spacing between of the Sun and inner 2 planets, this oscillation does not lead to interactions with other bodies in the system.
3.5.2 cjs10 and cjs13
The 2 most stable 5 planet configurations occurred in cjs10 and cjs13, with both systems classified as non-chaotic through runs 1, a and b. In the runs which vary eccentricity, cjs10 started to develop a weak 5:2 MMR between Jupiter and Saturn, causing mild chaos in run c (where eccentricities are increased by 150%). However, the chaos in this run was mild (MEGNO only rose to 3.703 and for this run was ). A step increase of 0.05 to Jupiter and Saturn’s eccentricity in run e was required to fully introduce chaos (maximum MEGNO values 198.9 and 186.4 for cjs10 and cjs13 respectively) to both of these systems. This excitation of the eccentricities of the giant planets drove an occasional 3:1 MMR between the second and fourth inner planets in cjs13, possibly contributing to this chaos. The most remarkable similarity between these systems is their mass spacing and distribution (summarized in Table 3). Both have similarly low values of (20.1 and 15.6). In fact, the planet spacing of both systems is somewhat reminiscent of a Titus-Bode Law series. For example, all orbital locations of cjs13 are within 6% of a Titus-Bode series beginning at the inner planet’s semi-major axis.
4. Discussion and Conclusions
We have presented an analysis of systems of terrestrial planets formed through direct numerical integration of terrestrial accretion, fully evolved to the present epoch. Our work aims to assess whether our solar system, and it’s inherently chaotic dynamics, is a likely result of planetary formation as we currently understand it. We report that roughly half of our systems display some form of chaos. By far, the most common source of this dynamical chaos is perturbations from Jupiter. Additionally, we find that systems in our sample with greater numbers of terrestrial planets are far more prone to chaos than those with fewer inner planets. Unfortunately, systems formed through numerical integrations (including those of Kaib & Cowan (2015) which are used for this work) still routinely produce Mercury and Mars analogues which are far too massive (Chambers, 2001a; O’Brien et al., 2006; Chambers, 2007; Raymond et al., 2009b). Consequently, we find it best to consider our solar system a 3 terrestrial planet system for the purposes of comparison in this work. This classification of the solar system works well with our results that 3-planet systems have an 50% chance of being chaotic, and a much lower probability when the giant planets are removed.
By varying the eccentricities of Jupiter and Saturn in 3 separate batches of simulations, we show that an eccentric system of outer planets can quickly introduce dynamical chaos and trigger instabilities in otherwise stable systems. This result confirms the findings in numerous previous works (e.g. Gomes et al., 2005; Morbidelli et al., 2005; Tsiganis et al., 2005). The inflation in eccentricity required to create such a chaotic system is surprisingly small. By varying the eccentricity of a batch of systems with Jupiter and Saturn on nearly circular orbits, we show that dynamical chaos quickly ensues. This sort of event is akin to a Nice Model-like instability (Gomes et al., 2005; Morbidelli et al., 2005; Tsiganis et al., 2005). We go on to show that in such an instability, the possibility of destabilizing the inner planets to the point where a terrestrial planet is lost by either collision or ejection is fairly high.
Additionally, we find that systems most immune to developing dynamical chaos tend to have centers of mass between 0.8–1.2 AU, though that range is by no means absolute. This is an interesting result, and likely related to Jupiter’s role in driving chaos in many of these systems.
We consistently identified systems throughout our suite of simulations which displayed many of the same chaotic dynamics as our own solar system. It is clear that chaotic systems such as our own are common results of planetary formation. The largest source of chaos in our own system, perturbations from Jupiter, is the most common source of chaos observed in our work. The solar system, however, is akin to only a small fraction (10%) of our simulations since removing just the planets beyond Jupiter turns our system non-chaotic. Additionally we show that late instabilities are common among these systems, and it is not far-fetched to imagine a late instability shaping dynamics within our own system. Finally, we find many systems with similar numbers of terrestrial planets, semi-major axis configurations, mass concentrations and chaos indicators ( and MEGNO) as our own.
Acknowledgements
This work was supported by NSF award AST-1615975. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. The bulk of our simulations were performed over a network managed with the HTCondor software package (https://research.cs.wisc.edu/htcondor/).
Appendix A Supplemental Data
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Agnor & Lin (2012) Agnor, C. B., & Lin, D. N. C. 2012, Ap J, 745, 143
- 2Àngel Jorba (2012) Àngel Jorba. 2012, UNESCO Encyclopedia of Life Support Systems, Vol. 6.119.55 Celestial Mechanics
- 3Brasser et al. (2009) Brasser, R., Morbidelli, A., Gomes, R., Tsiganis, K., & Levison, H. F. 2009, A&A, 507, 1053
- 4Brasser et al. (2013) Brasser, R., Walsh, K. J., & Nesvorný, D. 2013, MNRAS, 433, 3417
- 5Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
- 6Chambers (2001 a) —. 2001 a, Icarus, 152, 205
- 7Chambers (2001 b) —. 2001 b, Icarus, 152, 205
- 8Chambers (2007) —. 2007, Icarus, 189, 386