* Last updated 17 April 2011 *

- Why include relativity in orbit computations?
- An example of how much it can matter
- By how much does the longitude of perihelion precess due to GR?
- A pretty cheap and easy way to do it
- A better way of doing it
- Another good way of doing it
- Perhaps a better way of doing it?
- Warning
- Some more information about this subject
- But is it really general relativity?

** Why include relativity in orbit computations?: **
Quite a bit of orbital software (ephemeris generators, orbit
determination software, etc.) just deals with Newtonian mechanics:
you've got the sun, some perturbing objects, and that's it, all
working with the simple Newtonian inverse-square law. This is true
of my Find_Orb orbit determination software,
and of all the software I currently have released. The
Lowell Observatory ASTORB
database also neglects relativity; the documentation states
that this matters slightly for (1566) Icarus, and probably for other
objects that approach the sun closely. JPL and NeoDyS/AstDyS include
this effect. The Minor Planet Center usually includes it, if the
object is sufficiently affected by it.

The most famous example of a situation where "relativity matters" for orbit computations is the precession of the orbit of Mercury. From observations, it was determined that the longitude of perihelion (sum of the longitude of the ascending node and the argument of perihelion) was increasing at 574 arcseconds/century. Planetary perturbations can account for about 531 arcseconds/century, and until Einstein came along, the remaining 43 arcseconds/century was something of an embarrassment.

** An example of how much it can matter: ** In April 2005,
Rob Matson asked me if perhaps my Find_Orb
orbit determination software was neglecting general relativity.
He was having problems using it on an object with an orbit with a
small semimajor axis. (Perhaps an "NSO", or "Near-Sun Object".)
We quickly determined that the problem was elsewhere, but also that
Find_Orb indeed did * not * include general relativity. I
quickly fixed this. Rob suggested the object (137924) 2000 BD19 as a good
test case.

The data available for (137924) 2000 BD19 ranged from January 2000 to March 2005, with two precovery observations from February 1997. Click here to see a "pseudo-MPEC" for this object with general relativity ignored. As you will see, it leads to the precovery observations appearing to be in error by about a half arcminute. As a result, you could easily assume they were poorly timed or measured, and reject them.

However, if you look at a similar "pseudo-MPEC" but with general relativity correctly computed, you'll see that the precovery observations suddenly look perfectly reasonable. (Rob then suggested that we might be able to find this object on a plate from the 1950s, but it turns out that the object was just barely off the margin of the plate.)

** By how much does the longitude of perihelion precess due to GR? **
In general, the rate of perihelic precession will be given by

precession = 3.838"/[a^{2.5}* (1-e^{2})] per century = 3.838"/[P * a * (1-e^{2})] per century

where a = semimajor axis in AU, e = eccentricity, P = period in years. I am indebted to Rob Matson for pointing this out in a post on the Minor Planet Mailing List. Alternatively, as described here , the rate can be computed as (I've rearranged the formula a hair):

precession = 6 π GM / [a * c^{2}* (1-e^{2})]

which gives the perihelion shift per orbit for a planet with semimajor axis a and eccentricity e, orbiting a mass M (for Mercury, the mass of the Sun).

** A pretty cheap and easy way to do it: **
I came across this trick in
J. M. A. Danby's * Fundamentals of
Orbit Determination, * published by
Willmann-Bell. It mentions it as part of a "problem for the student",
and I do not completely follow the derivation. However, the actual
use of it is quite straightforward, and (as I will discuss) I've tested
it and found that it works (though you should probably use the
better way , or this yet
another way).

The basic idea is that the standard Newtonian gravitational acceleration,

accel = -GM / r^{2}

is modified slightly, as follows:

accel = -GM / r^{2}- alpha * h^{2}/ r^{4}

where

alpha = 3GM / c^{2}h = |(v x r)|^{2}= square of the angular momentum/unit mass

It would be perfectly OK to stop here, and in fact, the code I wrote
does so. But it's modestly educational to note that, if we set
v_{t} equal to the transverse component of the velocity, then
we can express this "modified acceleration" as

accel = -(GM / r^{2}) (1 + 3(v_{t}/ c)^{2})

in other words, the Newtonian acceleration is increased by a factor caused by the object's velocity.

I tested this by writing a little program to integrate the orbit of a simulated Mercury, printing out the running rate of the longitude of perihelion as it went (that is, the change in lon_per so far divided by the integration time elapsed so far). The program uses the method of Encke, meaning that it solves the two-body Newtonian problem analytically and integrates only the departure from this idealized orbit caused by perturbations. I included the relativistic term as "just another perturbation".

When I ran this program with no perturbing objects and Newtonian gravity, the elements remained exactly constant (as you would expect). Including the relativistic term caused some fluctuations in the rate of precession of the longitude of perihelion, but it eventually settled down to 43"/century. (The extra term causes both secular and periodic effects, which is why I had to integrate over a few hundred orbits to see the long-term effects.)

Including the perturbing objects Venus-Neptune, but no relativity,
got me the 531 arcseconds/century precession that everybody got in the
nineteenth century. Finally, including perturbations * and *
relativity got me the "real world" value of 574 arcseconds/century.

** A better way of doing it: ** I * thought * that the
above handled relativistic effects correctly to first order. I was
wrong. In February 2003, I received an e-mail from Werner Huget.
He had written a similar program and had used it to integrate the
motion of Mercury. He found a slight error in the mean anomaly,
amounting to 8" after a seven-year integration!

Some time ago, I had tested my Integrat software by going to the JPL Horizons ephemeris site and getting osculating orbital elements for (1566) Icarus for the epochs JD 2452400.5 = 6 May 2002 and JD 2380000.5 = 14 Feb 1804 (Horizons wouldn't go back much farther). The idea was that if Integrat could integrate the former elements to the latter elements, I'd have a lot of confidence that the program was functioning properly. I chose (1566) Icarus because its elongated orbit with a perihelion of only .187 AU is strongly affected by relativity.

With this test, I was able to replicate JPL's 1804 elements almost exactly, except for a 15" error in mean anomaly. Thanks to Werner Huget, I now know why I had that error, and have fixed it (the current version of Integrat replicates the 1804 orbit to within about .1 arcsecond).

The revised differential acceleration (amount to be added to the "normal" Newtonian acceleration) is

-GM accel = ---- {(4GM / r - v^{2})r+ 4(v.r)v} r^{3}c^{2}

where **r** = vector from the sun to the planet and **v**
= velocity vector of planet relative to the sun. **v**.**r**
is the vector dot product of the velocity and radius. Note that the
acceleration is no longer purely radial. With a bit of rearrangement,
you can show that for low-eccentricity orbits, the above formula
reduces to the simpler version shown above. But for even slightly
eccentric orbits (as Werner Huget found out with Mercury and I found
out with Icarus), the "simple version" isn't good enough.

A similar expression is given in the
* Explanatory Supplement to
the Astronomical Almanac, * p 281, except that some higher-order
terms are included in that book.

** Another good way of doing it: ** When this subject came up
on the Minor Planet
Mailing List, Aldo Vitagliano made an
interesting comment about a method that has certain advantages
over the above methods. In his
Solex software for integrating planetary motions, he
has an option to compute the revised differential acceleration as

-(GM)^{2}accel = ----- (6 / r - 9 / a)rr^{3}c^{2}

where a = semimajor axis of the orbiting object. The benefit of this
is that if you assume that a varies slowly enough, then you can treat it
as semi-constant and the above expression doesn't include the object's
velocity. That's helpful if you're using any of several second-order
differential equation integration algorithms that assume you don't know
the object's velocity. That is, these schemes can integrate equations
of the form x" = d^{2}x/dt^{2} = f(x, t), but not those
of the form x" = d^{2}x/dt^{2} = f(x, x', t).

Incidentally, it can be useful to know that 1/a = 2/r - v^2/GM. Using this, the above expression can be rearranged to yield

-9(GM)^{2}accel = ----- (2 / r - v^{2}/ GM)rr^{3}c^{2}

The result of this is * not * identical to that from the original
"better way" suggested by Werner Huget. Vitagliano's method results in a
purely radial force; Huget's adds a "sideways" term parallel to the
velocity vector. I think Huget's method is, strictly speaking,
more accurate; the difference is probably a periodic term of extremely
small magnitude.

** Perhaps a better way of doing it? ** Henrik Agerhäll recently
sent me a copy of some work he has done in this area, suggesting

-3(GM) accel = ----- (v.r)vr^{3}c^{2}

He mentions that for actual solar system cases, this produces essentially the same result as the above "better" method; but for an object much closer to the sun (or in some other situation involving very strong gravitational fields), it gives a better result than the above "better" method. There is an explanation here.

This method results in no relativistic acceleration for a circular orbit, and none for a non-moving object. I've tried it in my Find_Orb software, on (137924) 2000 BD19, and found that the residuals changed by a milliarcsecond or two.

** Warning: ** As mentioned above, I do not really follow all
of the reasoning behind the derivation of the first expression, and
don't have even a clue as to the derivation of the latter expressions.
For "practical" use in solar system work, it appears that any of
these expressions will suffice to handle relativistic effects.

However, they *do* yield slightly different results, and
they can't all be "right". I don't know which one (if any) actually
corresponds to the physical universe.

At least at some level, none of them do. You'll see that gravitational waves are ignored in all of them. For the solar system, such an effect is absolutely trivial. I recall hearing once that the Earth-Sun system radiates about 60 watts of gravitational waves.

In one famous case, there is a pair of pulsars that orbit each other so closely, and emit so much energy via gravitational waves, that their orbit is gradually contracting, by an amount that is measurable. (The fact that the increase in orbital period matches the energy loss from gravitational waves is, perhaps, the best evidence to date that gravitational waves really exist.)

If you wanted to investigate such behavior, the above methods would not be enough. As you can see, they conserve energy, so no gravitational waves can be involved.

** Some more information about this subject: **
On 28 Apr 2002, I got the following comments from Ken Shoemaker:

On your web page, you mention that you found a trick for incorporating some relativity, in Danby's book. What was missing was an explanation. I recall a nice discussion of such approximations in "Gravitation", by Misner, Thorne, and Wheeler... [The book] necessarily gets into differential geometry, but presents most topics in at least three different ways: (1) graphically, (2) with coordinates, and (3) more abstractly (no coordinates). I'd also suggest you check out John Baez's web pages, http://math.ucr.edu/home/baez/gr/gr.html.

** But is it really general relativity? ** This may sound a
bit out there in crackpot-land, but
this paper
makes a pretty good case that you don't need general relativity to
explain Mercury's precession or the aberration of starlight near the
Sun; special relativity will do the trick (with some not especially
hard-to-believe assumptions). The author does go so far as to say
that his alternative theory "disproves" GR, which is indeed
a stretch; it's more accurate to say that he has an alternative that
explains Mercury's precession (and the bending of light near the sun)
exactly as well as GR does. (I don't know if his alternative would
explain gravitational waves, though.) The paper makes some predictions
about the internal structure of Mercury, which may be verified (or
not) when MESSENGER goes into orbit around that planet in 2011.