Solar System Dynamics  (Project 63)  David Conner.  Newsgroup HET602A.  Supervisor Sarah Maddison.

 

Introduction.

 

This project investigates some of the possible causes of instability in planetary systems by running simulations on the Solar System simulator program on the Swinburne ‘supercomputer’ [1].  In section 1, I look at how we can construct a computer model from ‘first principles’, and test the Swinburne supercomputer model.  In section 2, a number of simulations will be run on the Swinburne computer model to investigate instability in planetary systems, with a view to determining what factors affect stability and estimating their relative importance.  No computer model is perfect, and in section 3 I investigate some characteristics of the Swinburne model itself.

 

Section 1. Constructing a model Solar System.

 

It is possible to construct a simulation program based entirely upon Newton’s laws.  Newton’s first law of motion [2] tells us that a body continues to move in a straight line unless acted upon by a force.  Every planet is acted upon by a force - the force of gravity, and this method calculates the force of gravity on each planet due to the Sun and the other planets in the system, and its subsequent effect on their motion.  Newton’s universal law of gravitation states [2];  ‘Two bodies attract each other with a force that is directly proportional to the mass of each body and inversely proportional to the square of the distance between them.’   

 

Let       F          =          gravitational force between two objects

m1        =          mass of first object

m2        =          mass of second object

r           =          distance between them

G         =          universal constant of gravitation, 6.64 x 10-11 newton m2 / kg2

 

Then                 F = G (m1m2) / r2

 

Force is a vector.  That is, it has both size and direction.  It can be shown [3] that the gravitational force experienced by mass 1 due to mass 2 is directed towards mass 2.  Conversely, the force experienced by mass 2 due to mass 1 is directed towards mass 1. This is as described by Newton’s third law of motion [2].  The properties of vectors [3] are such that we can add two forces together to get an equivalent single resultant force equal to

force 1 + force 2, as in fig. 1.

 

 

 

 

 

 

 

 

 


fig.1.  The addition of two forces to give an equivalent single ‘resultant’ force.

 

We can do the reverse by resolving a single force into two component forces as in fig.2.  All we need to do is ensure that the force vectors can be arranged into a triangle.

 

 

 

 

 

 

 

 

 


a).                                                      b).                                                            c).

 

fig. 2. A force can be resolved into two components .

 

We now use Newton’s second law of motion [2], which states that a force acting on a body results in the body undergoing an acceleration, that is, a change in velocity:

 

force = mass x acceleration          or, in symbols            F = ma             where

 

F          =          force (in newtons)

m         =          mass of body (in kilograms)

a          =          acceleration (in metres per second squared)

 

Force and acceleration are vectors.  Acceleration is in the same direction as the applied force producing it, and it may be a change in speed, a change in direction, or both.  Finally, we will use two equations of motion [3], where

 

            Dt = time interval

u = velocity at start of time interval Dt

            v = velocity at end of time interval Dt

            s = distance travelled in time interval Dt

            a = constant acceleration during time interval Dt

 

and

 

v = u + a Dt                  …equation (1)

 

s = ½ (u + v) Dt            …equation (2) 

 

 

  1. Constructing a one-planet system.

 

We can make the following reasonable, simplifying assumptions [3] [4].

 

1)      Gravity is the only force affecting the motion of the planet, and it obeys Newton’s universal law of gravitation.

2)      The Sun is much more massive than the planet and remains stationary.

3)      The Sun is in the plane of the planet’s orbit.

 

We define a Cartesian system of co-ordinates [4] with the Sun at the origin and an x-axis and a y-axis at 90° to each other.  The co-ordinates x and y give the starting position of the planet.  Using the properties of vectors we resolve the velocity vector of the planet, u, into two component vectors parallel to the x and the y-axes.  Call these ux and uy.

 

We next calculate the force on the planet due to the Sun’s gravitational field at that position, see fig.3.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


fig. 3.  A planet: position (x,y) ,velocity(ux,uy) and Sun’s gravitational force on it, F.

 

The distance of the planet from the Sun, r, is given by r = (x2+y2)1/2  by Pythagoras’ theorem (the superscript ½ meaning ‘square root’).  The force on the planet due to the Sun’s gravitational field, F, is directed towards the Sun and has a magnitude (‘strength’) given by

 

F = GMm / r2 where

 

 G = universal constant of gravitation = 6.67 x 10-11 newton m2 / kg2

 M = mass of Sun =1.99 x 1030 kg

 m = mass of planet in kg

 r   = distance of planet from Sun in metres

 

To calculate the x and y components of force, Fx and Fy, we calculate angle q (theta) (fig. 4) using the inverse trigonometrical function tan-1, where q = tan-1 (y/x), to find q.

 

Fx = F cosq and   Fy = F sinq

 

Using Newton’s second law of motion we now calculate the acceleration of the planet in the x and y directions, ax and ay caused by the Sun’s gravitational field:

 

ax = Fx/m          and       ay = Fy/m

 

Equation (1) now gives the components of the velocity at the end of a short time interval, Dt.

vx = ux + ax Dt       and         vy = uy + ay Dt

 

Equation (2) gives us sx and sy, the distances travelled in the x and y directions during Dt.

 

sx = ½ (ux + vx)            and           sy = ½ (uy + vy)

 

Hence, the new co-ordinates of the planet at the end of time interval Dt are

new (x co-ordinate)=old (x co-ordinate)+sx  and new (y co-ordinate)=old (y co-ordinate)+sy    These become the starting co-ordinates for the next time-step.  This process is repeated to give us the path of the planet, and is a simple version of the Runge-Kutta method [5] of calculating planetary orbits. (The specific method used here is essentially the modified Euler method [5].)   The importance of choosing a suitable value for Dt is explored in section 3.

 

  1. Upgrading to an n-planet system.

In a multi-planet system, orbits can be inclined to each other, that is, they may lie in more than one plane.  We therefore need to consider forces and motion in three directions, x, y and z, which are all at 90 degrees to each other, i.e. a standard 3-dimensional Cartesian co-ordinate system.  We first calculate the components of the gravitational force, Fpq, on planet p due to planet q.  Fig. 4 illustrates the position of p relative to q, as it is this that determines the size and direction of the force on p due to q.  The calculation of d is now in 3 dimensions:

 

By applying Pythagoras’ theorem we find that d = [(xp- xq)2 + (yp-yq)2 + (zp- zq)2]1/2 

 

Also, distance qr = d cos b, so cos a = (xq – xp) / d cos b    and     sin a = (yq – yp) / d cos b

 

where cos b = [(xp – xq)2 +(yp – yq)2]1/2 / d  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


fig. 4.  The gravitational force on planet p due to planet q.

We perform this calculation to determine the force on each planet due to each of the other planets, as well as the Sun.  These are added together to give the net force acting on each planet.  Once we have done this it is a matter of proceeding as in the previous section – remembering that we are now working in 3D instead of 2D.  

 

Problems with this method.

It is not possible – even in principle – to mathematically solve Newton’s laws to give us the positions of the planets at any particular time in all but the simplest systems [6] - there is simply not enough information available.  This is called ‘the many-body problem’ [7] [8].  We therefore proceed in the stepwise fashion outlined above to find what are called numerical solutions to our equations.  We must decide on the value of each variable (position, velocity) for each planet at the start of the simulation, and calculate new values at the ends of a small time interval, Dt.  These are used as new input data for the equations and we repeat the process for the next time interval, and so on, for as often as is required.  This may be hundreds, thousands, or millions of times, depending on the length of the simulation.

 

The equations of motion used assume constant acceleration during a given time interval but, as we see from Newton’s law of universal gravitation, the strength of the gravitational force varies with distance from the gravitating body.  Another error arises from assuming that the direction of the gravitational force throughout the interval is that at the start of the interval.  This is not necessarily the best average for the time interval, as can be seen from fig. 5.  A better average direction would be that at the middle of the interval.

 

 

 

 

 

 

 

 

 

 


 

 

 

 


fig. 5.  The force at middle of time interval as an average for the whole interval.

 

This is an application of the modified Euler method [5], where we calculate the acceleration at the mid point of a time interval and then apply it to the whole interval.

 

The number of calculations required in a simulation increases as the square of the number of the planets in the system [9].  Say we have n planets.  We need to calculate the forces on each planet due to the other (n-1) planets, and the one due to the Sun as well, a total of n forces on each planet.  This needs to be done for each of the n planets, so the total number of forces which need to be calculated is ((n –1 ) + 1) x n = n2.

 

This has a profound effect on the time required to perform the calculations for each time interval, and can tempt us to use longer time intervals to compensate.  But, we are not at liberty to freely select a value for the time increment, Dt.  As we will see in section 3, the success of this method depends critically on a suitable selection.  Too small, and the simulation would take an unrealistically long time to produce a result; too long, and the results are meaningless.  More sophisticated methods use a variable length time-step [10], to take advantage of the variable velocity of a planet around its orbit (the timesteps are adjusted so that the planet travels equal distances in each timestep).

 

Our method treats planets and the Sun as point masses.  That is, all the mass of an object is assumed to be at the exact centre of the body.  This has no significant effect if the distances between the objects are ‘large’ compared to their sizes and we can safely ignore tidal effects. 

 

Another characteristic of point masses is that they cannot collide!  It is possible for two points to approach each other far closer than real planets ever could. (See simulation No 13)

 

There is nothing ‘directional’ about time in any of Newton’s laws and, by giving the appropriate velocities negative signs (to indicate motion in the opposite direction), we can calculate motion backward in time just as easily as forward.  This would be useful, for example, for dating historic events that were recorded with reference to astronomical events.

 

An investigation of Swinburne’s supercomputer simulator.  

 

From running the simulations for this project, it appeared that at least three major simplifications have been made in writing this simulator:

 

1)      The planetary orbits are all in the same plane. 

 

2)      The major axes of all the orbits are on the same line, i.e. parallel to each other.

 

3)      The simulation assumes that all of the planets are at either perihelion or aphelion at the start of the simulation, that is, on a straight line.   (E-mail correspondence with my supervisor confirmed that this was to give the simulation time to run before the planets significantly interacted with each other gravitationally.)

 

None of these assumptions is valid for the real Solar System [2][11].

 

This model requires that we input values for the semi-major axis and the eccentricity of orbits.  The program then calculates the critical initial velocity of the planet that will give an orbit these values - we are automatically putting the planet into orbit!  If we had to input a value for velocity, there would be no guarantee that the planet would be in any sort of orbit. 

 

It also only allows short simulation times.  We can set the time-step to be as much as 1/20 of the inner planet’s period, and run the simulation for up to 100 000 of these time-steps.  However, if the inner planet is Mercury (period = 88 days), then the simulation time is only about 1205 years.  To investigate instabilities over such a short time scale requires the initial conditions to be totally unrealistic – the real Solar System could not be stable enough to produce them in the first place!  We can only investigate very unstable systems!

           

 

 

 

 

Section 2. Investigating instabilities in planetary systems using the Swinburne supercomputer model.

 

What factors might affect the stability of a system?

 

For the purposes of this project, I will define an ‘instability’ as any major change in an orbit (or orbits) during a simulation.  This may result in the ejection of one or more planets from the system, or it may simply result in a planet leaving one unstable ‘orbit and entering another orbit.  (In a similar way, a pencil stood on end is unstable, but it falls over and reaches a stable state.)

 

There are a number of factors that may affect the stability of a planetary system, including relative sizes of orbits, orbital eccentricity, relative orbital periods and planetary masses.  Each of these factors affects the net gravitational force experienced by a planet by varying the direction and strength of the individual gravitational forces acting upon it due to the other bodies in the system.  It is the strength and direction of the net force on a body that determines the effect on a body’s motion. 

 

We can now try to answer some questions about planetary dynamics.

 

Q) Is a planet perturbed more by another planet that is closer to, or further from, the Sun?

 

To test this, a massive planet was simulated 0.1AU inside the Earth’s orbit, and again 0.1AU outside the Earth’s orbit  The following parameters were input into to the simulator;

 

Other planet closer to the Sun.

 

Simulation number 1

Number of time-steps = 50000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.387                           0.206              

            2                      4000                            0.9                               0

            3                      1.000                           1.000                           0.107

 

Other planet further away from the Sun.

 

Simulation number 2

Number of time-steps = 50000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.387                           0.206              

            2                      1.000                           1.000                           0.017

            3                      4000                            1.1                               0

 

 

 

 

Simulation 1

The Earth (blue) is pulled into an unstable inner orbit for a number of revolutions before being put into an outer orbit again.  Its larger semi-major axis shows that it has gained some energy from the other planet, but not enough to be ejected.  (See section 3 for a discussion on energy and planetary orbits.)

 

Simulation 2

 

In this situation, the Earth (now green) gains so much energy that it ejected from the system. 

 

A) This simulation suggests that planet further away from the Sun is more likely to produce instabilities than one closer to the Sun.  This is possibly caused by the outer planet opposing the inward pull of the Sun’s gravitational field.  The affected planet will have most of its original velocity and will therefore be travelling too fast around what is, essentially, a reduced mass Sun to remain in a stable orbit.  See section 3 for a more detailed look at this phenomenon

 

[Note the use of a massless ‘Mercury’ as the inner planet.  Without this technique, one of the planets under investigation would itself be the inner planet, and could only move a maximum of 500 steps per orbit.  This would be too few for our purpose.  Again, see section 3.]

 

Q) Are similar results obtained if we increase the eccentricity of the massive planet?

 

Re-running simulation 2, with the eccentricity of planet 3 increased to 0.1 gives:

 

Simulation number 3

Number of time-steps = 100000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.387                           0.206              

            2                      1.000                           1.000                           0.017  

            3                      4000                            1.1                               0.1

           

 

A)    Yes.  (Running the simulation for longer we see that the Earth is eventually ejected)

 

 

 

Q) How does altering the mass of the perturbing planet affect the perturbations?

 

By altering the mass of the massive planet to 500 Earth masses, I obtain the following result.

 

Simulation 4

 

A)    Perturbations are smaller with a reduced mass planet (obvious, but it needed proving). 

 

Q) Keeping the product of the masses of the two planets constant, does altering the relative masses have any effect on stability?

 

A number of simulations were run with the following planetary masses (relative to Earth).

 

inner planet       outer planet

0.25          4000

1          1000

2          500

                                                4          250

5                    200

10        100

                                                20        50       

31.62   31.62

 

These were all remarkably similar, and stable, even the first system where the ratio of the masses is at a maximum, 1:16000.  This is shown in simulation number 5.

 

Simulation number 5.

Number of time-steps = 100000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.1                               0.206              

            2                      0.25                             1.000                           0

            3                      4000                            1.1                               0         

           

                   

 

Typical of the other ratios is the 1:62.5 ratio (masses 4 and 250 Earth masses respectively).

Simulation number 6 shows this:

                      

A)    The ratio of the planetary masses has relatively little effect on the stability of the system.

 

Q) Can we determine which has the greater effect, varying the mass or the distance of the perturbing planet?

A number of simulations were run to compare the relative effects of mass and distance.  The following two simulations are typical.

 

Simulation number 7

Number of time-steps = 4000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.387                           0.206  

            2                      1.000                           1.000                           0

            3                      4000                            1.1                               0

 

 

Simulation number 8

Number of time-steps = 4000

Number of inner orbit time-steps = 500

___________________________________________________________

        planet       mass / mass of Earth     semi-major axis / AU       eccentricity

            1                      0                                  0.387                           0.206  

            2                      1.000                           1.000                           0

            3                      1000                            1.05                             0

           

 

 

 

In each case, the magnitude of the gravitational force on the Earth of the massive planet is the same at the point of closest approach (from Newton’s law of universal gravitation).  However, the relative position of the two planets will also affect the direction of this force.  Kepler’s third law tells us that the Earth will overtake the other planet quicker in the first case because the planet is further out from the Earth and therefore moving more slowly. 

 

The more distant of the two massive planets will be more ‘to the side’ of the Earth’s track for any given separation, d (see fig. 6), and so tend to alter the direction of motion of the Earth more than a less massive planet closer in.  This effect is increased by it having the greater mass, and therefore the greater force on the planet, at any given distance.