N-Body Simulations of Chaos in the Orbits of Trappist System Planets

[charlesM]

TL;DR Summary

Using the "Solar system" C++ code of the Boost-ODEINT framework as a starting point I have simulated the orbits of the seven Trappist system planets. Chaos effects are seen on the scale of tens of years compatible with being observed by the Webb Telescope. The planets remain in stable orbits for millions of years. Comparable n-body simulations were performed for the Solar system planets reproducing known n-body effects there. Advice from professional astronomers would be appreciated.

As a retired physics professor with a long experience in complex simulation software for high energy physics experiments (e.g. the LHC) I revisited last July the n-body planetary simulations which I taught in an undergraduate physics course during the Spring 2017.

It was then that the discovery of the seven-planet Trappist system was announced. My renewed interest was prompted by news of the success of the Webb Telescope which will devote a significant effort to observing the Trappist dwarf star system of seven resonance-coupled planets. This project to simulate the Trappist system took several months of study. The main result of the study was that chaos plays an important role in the orbits of those planets, just as it does in the orbits of the inner planets of the Solar system, especially Mercury.

Unlike the Solar system chaos effects in the Trappist system should be observable over much shorter time spans, decades in fact. In order to do quality control checks on these results, the software used for the Trappist system planets was extended to study the Solar system n-body effects such as the precession of the perihelion of Mercury's orbit, the expulsion of Kirkwood gap asteroids, and the evolution of the eccentricities of the planets over ten million years.

The descriptions and results of the work can be found at the link https://www.dropbox.com/s/vu7x0oc8k2rbqj6/trappistSystemSimulations_11January2023.pdf?dl=0, with the last four pages 168-171 containing the summary.

As a non-professional in astronomy and well aware that flawed software can produce seemingly novel but false results I would appreciate comments from professionals in the field. Thanks.

 

[Drakkith]

Very nice! I wish I could give you some professional advice, but I'm certainly not a professional astronomer.

 

[Filip Larsen]

As an Engineer, with a bit of academic background and interest in simulation of chaos in mostly dissipative systems (i.e. not planetary motion), I am curios if you have analyzed if there are chaos present not just in the position of a planet for some initial configuration and state, but also in one or more of the shape parameters of the osculating orbit (e.g. the semi-major distance)?

Since chaos and periodic phase space orbits go hand in hand I am also curios if you have spotted periodic (resonant) solutions in your simulations? While periodic solutions are unlikely to emerge on "random" initial conditions, there are numerical tools that can assist in finding them.

And finally I am curious if you have analyzed or at least have an indication on how valid your results are when performed in a 2D space as compared to 3D? E.g. if you have chaos for 2D in a fairly large region of configuration space, can you then extend that region to 3D and show same chaotic measure?

By the way, for the numerical study of dissipative systems there exists a set of mathematical/numerical tools that really can speed up obtaining results, e.g. tracing out the regions of configuration space that exhibit similar behavior (e.g. stable/unstable manifold tracing), finding periodic solutions, etc. I assume analysis of conservative systems have a similar and somewhat overlapping set of tools that might benefit you.

 

[charlesM]

Very nice! I wish I could give you some professional advice, but I'm certainly not a professional astronomer.

Thanks for telling me that the post was interesting for you.

 

[charlesM]

As an Engineer, with a bit of academic background and interest in simulation of chaos in mostly dissipative systems (i.e. not planetary motion), I am curios if you have analyzed if there are chaos present not just in the position of a planet for some initial configuration and state, but also in one or more of the shape parameters of the osculating orbit (e.g. the semi-major distance)?

Since chaos and periodic phase space orbits go hand in hand I am also curios if you have spotted periodic (resonant) solutions in your simulations? While periodic solutions are unlikely to emerge on "random" initial conditions, there are numerical tools that can assist in finding them.

And finally I am curious if you have analyzed or at least have an indication on how valid your results are when performed in a 2D space as compared to 3D? E.g. if you have chaos for 2D in a fairly large region of configuration space, can you then extend that region to 3D and show same chaotic measure?

By the way, for the numerical study of dissipative systems there exists a set of mathematical/numerical tools that really can speed up obtaining results, e.g. tracing out the regions of configuration space that exhibit similar behavior (e.g. stable/unstable manifold tracing), finding periodic solutions, etc. I assume analysis of conservative systems have a similar and somewhat overlapping set of tools that might benefit you.

So far for the Trappist system planets I have looked at the chaos onsets only for the two orbital coordinates: radius and angle. However, it will be straightforward to check the semi-major axis time behavior when there are small changes made in the initial conditions. I'll post those results when they are done. For the Solar system planets I did look at the eccentricity shape parameter changing with time according to initial conditions, quite large for Mercury and much smaller for the other planets. I could do that too for the Kirkwood gap asteroid simulations, for which I did look at the orbital coordinate changes which are dramatic.

For the pre-Webb Trappist planet published orbit parameters, the semi-major axis and orbital times quoted values are quite precise: 0.86% for the semi-major axis numbers and 0.001% or less (one sigma) for the orbit times. My expectation, if there is significant chaos, that after several years the planets will not be where they are predicted to be for the occultation observations.

As for the extension from 2D to 3D that would be a numerical exercise for the Trappist planet system because I do not believe they can measure the respective inclination angles? For the Solar system that extension is on my to-do list since the ecliptic plane inclination angles are known.

Yes, it would be good for me to learn more about the tools available for chaotic analysis. Like the astronomy subject my knowledge of chaos derives mainly from teaching the undergraduate mechanics and computational physics courses. This included lectures on the damped and driven pendulum going in and out of chaotic motion according to the strength of the driving force, Poincare' plots, and Lyapunov exponents. But that is really only a superficial understanding.

 

[Filip Larsen]

I have looked at the chaos onsets only for the two orbital coordinates: radius and angle.

So, just so I understand, you classify the phase space trajectory pair for a particular initial condition (position and velocity for all planets) for a particular parameter set (number of planets and their mass ratios) as chaotic if the central radius or angle measure for the state of one of the planets decorrelate after some time, and non-chaotic if it does not (within the time limit)? The reason I ask is when reading through your result it was unclear to me if the Luapunov exponent of the trajectory flow was calculated or not.

Also, it seems that you use the term "onset of chaos" to mean that the trajectory pair at some point in time becomes uncorrelated, and at that point in time there is then "onset of chaos". When flipping through your result I was confused about this since "onset of chaos" for me means the boundary in parameter space where the system goes from having no chaotic trajectories to having them. For example, for the logistic map you mention, the onset of chaos would then be the parameter values from where chaotic trajectories begin to emerge.

While I obviously are not able to tell if a professional astronomer can dig right into your results without further explanation, I think it at least would benefit most general science readers on how to understand your results if there was a bit more precise description of the definitions and methods employed to produce those results. I guess I am a bit like, "ok, nice looking results, but where is the accompanying paper so I can understand what they mean?"

 

[charlesM]

I have competed a set of supplemental simulations for the Trappist System planet orbits, as I had indicated I would do. These simulations are the butterfly effect comparisons over 300 years for the eccentricities and the semi-major axis shape parameters of the elliptical orbits of the seven planets. The butterfly effect initiator was a pico-scale change in the velocity of Planet d, leaving all other initial coordinates the same for the comparison simulation. A 0.007 day time step was used. As a reminder a pico-scale change in the velocity coordinate is orders of magnitude smaller than what is empirically attainable.

The results of these comparison simulations are qualitatively the same as previously shown for the two position coordinates (radius and angle) of the orbits. The compared simulations have imperceptibly small differences until about 40 years at which time substantially large differences become visible. Moreover, all seven planets diverge at about the same time, although as before the inner two planets (b and c) show smaller effects as was the case for the position coordinates.

A new full-length (188 pages) version of this work is at the Dropbox link (also correcting a mistake on slide 12)

https://www.dropbox.com/s/9re03h9c4uj6327/trappistSystemSimulations_16January2023 copy.pdf?dl=0

The smaller (17 pages) set of supplemental slides is at the Dropbox link

https://www.dropbox.com/s/zk337hyuk...SystemSimulations_16January2023 copy.pdf?dl=0

This behavior of the shape parameters is not unexpected. Both parameters can be calculated from the values of the position coordinates and their time derivatives (velocities). So if the coordinates diverge at a given time it is likely that the shape parameters will too. Additionally, I don't believe that it is a surprise that all planet orbits diverge simultaneously. It is a tightly coupled system in stable equilibrium with a nominal 120-parameter potential energy function of the 16 Cartesian position coordinates of the eight masses in the center-of-mass frame, although that number can be reduced slightly because of the center-of-mass motion constraint. The coordinates re-adjust continuously to oscillate about the stable minimum point.

Some physics points for the general reader. In a simple two-body gravitational system the semi-major axis is equal to the product of Newton's Universal Gravitational constant times the central star's mass times the planet's mass all divided by twice the absolute value of the total energy in the center-of-mass system. For a bound orbit the total energy is negative. Hence the variations that are seen in the value of the semi-major axis, even in the first 40 years, represent exchanges of energy among the masses. This exchange of energy phenomenon in n-body gravitational systems was used by NASA in the "slingshot" or gravity assist trajectories of spacecraft such as Voyager 2 or ISEE-3 to gain speed or redirection with a minimum of fuel expenditure. This phenomenon is also the mechanism by which Jupiter can expel Kirkwood resonance gap asteroids which were originally put in bound orbits. Expulsion from the Solar System will also likely be the fate of the recently publicized "planet killer" 2022 AP7 asteroid because of close encounters with Jupiter long before the asteroid collides with an inner planet.

The final two slides of the supplemental presentation show the effect of turning off the inter-planet gravity force while retaining the gravity force of the Dwarf Star. The first of these slides shows that the divergence effect for the semi-major axis has been reduced to near imperceptible levels. However, there are still visible systematic variations in the value of the parameter itself, although nearly completely correlated between the two simulations. This indicates that there are still n-body energy exchanges occurring, but not directly. The reason for the persistence of the n-body effect in the absence of direct inter-planet forces is that each planet is exerting a reaction force on the Dwarf Star making it move slightly in the center-of-mass and thus changing the gravity force felt by the other masses. One way to minimize this Newton's Third Law effect is to reduce substantially all the planet masses. That calculation is shown in the last supplemental slide. The variation in the value of the semi-major axis disappears.

At this point I don't plan on making soon any long posts, although I will answer specific questions. I do have to study how chaos has been quantified in already published Solar System papers, apart from displaying the obvious large variations of Mercury's eccentricity over a long time span. I have also become acquainted within the past few days of another n-body simulation package for the Solar System which is said to have a much better symplectic integrator error correction. Duplicating the simulations with an independent code is a prudent action to address concerns about imperfect software.

Lastly, on the subject of imperfect software and on a subject nearer to my own field, I point out the recent article in Nature magazine (https://www.nature.com/articles/d41586-022-04545-z) on the disappearance of a tantalizing hint of a flaw in the Standard Model of High Energy physics. This supposed anomaly was a signature result of the LHCb experiment collaboration at the Large Hadron Collider accelerator in the CERN laboratory near Geneva. Much theoretical work had been done to understand the anomaly as a window into new physics and its vanishing is a disappointment. The details of how the LHCb physicists themselves, as part of the continual self-checking of any group, were able to make the effect go away are given in the Nature paper. Essentially it comes down to imperfect software for identifying pairs of electrons and positrons. The software was not wrong but it was just not as robust as had been thought. Too many pairs which were not electron-positron were being accepted, and this contaminant background fooled the experimenters initially.

 

[charlesM]

This is a major update to the present work simulating chaos effects in the orbits of the seven Trappist System planets. As I mentioned in my previous post I had just become acquainted with a newer n-body simulation code with better integration algorithms. That code is documented at the web site http://rebound.readthedocs.io/en/latest . The code has been developed by university researchers in the field and has already been used in two published works studying the Trappist System, although not the specific topics here. Given the good documentation and examples, it was a fast learning curve to enable me to duplicate my months-long studies with the modified Boost-ODEINT code previously used. The first conclusion is that the major results stated before about the effects of chaos in Trappist System planets have been confirmed with the REBOUND code. I need to point that any mistakes that I have made in using the REBOUND code or interpreting its results are solely my responsibility.

This REBOUND simulation work is contained in a 25-page presentation file located at the Dropbox link

https://www.dropbox.com/s/yszlju4spuh68rs/reboundTrappistSimulations_30January2023.pdf?dl=0

The revised complete presentation with 141 pages is at the link

https://www.dropbox.com/s/8pwq3x59jd3uwkf/trappistSystemSimulations_30January2023.pdf?dl=0

The revised presentation drops all the previous Solar System simulations. Those were included at the time as a "quality control" measure to verify that the Boost-ODEINT program would reproduce known Solar System n-body phenomena. That quality control function has been replaced by the use of the state-of-the-art REBOUND code. The REBOUND code is very well suited to continue the Solar System simulations.

There are some new additions to the simulations. For the "butterfly effect" of making small changes in the velocity of one planet to induce chaotic divergence at a later time, I have incorporated the use of a Gaussian (normal) distribution of precision errors in all the planet velocities to mirror what is the real physical situation. One doesn't know only one velocity imperfectly but rather one doesn't know all the velocities precisely. The end effect is the same. Depending on how much imprecision there is in knowing the initial velocities (and also the positions for that matter) then the predicted orbital positions and parameters (semi-major axis and eccentricity) will be observed to diverge drastically in a few years. The divergence will be earlier the less precisely known are the initial conditions.

A second addition is in the artificial changes in the mass of one or the other planets. These simulations lead to instability of the system within less than one million years. These instabilities were confirmed by the REBOUND software with the use of two more integration algorithms. The stability lifetimes varied in all cases, as one would expect since the different integration algorithms will be producing tiny numerical differences, an intrinsic butterfly effect of their own making. Looking at these results I realized that by reducing the size of the mass changes, one could check to see if the planetary system lived much longer. What sizes changes would be compatible with 10 MY lifetimes or 100 MY lifetimes? Such simulation studies are feasible on the scale of a few or many weeks. For the actual Trappist System seven planets, their masses vary over only a relative narrow range: from 0.4 to 1.6 times their average mass. Possibly the narrow range is a necessary condition to the system having GY lifetimes, along with their resonance coupling.

A final point is on quantifying the temporal character of the chaos. The REBOUND documentation site does have a subsection on "Chaos Detectors". I will be checking that section.

 

[Hornbein]

A miniature solar system! I had no idea. I can see why this interests you.

A question that might be easy for you to answer is, is it possible to have a solar system with two planes of orbit, presumably perpendicular?

Recently it was found that the binary system X-1 has something somewhat like that. X-1 is a big star orbiting a black hole. The orbital plane is perpendicular to the spin of the black hole. Not at all the same thing, but suggests that a protoplanetary disc can have two perpendicular planes of rotation.

 

[charlesM]

As to your first point about the interest in the Trappist System, there are many astronomers, astrophysicists, and astrobiologists interested in that system and who are awaiting stream of results from the Webb telescope for which this system will special observing time. As you may know, perhaps one or two of the planets may be in the habitable zone allowing liquid water. In turn, if the atmosphere and other conditions are right, then life could have developed on that planet.

As to the second point about perpendicular planes of orbit in a planetary system, I'll toss that question over to astronomers who specialize in planetary system evolution. The main physics constraint is that angular momentum must be conserved. So getting some of the planet orbits to have their angular momentum align with that of the central star, while one or more other planet orbits align in a perpendicular direction seems to me a stretch. But I don't have a research background in that field. However, it is straightforward to run a simulation where an orbit is placed initially in a plane perpendicular to the others for a non-binary star system and see how that evolves. In fact one of the constraints in believing a simulation output is to check that the total system angular momentum, like the total system energy, remains conserved to a high accuracy since only internal forces are acting. Whatever happens to the orbital planes, the total angular momentum about the central star cannot change.

 

[Hornbein]

As to your first point about the interest in the Trappist System, there are many astronomers, astrophysicists, and astrobiologists interested in that system and who are awaiting stream of results from the Webb telescope for which this system will special observing time. As you may know, perhaps one or two of the planets may be in the habitable zone allowing liquid water. In turn, if the atmosphere and other conditions are right, then life could have developed on that planet.

As to the second point about perpendicular planes of orbit in a planetary system, I'll toss that question over to astronomers who specialize in planetary system evolution. The main physics constraint is that angular momentum must be conserved. So getting some of the planet orbits to have their angular momentum align with that of the central star, while one or more other planet orbits align in a perpendicular direction seems to me a stretch. But I don't have a research background in that field. However, it is straightforward to run a simulation where an orbit is placed initially in a plane perpendicular to the others for a non-binary star system and see how that evolves. In fact one of the constraints in believing a simulation output is to check that the total system angular momentum, like the total system energy, remains conserved to a high accuracy since only internal forces are acting. Whatever happens to the orbital planes, the total angular momentum about the central star cannot change.

A barrier is that the evolution of binary systems is one of the great unsolved mysteries of astronomy. For example, the heavenly bodies in X-1 are closer together than Mercury is to our Sun. In the past they were much larger. So much so they would have intersected with one another. How can this be?
It seems to me that if there is a binary system so close like that then the two suns with perpendicular rotation would gravitationally be much like one star to distant planets. Each sun is a planar solar system and the planes are perpendicular. The idea is that the protoplanetary discs were rarified enough that they could have these intersecting rotations.

At any rate, I'm hoping you will be curious and try out a few systems like that. Maybe it will appear hopeless, maybe it will be surprisingly stable. I've developed an allergy to learning new software packages -- I've done this one or two too many times in my life -- so I'm hoping it will be so easy for you that you will try it as a lark.

 

[charlesM]

In my previous post I used REBOUND integrator simulations to compare with the prior Boost-ODEINT integrator simulations. These comparisons confirmed that chaotic effects will appear in the Trappist Dwarf Start seven-planet system on a time scale visible within the expected lifetime of Webb Telescope observations. On this topic I have a new, follow-up presentation “Chaos Manifestations in n-Body Integrators used by REBOUND to Simulate Planet Orbits in the Trappist Dwarf Star System”. This new 47-page presentation is largely self-contained, quoting previous results as necessary. It is available at the Dropbox link

https://www.dropbox.com/s/uc5fk1rkyryeu2x/integratorChaosREBOUND _18February2023.pdf?dl=0In the prior work it was seen that “butterfly effect” changes in the initial velocities of the planets would produce abrupt, chaotically unpredictable changes in the orbital coordinates, semi-major axes, and eccentricities of all seven planets 10 or 40 years later simultaneously, for velocity changes at the one-millionth or one-trillionth size, respectively. Present (2020) experimental errors for the orbital time of Planet d, for example, are at the six parts in one million level, while the error in the semi-major axis is at the 8 parts in one thousand level.

The REBOUND software system, for which references are given in the presentation, is useful for a comparative study of different integrator methods used to solve the coupled differential equations describing an n-body planetary orbits motion. Each of these integrators uses somewhat different algorithms to do the numerical solution of these equations. Hence they constitute their own intrinsic butterfly effect when being compared for chaotic systems. In non-chaotic systems, the small numerical differences in the algorithms will produce negligible differences in the outputs, or slowly varying differences depending on the innate accuracy (second order, fourth order, …, fifteenth order) of the integrator for a given time step size. In the presentation it is shown that successive 100-year comparisons of different integrator calculations in the Trappist systems all show abrupt, drastically increasing differences in the semi-major axis and eccentricity predictions simultaneously on time scales of a few years (for the least accurate) , 20 years, 30 years, or 60 years (for the most accurate) integrators. In other words, the different integrators manifest divergent results on about the same time scale as the velocity changes studied before.

Two sets of quality control checks have been done. In the first, the masses of the Trappist planets were all reduced by a factor of one million and the simulation comparisons were repeated. This has the effect of reducing the inter-planet force by a factor of one trillion while reducing the Dwarf Star gravitational force by only one million. Hence the Dwarf Star force becomes relatively one million times stronger than the inter-planet force. Those simulations, done for 300 years, show no divergent results even at highly magnified plotting scales. The second quality control check was to repeat the simulation comparisons for the Solar System planets, this time for a 10 million year span. To put the 100-year Trappist simulations and the 10 MY Solar System calculations into context, in 100 years the Planet d with an orbit time of 4 days will traverse about 9,100 orbits. Obviously the Earth will traverse 10 million orbits in 10 MY. In spite of the huge disparity in orbit numbers, the Solar System simulations showed no divergent behavior in that time range even at a highly magnified scale. Although there is known to be chaos in the Solar System planetary orbits, that fact wasn’t discovered until the 1980s even though astronomers had been doing observations and accurate calculations for centuries (precession of Mercury’s perihelion, predicted location of the then undiscovered planet Neptune). Chaos is not an easily observed effect in the Solar System.For astronomy purposes there was introduced twenty years ago a now well-recognized metric called the MEGNO statistic, standing for Mean Exponential Growth of Nearby Orbits, to quantify the time scale of chaos. Simplistically stated this dimensionless metric, symbolized as <Y>, is calculated over time as 2 asymptotically if the planetary system is non-chaotic. For a chaotic system, looked at over time scales tens of thousands of times longer than that of the longest orbit time, the time dependence of the metric will go as <Y>(t) ~ at + d where d is approximately 0 and the coefficient a is equal to Lambda/2. The parameter Lambda in inverse time units is the Lyapunov Characteristic Number (LCN). The LCN is a well-understood parameter for describing chaotic systems, though previously difficult to calculate for planetary systems before the MEGNO method was introduced. Fortunately, the REBOUND software system does provide a tool to calculate the MEGNO statistic over time. I have recently started to use that tool for the Solar System simulations over ranges greater than 10 MY. I’ll report on those results after I gain confidence in using the tool for both the Solar System and the Trappist Dwarf Star system.

 

[charlesM]

Following up on my previous post mentioning the MEGNO statistic as a measurement of chaos presence in planetary systems I am making available a new 17-page presentation titled "Calculating the MEGNO Metric and the Lyapunov Time Chaos Indicators for the Trappist Dwarf Star Seven Planet System". The presentation can be obtained from the Dropbox link

https://www.dropbox.com/s/ksea4o8q2wbzpui/megnoTrappistSystem_15March2023.pdf?dl=0

The result that I obtain for the Lyapunov Time of the Trappist System is 3.16 +/- 0.30 years, well within the lifetime of the Webb Telescope which will be observing this system. From the literature the definition of the Lyapunov Time is "an estimate of the time after which a system becomes chaotic and two nearby orbits start to diverge by a factor e". The factor e is the transcendental number base of natural logarithms, 2.71828... The Lyapunov Time can be derived from the linear time dependence of the MEGNO metric (Mean Exponential Growth factor of Nearby Orbits) for chaotic planetary systems. The MEGNO number can be repeatedly calculated in the REBOUND n-body simulation software package after a prescribed number of integration time steps have elapsed. Several examples of the MEGNO metric time dependence are given in the presentation, along with the linear fitting values used to extract the related Lyapunov Time.

For comparison the Lyapunov Time is also obtained for the four outer giant planets of the Solar System, in the absence of the inner four planets. That system Lyapunov Time was found to be 1.59 million years by this method, with an uncertainty to be determined. I am studying groups of adjacent Solar System planets.

I did a search on the Cornell arXiv site looking for "Trappist Lyapunov" and "Trappist MEGNO", coming up with one January 15, 2020 entry which was studying something related ("instability time distributions") but not directly the Lyapunov Time extraction nor any mention of the MEGNO metric.