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

m₁ = mass of first object

m₂ = mass of second object

r = distance between them

G = universal constant of gravitation, 6.64 x 10^-11 newton m² / kg²

Then F = G (m₁m₂) / r²

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 u_x and u_y.

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(u_x,u_y) and Sun’s gravitational force on it, F.

The distance of the planet from the Sun, r, is given by r = (x²+y²)^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 / r² where

G = universal constant of gravitation = 6.67 x 10^-11 newton m² / kg²

M = mass of Sun =1.99 x 10³⁰ kg

m = mass of planet in kg

r = distance of planet from Sun in metres

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

F_x = F cosq and F_y = F sinq

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

a_x = F_x/m and a_y = F_y/m

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

v_x = u_x + a_x Dt and v_y = u_y + a_y Dt

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

s_x = ½ (u_x + v_x) and s_y = ½ (u_y + v_y)

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

new (x co-ordinate)=old (x co-ordinate)+s_x and new (y co-ordinate)=old (y co-ordinate)+s_y 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.

2. 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, F_pq, 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 = [(x_p- x_q)² + (y_p-y_q)² + (z_p- z_q)²]^1/2

Also, distance qr = d cos b, so cos a = (x_q – x_p) / d cos b and sin a = (y_q – y_p) / d cos b

where cos b = [(x_p – x_q)² +(y_p – y_q)²]^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 = n².

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.

			closer orbit
	Massive planet = strong force off to side of track

fig. 6. The direction of pull of another planet: angle a is greater than angle b

A) A qualified ‘yes’. On balance, these effects lead to the first system being the less stable – the more massive planet at the greater distance is more destabilising.

[On a personal note, I found this result to be totally counter-intuitive. The effect of the changing direction of the force serves a reminder that Newton’s law of universal gravitation is a vector equation – such is the value of running simulations.]

Q) What effect does altering the eccentricity of the perturbing planet have?

I run simulation 9 with the inner massive planet eccentricity set to 0.2

Simulation number 9.

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.2

3 1.000 1.000 0.017

This is typical of several simulations with eccentric orbits - the blue orbit is the perturbed planet. Altering the eccentricity of an orbit can be very effective in bringing two planets close together and hence producing conditions that lead to significant perturbations.

A) Eccentric orbits are very likely to lead to instabilities.

The final factor to be considered here is resonances, that is, planets with periods that are simple ratios of each other.

Q) Can we say anything about which resonances are stable, and which are unstable?

For two planets, 1 and 2, with periods p₁ and p₂, and semi-major axes a₁ and a₂, Kepler’s third law states that

(p₁)² / (p₂)² = (a₁)³ / (a₂)³

This can be re-arranged to give a₂ = [(a₁)³ x (p₂)² / (p₁)² ]

We use the Earth as before, so we set a₁ = 1 AU and p₁ = 1 year.

To investigate the 2:1 resonance.

We set p₂ = 2 years which gives a₂ = 1.587401 AU.

Simulation number 10.

Number of time-steps = 10000

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.587401 0

The Earth’s orbit is severely disrupted, so a 2:1 resonance is not stable.

3:2 resonance.

Setting p₂ = 1.5 gives a₂ = 1.310371

In an attempt to reproduce the 3:2 stability, many simulations were run. The following simulation is typical, and suggests that this resonance is not stable.

Simulation number 11.

Number of time-steps = 100000

Number of inner orbit time-steps = 100

___________________________________________________________

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

1 0 0.02 0

2 0.001 0.9 0

3 0.001 0.95 0

4 0.001 1.0 0

5 0.001 1.05 0

6 0.001 1.1 0

7 400 1.310371 0

Here, the fourth planet at 1AU (in cyan), is in a 3:2 resonance with the planet at 1.310310 AU (in orange). It is not in a stable orbit.

This result does not appear to be consistent with information posted on the newsgroup by our tutor. Here it was stated that the Kirkwood Gap in the asteroid belt is due to a 2:1 resonance with Jupiter, while the Hildas family of asteroids resides at the 3:2 resonance (see also [12]).

3:1 resonance.

Setting p₂ = 3 gives a₂ = 2.080084

This is stable and is not shown here.

4:1 resonance.

Setting p₂ = 4 gives a₂ = 2.519842

This is stable and is not shown here.

4:3 resonance.

Setting p₂ = 1.333333 gives a₂ = 1.211414

Simulation number 12.

Number of time-steps = 10000

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.211414 0

This resonance is unstable in this system, the planet being ejected by the slingshot effect [13].

A) These simulations suggest that the 2:1, 3:2 and 4:3 resonances are un-stable, while the 3:1 and 4:1 resonances are stable.

In a resonant orbit, a planet is given a gravitational ‘push’ at the same part of each orbit. Individually, these don’t need to be strong, but the timing is important – rather like pushing a child on a swing. Gentle pushes, administered at the right time, can have a dramatic effect on the amplitude of the swing, whereas large pushes at the wrong times will not.

The simulator is not well suited to investigating resonances, probably because the simulation times are too short. Instabilities caused by resonances are not easy to produce accurately due to the over-riding need to use unrealistic masses and eccentricities.

Finally - is the following a real, physical, collision? (If so, wouldn’t the planets have disintegrated?) Or is it two mathematical points converging, that we should not interpret in terms of real planets? There is no easy way of telling…

Simulation number 13.

Number of time-steps = 7000

Number of inner orbit time-steps = 500

___________________________________________________________

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

1 0.1 0.387 0.206

2 4000 0.9 0

3 1.000 1.000 0.017

Summary of my findings.

The causes of instability I have investigated are listed in approximate order of decreasing effectiveness at producing instabilities in a planetary system:

Close approach by another planet – in this simulator this is often as a result of its orbit being given a large eccentricity at the start of the simulation.

The product of the masses of two interacting planets needs to be large.

The relative distance of the planets from the Sun is important – a planet is affected more by an outer planet than an inner planet (other factors remaining constant).

The ratio of the masses of the interacting planets is less important than their combined mass.

Resonance between orbits is not very significant in the short term

Section 3. Instabilities caused by the Swinburne model itself.

The stability of systems run on the Swinburne computer depends not only upon the physical system we are trying to model, but also on the model itself - in particular, length of time-step. The following simulations of the inner Solar System illustrate this:

Simulation number 14.

Number of time-steps = 100000

Number of inner orbit time-steps = 40

___________________________________________________________

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

1 0.055 0.387 0.206

2 0.815 0.723 0.007

3 1.000 1.000 0.017

The orbit of Mercury (red) is too eccentric, as is shown by the perihelion approaching the Sun and ‘colouring in the circle’. It is also rotating, or precessing, much faster than in reality. (Mercury’s orbits the Sun once every 88 Earth days. Dividing this into 40 steps gives each time-step a duration of 2.2 days. 100000 of these steps equals a little over 600 Earth years.) Using twice as many steps of half the length (i.e. to give the same total time) gives an inconsistent result, shown in simulation 15.

Here, it is the orbit of Venus that is very eccentric and precessing.

The system we are trying to simulate is exactly the same in both cases, yet the outcome of the simulation is affected by our choice of time-step length, Dt. These apparent differences are an artifact of the simulator itself. They cannot both be correct. To see why, we first consider the principle of the conservation of energy.

How does ‘conservation of energy’ relate to planetary motion?

The total energy of a planet in a stable orbit is constant [4]. This total energy consists of the sum of the planet’s kinetic energy and its gravitational potential energy. The kinetic energy of a planet (or any object) of a given mass is given by the expression

kinetic energy = ½ x mass x velocity²

It is not usual for a planet’s mass to vary significantly, and in this simulator they do not vary at all, so its kinetic energy depends upon the planet’s velocity only.

The gravitational potential energy of a planet is determined primarily by its position in the gravitational field of the Sun. We can picture the gravitational field of the Sun as a ‘potential well’, as in fig. 7.

fig. 7. Representing the Sun’s gravitational field by dent in a rubber sheet.

Imagine the Sun to be represented by a large ball bearing on a sheet of rubber. The depression in the rubber represents the Sun’s gravitational field. Two stationary planets of equal mass are represented by two smaller ball bearings on the sheet. If they were released from rest, they would roll down the slope as if they had been attracted by the ‘gravity’ of the central ball bearing.

‘Planet’ 1 is closer to the Sun than is ‘planet’ 2 and therefore experiences a greater attractive force due to the ‘Sun’s’ gravitational field. This is modelled by the rubber sheet being steeper at the position of planet 1 than at the position of planet 2. In other words, it feels a stronger force causing it to roll down the slope.

In order to move away from the Sun, planet 1 would need to ‘climb out’ of the depression in the rubber sheet. In doing so it would gain potential energy. It is lower down the slope than planet 2, indicating that it starts with less potential energy than planet 2 and needs to do more work in order to climb all the way out.

For a planet to orbit the Sun, it must roll around the inside of the Sun’s gravitational potential well rather in the manner of the ball in a game of roulette. It still experiences a force in towards the centre of the potential well, but it has a sideways motion that prevents it from falling directly in towards it. It remains in a state of ‘free fall’, forever travelling in a curved path around the Sun.

If we look down onto the potential well we might see the following:

a

a

	‘planet’ 1

fig. 8. Looking down onto the rubber sheet, showing the paths of two ball bearings, the arrows represent their velocity vectors.

In our analogy, the ball bearing representing planet 1 follows a circular path. Its velocity, represented by the arrow, will remain of constant magnitude as it travels around its orbit. So its speed and kinetic energy remain constant. Its distance from the Sun is also constant, so its potential energy also remains constant. Its total energy must therefore also remain constant.

Planet 2 follows an elliptical path. It is travelling fastest at a, slower at b, and slowest at c. It therefore has its greatest kinetic energy at a, and least at c. However, it has least potential energy at a, where it is closest to the Sun, and most at c, where it is furthest away. In rolling up the slope, it slows down and looses kinetic energy. Kinetic energy is converted into potential energy in such a way that the sum of kinetic energy and potential energy remains constant – a manifestation of the principle of conservation of energy [4].

Note that planet 2 at position a is travelling much faster than planet 1 is at a. They are at the same distance from the Sun, and so have the same potential energy, but planet 2 has more kinetic energy and can climb further out of the potential well than can planet 1. This results in the elliptical path. As it climbs out of the well, it gains potential energy and loses kinetic energy. As if travel further round and starts to fall back in towards the Sun it travels faster again, now gaining kinetic energy and losing potential energy. Again, the total energy remains constant. This variation of speed with distance is given by Kepler’s second law [2].

The more total energy planet 2 has, then the further it must climb out of the potential well to convert its kinetic energy into potential energy and the more elongated, or eccentric, the ellipse becomes. It is possible that the planet has so much energy that the well is not deep enough to hold it, and it escapes. This is the equivalent of reaching the ‘escape velocity’ for the Sun. The orbit would then be in ‘infinitely’ long ellipse, or parabola.

What does Dt have to do with conservation of energy?

In the simple model of part 1 the equations we used assumed that velocity and acceleration were constant. This is not the case in the real Solar System but, if we look at shorter and shorter time intervals, then this approximation becomes more useful, see fig.9.

fig. 9. Straight line approximation to an elliptical orbit. The elliptical orbit (solid line) is approximated by a series of line segments. Note that the second ‘orbit’ (dashed arrows), does not exactly follow the first orbit (continuous arrows).

The length of the time interval that we use is determined by the ‘inner orbit steps’ parameter [1]. Mercury’s period of revolution around the Sun is 88 days. Dividing this into 500 time-steps gives each step a duration of 4.224 hours. At an average speed of 48 km/sec. Mercury moves 730 000km. in each timestep – which is therefore the length of each line segment.

More sophisticated simulators might ensure that all the calculated line segments were tangents to the elliptical orbit. By doing so, they ensure that the total energy of the planet remains reasonably constant over time [6]. That is, it remains in a ‘stable’ orbit. The Swinburne simulator does not appear to do this. Running the simulator with different time-steps produces a semi-stable system or one that ‘explodes’, never one that ‘implodes’.

Why does the Swinburne model never implode?

In part 1 we saw that increasing distance represents a gain in potential energy. This must come from the conversion of an equal amount of kinetic energy – the planet must slow down by an appropriate amount. For a real planet to be ejected it must gain energy from another planet via a gravitational interaction. However, this model ejects planets from single planet systems – see simulation No 16 (page 27). Where is the planet getting this ‘energy’ from?

Assuming that the simulator calculates accurate initial velocities for the start of the simulation (based on our values of a and e) then it is being gained as the simulator runs, it does not make corrections to keep the total energy constant (as in fig. 10). At any instant of time, the planet is moving in the direction of the tangent to its position in its orbit (fig.10). The simple method in part one assumes this to be the true path for a time Dt, (fig 10b). The Runge-Kutta method improves on this by making certain corrections, but its improved estimate for the direction is still directed towards the outside of the true orbit (fig 10c).

fig. 10. Straight-line approximations point to the outside of the orbit.

The planet follows its real orbit from ‘p’ to ‘a’. Our simple model approximates this by the straight line segment ‘p-b’. This places the planet outside its actual orbit. Improved models approximate the true path with line segments such as ‘p-c’. This is closer, but still outside the actual orbit. So, whatever kinetic energy the planet may have had while traversing the straight-line segment, its new position is always further away from the Sun than it should be. Its assumed constant velocity means that it has not slowed down the correct amount to compensate for this gain in potential energy. It is therefore travelling too fast for its new position, and recedes from the Sun. It may eventually be lost from the system. This is illustrated in fig.11.

fig. 11. A planet moving away from the Sun as a result of errors introduced by approximating its elliptical orbit with straight-line segments.

It can be seen that the net result of these accumulated errors is that the planet spirals away from the Sun – as is often seen in the Swinburne simulator. This effect is greater for small orbits, with their large curvature, and less for large orbits, which are closer to being straight lines themselves. This effect is seen most readily in the orbit of Mercury. If we don’t use a hypothetical planet inside the orbit of Mercury, then it is not possible to divide Mercury’s orbit up into more than 500 segments. Mercury’s orbit is also very curved. Each of these increases the error in the distance from the Sun at each step, see fig.12.

Errors produced at each step.

Errors accumulating around an orbit.

fig. 12. Small orbits produce large errors, large orbits produce small errors.

Another consequence of using a straight-line approximation with a constant acceleration is that the direction of the accelerating force is not accurate for the whole time-interval. If we assume the value of acceleration at the start of the interval applies throughout the rest of the interval, then we produce an error that behaves a little like tidal acceleration, see fig.13.

fig. 13. Spurious tidal force produced by errors in the direction of the calculated force.

‘Tidal acceleration’ is inconsistent with point masses. It is demonstrated in the following simulation of Mercury. Note that Mercury is speeding up as it moves away from the Sun.

Simulation number 16.

Number of time-steps = 10100

Number of inner orbit time-steps = 20

___________________________________________________________

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

1 0.055 0.387 0.206

When run for 100000 steps, Mercury has travelled to a distance of 730000 AU.

(An aside; this diagram is not very clear. The simulator does have a black and white printout facility, but it ‘joins the dots’ and would not show this effect.)

These spurious results defy the principle of conservation of energy [4].

There were other numerical errors caused by the model itself, such as errors due to choice of run-time. They are also numerical instabilities, but are outside the main scope of this project. Examples have been included as appendix 1 for the sake of completeness.

An error (?) with the program itself is illustrated in appendix 2. The inner orbit is not always divided up into the chosen number of time-steps. ‘Inner orbit steps’ is not an independent variable for us to choose freely. We do not necessarily know its actual value in our simulations.

Instabilities – physical or numerical?

When real, physical, planets interact gravitationally; the total energy of the planetary system is conserved. Whatever perturbations, collisions or tidal effects there may be, energy is simply re-distributed among the planets. Whole branches of mechanics [14] have been developed which directly calculate energy (Lagrangian mechanics) and momentum - the product of mass and velocity - (Hamiltonian mechanics), rather than position and velocity. They inherently conserve energy, and are much more successful in preventing simulations from running out of control due to numerical instabilities of the type we have been discussing. Instabilities that they display are much more likely to be instabilities in the real planetary system rather than artifacts of a necessarily less than perfect numerical computer model.

In conclusion…

Applying the principles discussed in parts 2 and 3, it is not too difficult to construct unstable systems at will. Some instabilities will be physical, others will be numerical. (As implied throughout this project it is, strictly speaking, the stability of the computer model that we are investigating.) The principles we have discussed would apply to the real Solar System, but not necessarily to the extent suggested by the model

Using my findings, an exploding planetary system – where all the planets are ejected - was constructed. The system, initially the size of the orbit of Mars, expands to over a light year in diameter. Note that:

1) The orbit of the innermost planet is divided into only 20 time-steps. This is where most of the ‘energy’ for the explosion comes from (i.e. this is a numerical instability).

2) 5 of the orbits are eccentric, some extremely so. This distributes that extra ‘energy’ among all of the planets.

3) There are two massive outer planets, destabilising the inner orbits even more.

4) Some of the products of planetary masses are large (those involving planets 9 and 8 and, to a lesser extent, 1), giving some large individual perturbations.

5) Planets 4 and 9 are initially in the unstable 2:1 resonance.

6) The semi-major axis of planet 7 is critical. The intended value of 1.3 results in ‘only’ two ejected planets. The extra 0.02 AU (the result of a typographical error!) pushes the system over the edge into total chaos.

But that’s another story...

Simulation number 17.

Number of time-steps = 100000

Number of inner orbit time-steps = 20

___________________________________________________________

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

1 10 0.7 0.8

2 1 0.8 0

3 1 0.9 0

4 1 1.0 0.5

5 1 1.1 0

6 1 1.2 0.5

7 1 1.32 0.1

8 1500 1.4 0.8 9 4000 1.5874011 0

References.

[1] Numerical Experiments in Planetary Dynamics. http://www.swin.edu.au/astronomy/dynamics/

[2] W. J. Kaufmann III and R. A. Freedman, Universe, (Freeman and Co., New York, 5^th edition 1999).

[3] Marcelo Alonso and Edward Finn, Physics, (Addison-Wesley Publishing, 1970)

[4] M. Zeilik and S. A. Gregory, Introductory Astronomy and Astrophysics (fourth edition, 1998), (Saunders College Publishing).

[5] Module 5. Particle Applications - Pipeline Computing http://www.npac.syr.edu/EDUCATION/PUB/hpfe/module5/index.html

[6] XSTAR. An Introduction to the N-Body Problem. http://www.midwestcs.com/xstar/documentation.html

[7] The restricted three body problem. http://www.physics.cornell.edu/sethna/teaching/sss/jupiter/Web/Rest3Bdy.htm

[8] Dr. Robert Maddison, A Dictionary of Astronomy, (Hamlyn Publishing, 1980).

[9] N-Body / Particle Simulation Methods http://www.maths.napier.ac.uk/gavin/nbody.html

[10] Three Bodies in Gravitation.

http://astro.u-strasbg.fr/~koppen/body/ThreeBodyHelp.html

[11] The Handbook of the British Astronomical Association 2001 (ISSN 0068-130-X)

[12] News on HILDA ASTEROIDS. http://www.rzuser.uni-heidelberg.de/~s24/hilda.htm

[13] The Slingshot Effect

http://www.dur.ac.uk/~dma0rcj/SL/summary.html

[14] Lagrangian and Hamiltonian mechanics -- A short introduction http://alamos.math.arizona.edu/~rychlik/557-dir/mechanics/

Appendix 1.

The following system was run for two different runtimes, the results were mutually exclusive.

Simulation number 18.

Number of time-steps = 100000

Number of inner orbit time-steps = 50

___________________________________________________________

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

1 0.055 0.387 0.206

2 0.815 0.723 0.007

3 1.000 1.000 0.017

4 0.107 1.524 0.093

5 317.83 5.203 0.048

6 95.16 9.554 0.056

Increasing the total number of steps to 51 sends Mercury out at almost 90 degrees to this direction, as in the following simulation, No. 19.

The cause of this type of anomaly is obscure.

A possible clue comes from running short simulations, where Mercury would be expected to complete only one or two complete orbits. The simulator input suggests that we can choose a specific number of inner orbit steps and, independently, run the simulator for a specified number of those steps. However, it appears that these two variables are not truly independent, choosing one actually alters the value of the other in a way that we cannot control. We therefore do not actually know their value accurately, and may be making invalid assumptions about them.

This is demonstrated in appendix 2.

Appendix 2.

The input variable ‘inner orbit steps’ is supposed to determine how many steps the innermost orbit is divided up into. However, setting this variable to 50, and the ‘total number of steps’ to 50 produced the following plot:

This should have had Mercury completing one orbit, instead of just 2/3 of one. Reducing the inner orbit steps to 20 produced the following: just one orbit instead of the expected 2 ½ .