Precession of Planetary Orbits

Student Projects in AP1603 investigating classical motion.

All projects must have…
Name:____________
Email:_____@columbia.edu

Introduction

During our classroom discussions, we've investigated the motion of planets orbiting a massive star. In particular, we've shown how Mathematica can be used to visualize elliptic orbits, to demonstrate the conservation of angular momentum, and to prove Kepler's "harmonic law" of planets. Kepler's laws result from the inverse-square dependence of the gravitational force on the separation between massive bodies. What happens when the force does not have a perfect inverse square dependence?

In this research project, you are to investigate the stability of planetary motion when the radial-dependence of the gravitational force  changes from the ideal 1/(r^2). Nature produces these perturbations in two ways. First, the each planet accelerates according to a total force: the vector-sum of the forces from the Sun (very large) and the other planets (much smaller, but still important). The motion of three or more interacting objects can not be integrated analytically, and computer simulation of the dynamics of three (or more) objects shows the possibility of orbit instability and chaos. Secondly, near extremely massive objects, Newton's law of gravity becomes inaccurate, and we need to include corrections due to Einstein's theory of general relativity. When the gravitational force deviates from a perfect inverse-square dependence, elliptical orbits precess.

The following  two sections present related suggestions for your study and research using Mathematica. Take one of these as a starting-point, and conduct "computational experiments" and analysis to explore how elliptical orbits precess. Re-examine the notebooks discussed in class, copy relevant expressions, and modify them during your investigations. 

Your goal is to produce a complete Mathematica notebook containing a summary of your investigations. This will most likely be an answer to one or more questions concerning planetary orbits. It should contain Mathematica expressions and graphics which illustrate your solution. Your notebook should be concise and easy to read. Try not to include a large number of repeated expressions. Instead, generate a table or graphic of your results. In all cases, format your notebook and include textual comments and descriptions. Your notebook need not be long. It should be interesting and arrive at a definite conclusion.

Modification of the Inverse-Square Law

The inverse-square dependence of the gravitational force (and the Coulomb force between two charged objects) is related to  the (r^2) radial dependence of the area of a surface enclosing one object. The force of gravity emanates from an object, and the force is directed towards that object. The total gravitational "field lines" or "lines of force" are conserved as the enclosing surface moves to different radii. Since the surface area scales as (r^2), the gravitational force must vary as 1/(r^2). In this way, the product of the inward gravitational force and the surface area of a sphere at radius r is independent of r. (At least this is true in the classical picture of a Euclidean universe.)

An easy way to explore the consequences of the inverse square dependence of gravity is to examine what happens when gravity behaves differently. For example, consider the orbit of a planet when the force of gravity is given by the equation
	F= 4(π^2)m / (r^β)
where β is not precisely equal to 2. (In this equation, time is measured in years and radius in AU, the average distance between the Sun and the Earth.)

Explore what happens to the planet's orbit as β changes. When β = 3, planets are eventually ejected from the solar system. For β = 2.1, elliptic orbits precess. (What about circular ones?) Is there a simple relationship between orbit precession and the value of β? 

[For the mathematically inclined, the precession rate can be estimated easily for nearly circular orbits. First, express the equations of motion into polar coordinates, (r, θ). The radial motion depends upon the inward gravitational force and the outward centrifugal force. For circular orbits, these two "forces" balance exactly and (∂_t)r ~ 0. For elliptical orbits,  (∂_t)r = 0 only at the minimum (apogee) and maximum (perigee) distance from the sun. Then, using the conservation of angular momentum (since the force is always central), the rate of precession can be shown to be equal to difference between the angular frequency at perigee and the frequency of oscillation of the radius about the average (circular) radius of the orbit. For nearly circular orbits, the precession rate is (1- (Sqrt[3 - β])) times the orbit rate. This requires an understanding of how to linearize the equations of motion for small deviations from a purely circular orbit.]

The Precession of Mercury

Mercury (the inner most planet) and Pluto (the outermost) have orbits which deviate most from circular. Since Mercury is close to the Sun, the variations to the usual gravitational force caused by Einstein's general theory of relativity cause the elliptical orbit of Mercury to precess. The precession of Pluto's elliptic orbit is imperceptible since Pluto is much further from the Sun than Mercury.

In this problem, you are to investigate the precession of Mercury due to the first (and most significant) correction to Newton's law of gravity due to Einstein. Let the force of gravity be
	F= 4(π^2)m (1 / (r^2) + α/(r^4))
where α ≈ 1.1 × (10^(-8)) when radius is measured in units of AU and time in units of years. The correction, α/(r^4), is very, very small. Mercury's orbit precesses once every 230,000 years! Therefore, in this problem you will need to examine larger values of α. The precessional rate will depend on the value of α. When α is large, the orbit precession rate will be larger than when α is small. (See note below...)

In order to calculate Mercury's orbit, you will have to initialize the orbit with realistic values for its position and velocity. Mercury's semimajor radius (the longer radius of its ellipse, a) is about 0.39 AU, and the eccentricity is e = 0.206. (Remember, the eccentricity relates the semiminor radius, b,  to the semimajor radius, a, by b = a (Sqrt[1 - (e^2)]). The distance from the Sun to Mercury is (1 + e) a = 0.47 AU.)  The velocity of Mercury at its greatest distance from the Sun is about 8.2 AU/year. This value can actually be determined by consideration of conservation of total energy (kinetic plus potential) and conservation of angular momentum. 

Once you've initialized Mercury's orbit, examine its orbit for various values of α. Carefully consider how to let Mathematica determine the changing angle of the elliptic axes of the orbit over time. See if you can develop a relationship between the value of α and the rate of precession of the ellipse. This relation is linear for small values of α. By calculating the precession rate for values of α between 0.001 and 0.01, you can perform a linear extrapolation to α ≈ 1.1 × (10^(-8))and compare your estimate of Mercury's precession to observations.

The Precession of Moon

Because the Moon is attracted to both the Earth and the Sun, lunar motion is fundamentally a three-body problem. The mathematics of three-body dynamics is complex, and the computation of the precession of the Moon's elliptical orbit was one of the most challenging of 18th-century physics. Newton struggled to explain the lunar orbit, but failed. He calculated a precession period of 18 years, which is twice as long as observation. 

Using computational methods, the rate of precession of the lunar orbit can be understood. The key to the correct solution is to allow both the Earth and the Moon to execute orbital motion about each other as the combined Earth and Moon system orbits about the Sun. The Sun can be assummed immobile, but the Earth and Moon are fully dynamic. You can see Newton's error by neglecting the small displacement of the Earth due to the Moon's orbit.

The relative parameters of the the three body system are: Earth: (mass = 5.97 × (10^24) kg, eccentricity = 0.0167), Moon: (mass = 7.348 × (10^22) kg, eccentricity = 0.0549). There is a 6.68 degree tilt of the lunar orbital plane to that of the Earth-Sun orbital plane. I recommend that you ignore this relatively small inclination.

Other Ideas for your Project

You are encouraged to investigate any research topic related to the solution of coupled nonlinear ordinary differential equations. Please contact Prof. Mauel (mauel@columbia.edu) prior to starting your research project