Including relativistic effects in orbit computations

Last updated 17 April 2011

Italian version

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"/[a2.5 * (1-e2)] per century
           = 3.838"/[P * a * (1-e2)] 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 * c2 * (1-e2)] 

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 / r2

is modified slightly, as follows:

accel = -GM / r2 - alpha * h2 / r4


alpha = 3GM / c2
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 vt equal to the transverse component of the velocity, then we can express this "modified acceleration" as

accel = -(GM / r2) (1 + 3(vt / 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

accel =  ---- {(4GM / r - v2) r + 4(v.r)v}

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

accel =  ----- (6 / r - 9 / a) r

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" = d2x/dt2 = f(x, t), but not those of the form x" = d2x/dt2 = 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

accel =  ----- (2 / r - v2 / GM) r

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

accel =  ----- (v.r)v

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, 

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.