Using N-body simulations to investigate solar
system stability after the insertion of a planet between the orbits of Mars and
Jupiter.
Graeme Ing, HET606, 2005 S2
Project 63: Planetary Dynamics
Project Supervisor: Adam Deller
Abstract
In this paper I present a summary of the Equations of Motion as they
pertain to the orbits of the planets, and of N-Body simulation, with a
discussion of numerical errors inherent in mathematical modeling of this type.
Further, I use the scenario of a planet inserted between the orbits of Mars and
Jupiter, according to Bode’s Law, to investigate stability of a planetary
system. Multiple simulations were run, varying the eccentricity, orbital
semi-major axis and mass of the inserted planet and the resulting data analyzed
with the help of visual data plots.
1. Introduction
This paper details a series of simulations designed to investigate the
effects of introducing a planet between the orbit of Mars and Jupiter, upon the
stability or otherwise of our solar system. My basic premise is: “What if a
planet existed at the same distance from the Sun as the Asteroid Belt, in
accordance with Johann Bode’s Law?” I elected to name this object “Planet V”.
To investigate system stability, I intend to
run numerous simulations to vary the mass, orbital eccentricity and orbital
distance from the sun. These simulations of multiple planetary bodies
constitute a series of “N-body problems”. It is therefore necessary for this
paper to outline the N-Body problem and the equations of motions that will
constitute the simulations, as well as discuss the numerical errors inherent in
mathematical simulations.
I will present and discuss a
summary of the results gathered from the simulations, concentrating on how the
introduction of our mythical Planet V and the variance of its mass and orbital
properties affected solar system stability.
2. Methodology
2.1 Bode’s
Law
In 1772, Johann Bode announced a mathematical series that predicted the
distance of each planet from the Sun. Though he published his “law” before the
discovery of planets beyond Saturn or of the asteroids, it remained remarkably
accurate for Uranus and Ceres, the first asteroid to be observed.
Table 1 compares Bode’s
predicted distances from the Sun (in Astronomical Units, or AU) to actual mean
distances.
|
|
Mercury |
Venus |
Earth |
Mars |
Ceres |
Jupiter |
Saturn |
Uranus |
|
Pluto |
|
Bode |
0.4 |
0.7 |
1.0 |
1.6 |
2.8 |
5.2 |
10.0 |
19.6 |
38.8 |
77.2 |
|
Actual |
0.4 |
0.7 |
1.0 |
1.5 |
2.8 |
5.2 |
9.5 |
19.2 |
30.1 |
39.5 |
Table 1: Mean planetary distances from the Sun in AU
(EncAA-IoP)
In the simulations to follow, I shall be investigating the effects of a
planet inserted at a distance of 2.8 AU from the Sun, as predicted by Bode’s
Law.
2.2 Equations
of Motion
Equation 1 details the basic equation of motion,
Eq. 1
Eq. 2
Acceleration (a) is the rate of change of velocity over time. Equation 3a is a discrete form of this analogue relationship, more suited to a mathematical model, where Δ is a small but discrete change. Similarly, velocity (v) is the rate of change of position over time, discretised to obtain Equation 3b. (WebWolfram)
Eq. 3a 
Eq. 3b
Eq.
4
Rearranging Equation 3a and substituting in Equation 4 produces Equation
5, providing a simple definition for a small change in velocity v. Equation 6
represents a small change in position, a rearrangement of Equation 3b.
Eq. 5 
Eq. 6
These two equations form the basis for the mathematical model of planetary movements. The simulation integrates over tiny changes in time, Δt, and calculates the corresponding changes in velocity v and position r, applying these changes to the previous values to obtain new starting values for the next iteration of Δt.
One major issue with Equation 2 is that as bodies become close to one another, the force between them increases. If the distance r becomes zero then the force is infinite. The effect of this relationship upon the mathematical model is that large errors occur when bodies become very close. In fact, if Δt is large enough it is possible for two bodies to move right through one another without colliding. Because of this, Δt is said to control the stability of the system. The solution is to use a variable time step Δt, where Δt can be larger when bodies are far apart, but is suitably smaller when the bodies are close.
2.3 N-Body
simulation
Simulating the equations of motion for one or two bodies is straightforward. With three or more bodies, it is necessary to consider the forces acting upon each body in turn by every other body in the system. For example, Equation 7 states the force acting upon just one of four bodies, F1. In addition, F2, F3 and F4 must be calculated using similar equations and then applied for each time step Δt.
Eq
7
This is known as N-Body simulation. In general, there are N(N-1) calculations per time step Δt, where N is the number of bodies. (E.g. for 4 bodies, there are 12 calculations.) Equation 8 states the generic form of Equation 7.
Eq.
8
The simplest and most common N-Body simulation method is the
Particle-Particle method (Graps 2000; Ricker 2005). This method can be broken
down into four basic steps.
· Calculate all the forces upon each particle due to every other particle, using Equation 8.
· Integrate the equations of motion (Equations 5 and 6), substituting in the forces calculated in the step above.
· Increment the time counter by Δt.
· Repeat for the intended duration of the simulation
2.4 Numerical
Errors and instabilities
Mathematical simulations are prone to numerous errors that affect not
only the results, but also the stability of the simulation. It is very
important to consider the effects of these errors and seek to minimize them.
The most basic errors to consider when
performing a computerized simulation are the accuracy and resolution of the
floating-point math. Resolution is defined by the number of storage bits
allocated to the fractional part of the floating-point number. The resolution
of the exponent part is not relevant to the scale of numbers encountered in a
solar system simulation. Floating-point accuracy errors are outside the scope
of this paper except to understand that they will lead to rounding errors, even
in just a few decimal places.
I must also be cognizant of the accuracy of
input data. Expressing the Gravitational Constant to 1 decimal place is going
to result in a significantly less accurate simulation than expressing it to 4
or more places. Given the effect of rounding errors, I have chosen to limit all
input data to 5 decimal places. (The default data accuracy for the SWIFT
simulator that I shall use is 5 decimal places in any case.) After numerous test
runs, I noticed regular inaccuracies in the 3rd or 4th
decimal place, so 5 places is likely to be more than sufficient.
Other errors occur by the necessity to discretise
the equations of motion into a mathematical model. The equations of motion are
continuous or analogue relationships, which, as described in section 2.2, have
been converted into equations whose values depend upon a minute but discrete
interval of time, Δt,
referred to as the Integration Time Step.
The accuracy of the model is dependent upon the size of the integration time step.
If too large, the model lacks overall resolution; if too small, the run-time of
the simulation lengthens extremely quickly. A small integration time step is
more critical in the inner solar system, when bodies are closer. Best results
are obtained by utilizing a variable integration time step, reducing its value
as bodies become very close. It is not known whether the SWIFT simulation
software employs a variable time step, so this factor is outside my control.
It is important to record data as often as
possible to prevent simulations becoming too coarse. In SWIFT, the Output Time Step
parameter controls the frequency of data output.
The total running time of the simulation is
controlled by the Total Integration Time parameter. It is important to remember
that the longer the simulation runs, the more rounding and other numerical errors
will accumulate. Ideally, the total integration time is a large multiple of the
longest planetary orbital period, to allow each planet multiple orbits in each
run, otherwise results might be skewed by outer planets that remain
significantly apart or close together, e.g. Neptune and Uranus on the same or
opposite sides of the Sun for the entire run. There are two simple ways to
minimize this effect: a) Perform a large number of runs with random starting
positions and velocities of each planet, and b) Continue each run using the
positions and velocities from the end of the previous run.
Finally, to smooth out any rogue errors or
positional coincidences, I will make a number of runs to create a statistical
mean.
2.5 Simulations
For this project, I chose to run three different simulation
experiments, investigating the effect of distance from the Sun, in the form of
semi-major axis of the orbital ellipse (SMA), orbital eccentricity and mass of Planet
V. To investigate these factors, I chose to simulate the same planetary system
in all cases, varying a single variable in each.
First, I decided to reduce
the solar system to 6 planets, namely Venus, Earth, Mars, Planet V, Jupiter and
Saturn. My reasons for doing so were:
·
Test
simulations suggested that the effects on or from Mercury and planets beyond
Saturn were minimal. In the case of the larger planets Uranus and Neptune, I
made this assumption knowing that specific alignments of major gas giants may
have an effect on some simulation runs, but would likely have a statistically
insignificant effect.
·
Reducing
the number of planets significantly decreased simulation run times as per
N-Body simulation theory. Six planets require 30 calculations per time step,
compared to 90 for 10 planets.
·
By
excluding planets beyond ~10 AU, the graphical results provided data that are
more granular, allowing easier analysis of data plotted in the range 0.7 AU
(Venus) to 5.2 AU (Jupiter) - my chief area of interest. Using the full solar
system, graphical results would have extended to 40 AU (Pluto).
I then carefully selected the following operational parameters to use
for all simulations:
Total Integration Time:
100,000 years. For each simulation, I made 4 continuations for a total
run of 500,000 years. My assumption was that a planetary system would become
unstable in this period, though I understand that more subtle instabilities
might not show for millions of years.
Output Time Step: 50 years. I chose this value as a compromise between creating sufficient
data points yet not producing too much data for the SWIFT software to output.
SWIFT has a limitation of 2000 data points per body, yet I wanted to run each
simulation for hundreds of thousands of years.
Integration Time Step: 0.1
years. For accuracy, I wanted several time steps per Earth year. For the
purposes of this project, I did not believe it necessary to make this value
smaller.
Table 2 summarizes the initial data I used for each planet.
|
Planet |
Ecc |
S.M.A. |
Incl. |
Mass |
|
Venus |
0.00677 |
0.72333 |
3.39471 |
0.00256 |
|
Earth |
0.01671 |
1.00000 |
0.00005 |
0.00315 |
|
Mars |
0.09341 |
1.52366 |
1.85061 |
0.00034 |
|
Planet V |
Variable |
Variable |
0.00000 |
Variable |
|
Jupiter |
0.04839 |
5.20336 |
1.30530 |
1.00000 |
|
Saturn |
0.05415 |
9.53707 |
2.48446 |
0.29942 |
Table 2: Planetary Input Data
3. Simulation
Results
3.1 Determining
an unstable system
Before I could interpret my results, I had to define an unstable
planetary system. I decided upon the following:
·
A
runaway or erratic eccentricity, as demonstrated in Figure 1. An erratic
eccentricity is a good indicator of a planetary collision, ejection from the
system, or significant and unstable migration toward/away the Sun.
·
A
runaway inclination, as demonstrated in Figure 2.
·
A planet
crossing the orbit of a neighbour and moving to within 50% of the next planet,
as demonstrated in Figure 3. Whilst orbit swapping or crossing is not itself
indicative of an unstable system (c.f. Pluto and Neptune), I assumed that the
migration of an inner solar-system planet across almost two orbits was a likely
candidate for a collision or major instability over time.
·
The
complete ejection of a planet from the system, as demonstrated in Figure 4.
Any of these events would result in the classification of an unstable
system, even if the resulting planets settled into stable orbits for the
remainder of the simulation.
Figure 1: Runaway Eccentricity Figure 2: Runaway Inclination
Figure 3: Orbital crossing Figure 4: Ejected planet
3.2 Eccentricity
In this set of simulations, I placed Planet V at a distance (SMA) of
2.8 AU, according to Bode’s prediction, and chose an arbitrary mass of 0.006 MJup or about twice the mass of the Earth. Table
3 shows the percentage of stable systems after running each value 10 times.
Figure 5 plots the stability curve:
|
Eccentricity |
Stability |
|
0 |
100 % |
|
0.2 |
100 % |
|
0.25 |
100 % |
|
0.3 |
70 % |
|
0.4 |
30 % |
|
0.5 |
20 % |
|
0.6 |
10 % |
|
0.7 |
0 % |
Table 3: Eccentricity results Figure 5: Eccentricity plot
The results obtained from this set of simulations were as I expected. A
small mass in a perfectly circular orbit part way between Mars and Jupiter was
100% stable in all simulations. Instability began to occur with an eccentricity
greater than 0.25. All solar system planets possess an eccentricity
significantly less than 0.1 except for Mercury (0.2) and Pluto (0.25), so their
orbital stability is consistent with my data. As Figure 5 demonstrates, the
plot of stability against eccentricity is close to a linear relationship.
Since the mass of Planet V
was very small in comparison to that of Jupiter and Saturn, I fully expected to
see turmoil amongst the inner planets. Pre-simulation, I predicted that my data
would mostly present ejections of Planet V and/or Mars, since they were the two
planets most affected by the increase in eccentricity. Figure 6 demonstrates
the result of one simulation that met my predictions – the ejection of Planet V
from the system (cyan). Figure 7 demonstrates a Mars ejection (blue).
Figure 6: SMA plot – Ejection of Planet V Figure 7: SMA plot – Mars ejection
The second example is interesting in that it also demonstrates Planet V
(cyan) crossing inside the orbit of Mars, whilst dragging both Venus (red) and
Earth (green) outwards before finally upsetting Mars’ (blue) orbit entirely.
Figure 8 shows the
eccentricities of the planets in the same simulation as Figure 7. It demonstrates
turmoil amongst the terrestrial planets (as I predicted), and notice the
runaway eccentricity of Mars (blue) prior to its ejection. This Figure also
shows the stability of Jupiter (magenta) and Saturn (yellow), both too large to
be affected by the inner planets.
Figure 8: Ecc plot – Inner Planet Turmoil Figure 9: SMA plot – Saturn ejection
In a couple of simulations, the gas giant
planets were significantly perturbed in their orbits, much to my surprise.
Figure 9 demonstrates one such case where Earth, Mars, Planet V and Saturn were
ejected from the system. Note that the ejection of Mars (blue), pushed Saturn’s
orbit (yellow) outward and then again at the ejection of Earth (green). Finally,
Saturn appeared to achieve solar system escape velocity entirely.
3.3 Semi-Major-Axis
In this set of simulations, I retained the mass of 0.006 MJup, and locked the eccentricity of Planet V at
0.0. Starting at a distance (SMA) of 2.8 AU, I intended to research stability
as I moved the planet closer to Mars or Jupiter. Table 4 shows the percentage
of stable systems after running each value 10 times. Figure 10 plots the
stability curve:
|
SMA (AU) |
Stability |
|
1.58 |
0
% |
|
1.65 |
30
% |
|
1.80 |
50
% |
|
2.80 |
100 % |
|
3.50 |
100 % |
|
4.00 |
80
% |
|
4.05 |
40
% |
|
4.10 |
10
% |
|
4.25 |
10
% |
|
4.50 |
0
% |
Table 4: SMA results Figure 10: SMA plot
A bell curve is exactly what I would have expected from this
simulation. There is a broad orbital range midway between Mars and Jupiter
where Planet V has a stable orbit, between about ~2.5 AU and 4.0 AU. As the
planet’s SMA approaches Mars or Jupiter, its stability drops off.
My results showed that as Planet V, with a
mass of two Earths, approached Mars orbit, the system remained stable to within
~0.1 AU and in most cases, Mars was ejected from the system. This is likely due
to the small masses of both planets allowing them to approach closely, yet the
much smaller Mars was less likely to be able to retain its orbit, resulting in
its ejection more often than the more massive Planet V. At the other end of the
bell curve, Planet V was unable to approach within 1 AU of Jupiter without
being ejected from the system due to the enormous strength of Jupiter’s
gravitational attraction.
Figure 11: SMA plot – Ejection of Planet V Figure 12: SMA plot – Mars Ejection
Figure 11 demonstrates a swift ejection of the Planet V (cyan) with a
starting orbit of less than 1 AU from Jupiter. Figure 12 demonstrates erratic
behavior of Mars (blue), explained by the eccentricity of Mars becoming
excessive enough to first throw it out beyond the orbit of Saturn, bring it
back inside of Jupiter and then eject it entirely from the system. Notice also,
how the ejection of Mars destabilizes the orbits of Earth and Venus, bringing
them to within 0.1 AU of each other. Had I continued this simulation the
results were likely to be a collision or ejection of one or both planets.
Figures 13 and 14
demonstrate disturbances in Mars’ orbit due to the proximity of the SMA of Planet
V. In the first case Mars moves slowly inward to cross the orbit of Earth
whilst the in the second, Mars spirals outward to a new, stable orbit just less
than 3 AU. In both cases, the orbits of the other terrestrial planets remained
stable. I found it very intriguing to see such a smooth “orbit swap”.
Figure 13: SMA plot – Mars moves inward Figure 14: SMA plot – Mars moves outward
3.4 Mass
In the final set of simulations, I repositioned Planet V at a distance
(SMA) of 2.8 AU, with an eccentricity of 0.0 and investigated the effect of
altering planetary mass. Table 5 shows the percentage of stable systems after
running each value 10 times. Figure 15 plots the stability curve:
|
Mass MJup |
Stability |
|
1 |
90 % |
|
2 |
80 % |
|
2.5 |
60 % |
|
2.55 |
80 % |
|
2.575 |
40 % |
|
2.6 |
20 % |
|
2.75 |
20 % |
|
3 |
20 % |
|
4 |
20 % |
|
5 |
0 % |
Table 5: Mass results Figure 15: Mass plot
Results for these simulations displayed two obvious quirks: The
instability knee fell in the range ~2.5 MJup to ~2.7 MJup yet data for 2.5 MJup and 2.55 MJup was not definitive. The data tail between
2.6 MJup and 4 MJup came out flat, whereas I would have expected
to see a smoother drop off between 2.55 MJup and 2.75 MJup followed by a shallow tail out to 5 MJup. My best assessment of the cause for these
quirks is the limited number of runs (10) per case. A sample of 100 per should
create a smoother, more accurate curve.
At larger masses, most instability
resulted in the ejection of Mars or Earth from the system, for neither planet
could survive the pull of a planet several times the mass of Jupiter. Figure 16
demonstrates a runaway eccentricity for Mars (blue) and a highly erratic
eccentricity of Earth (green). Figure 17 demonstrates an ejection of both
planets. In both cases, the mass of Planet V was 4 MJup.
Figure 16: Eccentricity plot Figure 17: SMA plot
Figure 18 demonstrates an interesting result against a mass of 4 MJup. Once again, Mars and Earth are ejected, but
this time Jupiter (magenta) is pushed out to take over Saturn’s orbit (yellow),
and the latter is sent into a highly eccentric orbit that results in its total
ejection. The resulting system is a stable one, consisting of Venus, Planet V
and Jupiter.
Figure 18: SMA Plot – Jupiter moves outward Figure 19: Jupiter and Saturn Cross over
Figure 19 demonstrates a similar event, except in this scenario Jupiter
(magenta) moves outward about 1 AU and Saturn (yellow) is pulled inside of
Jupiter’s orbit to ~5 AU. The future stability of these two gas giants with
such close orbits is uncertain. Venus is the only inner planet to remain in the
system.
In several simulations with
a mass of 2.5 to 3 MJup,
Earth and Venus swapped orbits in a stable manner. Figure 20 demonstrates one
such case
.
I fully expected the
introduction of a significant mass inside of Jupiter’s orbit to destabilize the
inner planets, but I did not predict the interesting behavior of Jupiter and
Saturn in several simulations. I was also
Figure 20: SMA plot – Venus/Earth Orbital swap
surprised that every simulation I ran with a mass less than 1 MJup produced a totally stable system and that it
was even possible to maintain a stable system with a mass as high as 3 or 4 MJup. An interesting additional set of
simulations would be to move Planet V inward or outward by up to 1 AU and plot
the results. I would expect to see 0% stability with much smaller masses.
4. Conclusion
The objective of this paper was to investigate the effects of a planet
inserted between the orbits of Mars and Jupiter upon a reduced model of the
solar system. In the course of this investigation, I used the Swinburne SWIFT
simulation software to model the outcome of a variety of masses, orbital
eccentricity and distance from the Sun (orbital semi-major axis).
The results provided a solid
demonstration of how altering these parameters affected the stability of the
solar system. In most cases, the simulation results confirmed my expectations,
but a number of scenarios occurred that both surprised me and added to my
knowledge; for example, how often it is possible for two planets to swap orbits
in a stable manner. My results also demonstrated how easy it is to destabilize
a planetary system and cause the ejection of one or more planets.
During the course of the
hundreds of simulation runs, I also witnessed first hand the pseudo-random
effects caused by the initial planetary positions and velocities, and the
various numerical errors that can occur in mathematical modeling.
A suitable follow-up to this
project would involve: - Running simulations for more than 500,000 years;
developing a better result set than 10 runs per parameter change (perhaps 100) to
reduce the “random noise”; and perhaps investigating the effects of multiple
planets in the Mars-Jupiter space.
References
EncAA-IoP: Encyclopedia of
Astronomy & Astrophysics, Institute of Physics
Table 1: “Bode’s Law”, 2000
http://eaa.iop.org/full/eaa-pdf/eaa/4479.html
Freedman,
R., Kaufmann, W. 2001 “Universe” sixth edition
Graps,
A. 2000, “N-Body / Particle Simulation Methods“
http://www.amara.com/papers/nbody.html#pp
Ricker, P. 2005, “Computational Astrophysics: Lecture 4: Particles and
Gravity”
http://academic.sun.ac.za/summerschool/2005/documents/lecture_notes/ricker/ricker_lecture_4.pdf
Swinburne SWIFT Planetary Dynamics Simulator
http://astronomy.swin.edu.au/staff/maddison/astro/teaching/dynamics-how.html
WebWolfram: “World of Physics”
http://scienceworld.wolfram.com/physics/topics/DynamicsandKinematics.html