Faking relativity

As you know, the orbit of Mercury about the Sun cannot be explained entirely by Newtonian physics. In an ideal Newtonian system, the orbit of a small body around a larger one is a stationary ellipse. In the case of Mercury, it has been observed that the orbit has a noticeable precession: the orbit itself rotates.

Some of this can be explained by applying Newtonian mechanics to the intricacies of the solar system — such as the effects of other planets and the oblateness of the sun. But there was always an extra component to the precession. One of the major early confirmations of Einstein’s theory of relativity is his accurate prediction of this anomaly based on Mercury’s proximity to the immense solar mass.

As noted before, JPL Horizons claims to simulate the solar system using “all known physics”. Can I achieve the same with my simplistic galaxy program? The challenge is much reduced if you aim instead to implement my knowledge of physics!  Fidelity to the latter means taking into account the following:

1. I know about the inverse-square law of gravity.
2. I have access to the position and velocity vector data from JPL Horizons.
3. I heard once of the precession of Mercury’s orbit and that it had something to do with relativity.
4. I know how to use Wikipedia.

The problem with 4. is that there is a wealth of knowledge on Wikipedia, and the more knowledge I have the less feasible the project becomes.  To ameliorate this I’ll add the disclaimer that a pretty picture (or animation) is just as much a goal as an accurate simulation is.  It is realistic to explore physics as a dilettante, but I have no delusions of contributing anything professional in the area.  I acknowledge that I am picking and choosing my physics and apologise in advance to any serious physicists who end up reading this. ;-)

A clock on top of a mountain ticks more quickly than one at sea level.  The lower clock is deeper in the gravity well and is subject to time dilation.  The formula describing the effect is k = sqrt(1 – 2GM/rc^2); 0 <= k <= 1 is a dilation factor for a distance from the Sun of r.  The time dilation effect between two altitudes r_1 and r_2 is the ratio k_1/k_2.  This equation is easy to implement, so I picked it as the starting point for relativistic correction in the simulation.  I added a hack to the apply_force function.  It checks for the Mercury update and calculates k based on Mercury’s distance from the Sun.  It’s particularly messy for two reasons: 1. it’s a special case for one planet only, and 2. it’s implemented in entirely the wrong part of the code.

if (i == 1)  /* First planet is being processed. */
{
    STAR *mercury = s;
    STAR *sun = g->stars[0];
    #define GRAVITY 6.67428E-11
    #define SPEED_OF_LIGHT 299792458.0
    double d2 = get_distance2(mercury, sun);
    double schwarzschild_radius = 2 * GRAVITY * sun->mass / (SPEED_OF_LIGHT * SPEED_OF_LIGHT);
    double k = sqrt(1.0 - schwarzschild_radius / sqrt(d2));
    ts = ts * k;  /* time step to use when updating velocity and position */
}

I’ve separated out the Schwarzschild radius of the Sun r_0 = 2GM/c^2 into its own variable.  (This is constant per body and could be cached.)  Note that at an altitude lower than the r_0, k is imaginary — this appears to correspond to the case of being inside a black hole!  In fact that’s where this formula is too simplistic.  The Schwarzschild radius of the Sun is approximately 3 km.  But this radius forms an event horizon only if all the mass of M is within that radius.  All our physics in this project assumes point particles, so the equation is adequate, except that we cannot expect the simulation to behave like reality when bodies are close enough for this oversimplification to matter.

The next question: how do we apply k to the position and/or velocity of Mercury?  First I tried applying it to both:

/* N.B. vector_add_scaled(v1, v2, f) implements v1 = v1 + f*v2. */
vector_add_scaled(s->vel, forces[i], ts / s->mass);
vector_add_scaled(s->pos, s->vel, ts);

If I had thought about it, I would have realised that exactly the same ellipse is being described, it’s just being drawn more slowly.  :-/ A different time step is being used for that planet, but that is all.

The next crazy attempt is to apply k to just one of the state updates. From an observer’s point of view, everything happening on Mercury happens more slowly.  This includes the application of gravitational force on it.  So: apply the dilation effect only to the velocity update, and assume that the position update can be done without dilation.  This does not result in visible precession — except when we add a small fudge factor:

#define FUDGE_FACTOR 1000000
double k = sqrt(1.0 - FUDGE_FACTOR * schwarzschild_radius / sqrt(d2));
adj_ts = ts * k;
vector_add_scaled(s->vel, forces[i], adj_ts / s->mass);
vector_add_scaled(s->pos, s->vel, ts);

This seems to be what we want:

Click to show animation

We parted ways with accuracy quite early in this relativistic attempt.  A reconciliation (which is not necessarily a priority!) would begin by removing that fudge factor and attempting to adjust the precession to a similar magnitude to that observed in real life — most of which is not relativistic at all:

Sources of the precession of perihelion for Mercury

Amount (arcsec/Julian century)	Cause
5028.83 ±0.04	Coordinate (due to the precession of the equinoxes)
530	Gravitational tugs of the other planets
0.0254	Oblateness of the Sun (quadrupole moment)
42.98 ±0.04	General relativity
5603.24	Total
5599.7	Observed
-3.54 (-0.0632%)	Discrepancy

(From Wikipedia.)  I’m not sure how much the other factors from that table would manifest in the simulation, but relativity accounts for 43 arc seconds per (Earth) century, which amounts to one precession cycle every 3 million years or so.  The video shows one cycle in about 2.5 years. If we remove the fudge factor, and compare the precession amounts with and without the relativistic adjustment, we will get a better idea of how far from reality we are. But that is work for a future post…