Exploring the orbital evolution of planetary systems

The aim of this paper is to encourage the use of orbital integrators in the classroom to discover and understand the long term dynamical evolution of systems of orbiting bodies. We show how to perform numerical simulations and how to handle output data in order to reveal the dynamical mechanisms that dominate the evolution of arbitrary planetary systems in timescales of millions of years using a simple but efficient numerical integrator. Through some examples we reveal the fundamental properties of planetary systems: the time evolution of the orbital elements, the free and forced modes that drive oscillations in eccentricity and inclination, the fundamental frequencies of the system, the role of the angular momenta, the invariable plane, orbital resonances, and the Kozai–Lidov mechanism.


Introduction
With few exceptions, astronomers cannot conduct experiments and are instead limited to observe the universe. The laboratory for the astronomer usually takes the form of computer simulations. This is the most important instrument for the study of the dynamical behavior of gravitationally interacting bodies. A planetary system, for example, evolves mostly due to gravity acting over very long timescales generating what is called secular evolution. This secular evolution can be deduced analytically by means of the theory of perturbations, but can also be explored in the classroom using precise numerical integrators. Some facilities exist to visualize and experiment with the gravitational interactions between massive bodies [1][2][3]. All of them can show the dramatic dynamical effect due to strong gravitational interactions in a short timescale. However, in general, they are not devised to study the slow changes that planetary systems exhibit in timescales of millions of years or longer. On the other hand, astronomers have several precise orbital integrators to study the long term dynamics of planetary systems [4,5], but their handling require some expertise because they are not thought to be used for educational purposes. Despite the existence of excellent examples halfway between those extremes [6,7], we prefer to present a simpler numerical integrator whose code can be easily examined and modified. The tool we will use here, Orbital Evolution (ORBE), was designed to be precise and fast in the computation of long term interactions, and very simple to use: it contains only one plain text input file with initial conditions and only one plain text output file with the time evolution of the orbital elements of the interacting bodies. Both files can be edited with the preferred plain text editor and the results can be plotted with the preferred graphic utilities. The code, usage details, and several examples can be found on its website [8]. This paper is organized as follows: first, in section 2 we will introduce a brief explanation on orbital elements and orbital integrators, and then in section 3 we will show how to discover the fundamental properties of planetary systems through some numerical experiments. Finally, in section 4 we present our conclusions.

Orbital elements and orbital integrators
In orbital simulations of the solar system the reference system is usually composed by the plane xy defined by the orbit of the Earth as it was in the year 2000, with the x axis pointing to the J2000.0 dynamical equinox, or Aries point. The z axis is orthogonal and directed according to the sense of the Earth's orbital revolution around the Sun. The center of the The angular orbital elements i, Ω, and ω. The spatial orientation of the orbital plane is defined by i and Ω, and the orientation of the axis of the ellipse in the orbital plane is defined by ω.
system is usually the Sun, which is not inertial due to planet attractions, but this is not a problem for representing relative motion. Of course, the equations of motion are first written in an inertial frame with the origin in the barycenter of the system, and then the origin is changed to the Sun [9]. Cartesian positions and velocities are not appropriate for the representation of orbital evolutions because they vary quickly according to the revolutions around the central star. Instead, astronomers use orbital elements, which are slowly varying parameters. Orbital elements emerge assuming that the orbit is instantaneously a conic section. The semimajor axis, a, expressed in astronomical units and eccentricity, e, define the size and shape of the orbit, respectively. - are the perihelion and aphelion, or minimum and maximum, distance to the Sun, respectively. The inclination,     i 0 180 , and longitude of the ascending node,    W  0 360 , are the two angles that define the orientation of the orbital plane with respect to the reference plane xy. The argument of the perihelion, 360 , or the longitude of the perihelion, v w = W + , define the orientation of the axis of the ellipse in the orbital plane. Finally the mean anomaly, 360 , defines the position in orbit; for example, with =  M 0 being the perihelion and =  M 180 being the aphelion. Figure 1 shows the angular orbital elements w W i, , with respect to the reference system xyz. It is possible to show that in the case of a regular evolution a, e, and i have small amplitude periodic oscillations and are related to the action canonical variables; meanwhile, v W M, , also includes variations proportional to time and related to angle canonical variables [10]. In the numerical integrations, for each object we must specify the set of initial elements , which is equivalent to , , , , , x y z , and the mass m in units of solar masses  ( ) M . Burns [11] presented a useful analytical deduction of the equations for the time evolution of orbital elements due to different kinds of perturbations when expressed analytically. Our approach is completely numerical and restricted only to gravitational perturbations, which in fact dominate the long term evolution of planetary systems.
To obtain the position and velocity at time starting from position and velocity at time t i for a particle under acceleration a  i , elementary numerical algorithms usually implemented to resolve general problems are roughly of the type = + D For the case of orbital dynamics, Encke [9] proposed a very clever alternative taking into account that the total acceleration acting on each particle, a  , is in fact the addition of a large acceleration due to the Sun, a  s , plus a small acceleration or perturbation, a  p , due to other bodies and other eventual non-gravitational effects, like the forces generated by comet emissions or the several forces due to solar radiation acting on the surfaces of the bodies [12]. The heliocentric part of the acceleration does not need to be numerically integrated because it has an analytical solution given by the well-known two-body problem. We only need to integrate the perturbation given by a  p , which is of the order of 1000 times smaller than the heliocentric acceleration, and then we can use timesteps 1000 times greater, consequentially speeding up the numerical integration. A rough scheme for Encke's method is as follows. We call S i the heliocentric, non-perturbed part of the solution defined by the initial conditions The new position at + t i 1 is given by S i evaluated at + t i 1 : However, the velocity at + t i 1 corresponds to Now, the new heliocentric solution at + t i 1 can be defined as and the procedure continues. Note that the only step that is numerically integrated is that in equation (5) because S is an analytical expression. Another very clever idea is to construct the algorithm such that in each step of the numerical integration energy is preserved. The leapfrog scheme [13], for example, is one such method for accomplishing this. These are the key ideas that allowed astronomers to investigate planetary systems by means numerical integrations of billions of years. We can be confident in these results because using different algorithms elaborated by different researchers the main results are essentially the same. The ORBE integrator is a very simplified version of the code Evorb [14], which was presented by Brunini and Gallardo and follows the aforementioned ideas, also uses Encke's method and integrates the perturbations by a leapfrog scheme with a timestep of the order of 1/40 of the smallest orbital period, which is a known good compromise between computation velocity and precision [15]. In its present state ORBE only takes into account the Newtonian gravitational interactions, and no General Relativity corrections or other forces are considered. Another limitation is that it handles events of very close encounters between two orbiting bodies with low precision.

Secular evolution
One of the most impactful results of the planetary theory developed by Lagrange and Laplace in the eighteenth century was the analytical proof of the stability of our planetary system. They showed that despite the complexity of the mathematical problem, the planetary orbital eccentricities and inclinations have only bounded amplitude oscillations while the semimajor axes are almost constant in time [10]. This can be illustrated by integrating the solar system by 1 million years (Myr) with ORBE, using the initial conditions given in table 1 and plotting   a,e, and i versus time using the output data. Two centuries after the analytical prediction we can verify it numerically in a few minutes. Figure 2 shows the time evolution of the semimajor axes of the four terrestrial planets obtained from the same simulation. No secular trends are observed, although there are high frequency oscillations with amplitudes of a few 10 −6 au that cannot be discerned in this scale but can be observed by looking at the output data file. Figure 3 shows the time evolution of the eccentricities that seem to be the result of various oscillations. The inclinations show similar behavior to eccentricities. Note that at present Earth's eccentricity is smaller to that of Mars, but this situation was not always so in the past and will not always remain this way in the future. The variations in orbital eccentricity have a known impact in the insolation that a planet receives from the Sun, which is an external cause of climate variations [16].
The time evolution of the system is not altered if we add the dwarf planets Ceres or Pluto or any asteroid, transneptunian object, or comet. If we drop Mercury instead, the outer solar system is not greatly affected, but changes in e(t) and i(t) for Venus and Earth clearly appear after~10 5 years of orbital evolution. Do the experiment by yourself.

Angular momentum
It is clear from figure 3 that the eccentricities of Earth and Venus evolve in opposed phases. This is more evident for the Jupiter-Saturn pair shown in figure 4, which used the output of the same numerical integration performed earlier. As we explain below, this is because these pairs are coupled such that the angular momentum of these subsystems tend to remain constant.
The orbital angular momentum for a planet with mass m is a vector orthogonal to its orbital plane and by definition is the vectorial product =   L r mv, which written in terms of the orbital elements becomes [10]    Since a J and a S remain constant, an increase in e J must be correlated to a decrease in e S . This effect is seen in figure 4. If this reasoning is correct and we run a numerical integration with a 'Saturn' in a retrograde orbit, that is orbiting contrary to Jupiter with~ i 180 S , we will obtain e J and e S evolving in phase. The result of this experiment is shown in figure 5. Planetary systems composed of only two planets orbiting a star always exhibit coupling between eccentricities.

Forced and free modes, fundamental frequencies, and chaos
To understand the time evolution of the eccentricities shown in figure 3 we will make a simple experiment considering only Jupiter and a fictitious particle (m=0) with initial arbitrary asteroid-like orbital elements given by     ( ) 2.3 au, 0.035, 0.5 , 90 , 270 , 0 and integrate it by 40000 years. The time evolution of its eccentricity is shown in figure 6. It seems to be a kind of sinusoid, but its true meaning is revealed when we plot the data in the sin , thus obtaining figure 7. Here, it is evident that the instantaneous, also called osculating, e(t) is the composition of a circulation with almost constant amplitude (ẽ 0.012 free ). This is called free mode, and circulates with a proper frequency around a fixed point at (  constant forced mode and a free circulation [10]. Experimenting with the same asteroid but with different perturbing planets it is possible to explore how the free and forced modes arrange to provide the same initial osculating eccentricity.  If we consider a system of perturbing planets instead of the simple system composed by Jupiter alone, we will verify that the forced mode is a composition of several vectors each with a different frequency. If we eliminate Jupiter from the solar system, the component of the forced mode due to Jupiter vanishes and Earth's eccentricity for example will have a smoother time evolution. Figure 8 shows the time evolution of Earth's eccentricity in the actual solar system and in an hypothetical solar system in which Jupiter was eliminated. A spectral analysis of the time evolution of the variables q(t) or p(t), and k(t) or h(t) corresponding to a member of the solar system will show the presence of a set of well-known frequencies, called fundamental frequencies, of the solar system and are denoted as f i for the variables (q, p) and g i for the variables (k, h). For example, in figure 9 we show the power spectrum of k(t) for the Earth. Three relevant frequencies in units of yr −1 appear with large amplitude:´-3.3 10 6 ,´-5.7 10 6 , and´-13.5 10 6 . They correspond to the fundamental frequencies g 5 due to Jupiter, g 2 due to Venus, and g 3 due to Earth [10]. There is another less relevant frequency at´-21.8 10 6 yr −1 , known as g 6 , which corresponds to Saturn. The existence of these well-defined frequencies is a strong indicator of the stability of the system. On the other hand, a planetary system with poorly defined frequencies or varying frequencies is an indicator of its chaotic nature [17]. Several sophisticated tools have been developed to diagnose chaos in the dynamics of the planetary systems. Although not the focus of this paper, we can have a hint at the magnitude of the chaos in a planetary system by looking at the time evolution of the semimajor axes or analyzing how well-defined the fundamental frequencies are [18]. In fact, our planetary system has its fundamental frequencies well defined in timescales of millions of years, but in timescales of billions of years, small variations have been detected [17], which means that in long timescales the solar system is chaotic but not necessarily evolving to disruption [19]. . This plot describes the evolution of the spatial directions of the planetary angular momenta. The trajectory for each planet is labeled with the planet's starting letter. The cross indicates the direction of the total angular momentum of the system, which is given by   i 1.6 and W   108 .

The invariable plane
As we have explained in section 2 it is a generalized practice to take reference plane xy as coincident with Earth's orbital plane for the year 2000. That choice is arbitrary from the physical point of view because the Earth's orbit is not dynamically privileged with respect to other orbital planes. There is another plane with a clear physical meaning: the invariable plane, which is defined as the plane perpendicular to the total angular momentum of the system  L T : where  L i for each planet is given by equation (7). Despite its large mass the contribution due to the Sun can be neglected because it is very close to the barycenter and moving with very low velocity. If we consider the system as isolated,  L T and the invariable plane will be constant. According to equation (7) the projection of the angular momentum unity vector for each planet on the reference plane xy is W W ( ) i i sin cos , sin sin , where we can take i instead of i sin for low inclination orbits as is the case for the solar system. If we plot the orbital states in sin for all the giant planets along 2 Myr we obtain figure 10. This figure shows that all planets have their angular momentum unity vectors oscillating around a fixed direction marked with a cross, which is defined by the  L T of the system. An analogy can be found between the oscillations of the  L of a given planet around  L T and the precessional motion of a spinning top around the vertical. The invariable plane for planetary systems has its analog in the horizontal plane for the physics laboratory. Traditionally we use the Earth's orbital plane as the reference plane, but sometimes for some particular dynamical studies the use of the invariable plane is preferred.

Orbital resonances
At present there are around 700 000 asteroids with determined orbital elements in public databases [20], but their distribution between Mars and Jupiter is far from uniform and there are profound gaps and areas of higher concentrations. This peculiar distribution in the asteroids' semimajor axes must be due to their different dynamical evolutions. Let us compute the orbital evolution of a set of 21 fictitious particles, with m i =0, distributed with initial semimajor axes between 2.4 and 2.6 au and initial e i =0.1, =  i 5 i , and arbitrary w We integrate by 1 Myr and take snapshots of the system at each 1000 yrs, plotting (a, e) for all particles. Thus, we obtain figure 11. Particles withã 2.5 au i are affected by a strong dynamical mechanism, a mean motion resonance, which means that these particles have their orbital periods commensurable with the orbital period of a planet, in this case Jupiter. The obvious effect of this resonance is to excite the eccentricities and destabilize the orbits. By Kepler's third law we know that the orbital periods T verify Then, in the case of a particle with initial = a 2.5 au we have , and it is said that the particle is in 3:1 resonance with Jupiter. This commensurability in the orbital periods generates a non-random spatial distribution of the perturbations by Jupiter making the orbital evolution of the resonant particles different from those that are non-resonant. There is a profound gap atã 2.5 au in the semimajor axes distribution of asteroids generated by this resonance, and it is not the only one [10]. The solar system is full of mean motion resonances but only some of them have dynamical relevance [21]. Some resonances generate unstable dynamics but there are also resonances that provide a very stable orbital evolution. The term resonance is usually associated with a system that is approaching disruption, but in orbital dynamics they are configurations of equilibrium that are sometimes unstable and sometimes stable. The Hildas, for example, is a stable dynamical family of asteroids trapped in 3:2 resonance with Jupiter. In general the stability depends on the degree of eccentricity growth. Large changes in eccentricity like in the resonance 3:1 make the asteroid approach other planets, which, by gravitational pull or even collisions, remove the asteroid from the resonance and thus generates the observed gap in the distribution of semimajor axes. A very interesting case of mean motion resonance are the co-orbitals, that is, objects with the same orbital period or trapped in the resonance 1:1. There are several examples of co-orbitals in the solar system, and quasi-satellites are probably the most interesting [22]. They seem to be satellites of a planet but are just revolving around the Sun with the exact same planetary orbital period, sometimes ahead of the planet and sometimes behind, generating a relative trajectory with respect to the planet that seems to be a satellite-like orbit. This happens not by chance but is one of the possible configurations of equilibrium in the resonance 1:1. There are precise methods to identify when an asteroid is locked in resonance but a useful indicator is the ratio ( ) a a P 1.5 according to equation (10), where a P is the semimajor axis of the planet and must remain constant and very close to a simple fraction over some thousand years at least.
Resonances can also occur as a commensurability between the proper frequencies of the particle and the fundamental frequencies of the planetary system given by the time evolution of the variables q(t) or p(t), and k(t) or h(t). This situation generates a secular resonance and large orbital changes usually occur [18]. Since the proper frequencies of a particle depends not only on its semimajor axis but also on its eccentricity and inclination, every planetary system has dangerous routes in the space ( ) a e i , , dominated by secular resonances, which could generate drastic orbital changes in approaching asteroids, comets, or meteoroids.

The Kozai-Lidov mechanism
Non-resonant low eccentricity and low inclination orbits in general will undergo typical secular time evolutions of e and i with small amplitude oscillations without links between them, but for moderate values of eccentricities and inclinations there appears to be a clear link between e, i, and ω. This was first studied by Lidov and Kozai [23]. Its origin is in the quasiconservation of the component L z of the angular momentum of the body, and imposes - ( ) a e i 1 cos constant 2 according to equation (7). In particular, L z is strictly constant for a particle perturbed by an ideal planetary system composed by coplanar and zero eccentricity orbits [24]. Indeed, in this idealized case the study of long-term evolution can be conducted as originally proposed by Gauss, replacing our planetary system by a system of concentric coplanar rings with each one corresponding to one planet, and with the radius and mass corresponding to the semimajor axis and planet mass, respectively. This axisymmetric distribution of mass generates a gravitational force,  F, on a body located at  r that is not central but is always coplanar with  r and the z axis. As the rate of change of the angular momentum equals the applied moment, and it follows that In the real solar system the L z of asteroids or comets is quasi-constant. This mechanism is responsible for very large excursions in e and i, and in some cases, driving to collision with the central star. Such objects are known as sungrazers [25]. Thousands of sungrazer comets are known and some asteroids also have this type of orbit. In figure 12 we show the obtained ( ) ( ) e t i t , for a fictitious particle with initial circular orbit but with =  i 70 under the perturbation of the four giant planets. Its semimajor axis remains almost constant but its eccentricity suffers large variations correlated with changes in inclination. It is convenient to follow the evolution of the perihelion distance given by -( ) a e 1 because when it reaches the value~  R 0.005 au a collision with the Sun is certain.

Conclusion
By means of numerical simulations using a physical model for a planetary system we will not discover new fundamental laws of physics but we can discover new physical properties of complex dynamical systems. These properties can be deduced by sophisticated perturbative theories but can also be observed by means of a well-conducted analysis of the output of numerical integrations. We have shown some of the basic dynamical properties of planetary systems using our solar system but these properties can also be found in extrasolar [26] and fictitious systems. Several experiments can be performed with a code like the one we have presented here in order to explore and reveal the dynamics of planets and minor bodies. Some examples can be found on the ORBE website; the initial conditions for the examples presented here can also be found there. We remark that as ORBE is only for educational purposes, serious research must be conducted using other professional integrators. A planetary system is usually imagined as a static collection of ellipses but very rich dynamics emerge when we consider their mutual interactions. Regular oscillations with well-defined frequencies or chaotic evolutions that sometimes lead to catastrophic ends are two welldifferentiated dynamical evolutions. Since the observed planetary systems are in general the result of billions of years of orbital evolution it is largely more probable to find stable planetary systems than unstable ones. Nevertheless, unstable planetary systems can be found if they are in the first steps of their orbital evolution, like those that exhibit traces of the accretion disk. It is exciting that we can explore them in the classroom with simple numerical integrators.