Series Solution for Relativistic Orbits


Solving the field equations of general relativity for the orbit of a small test particle in Schwarzschild spacetime leads to the non-linear differential equation



where m is a measure of the mass of the gravitating body, h is a measure of the angular momentum, and u = 1/r. The parameters m and h are both normalized to have units of length. The second derivative of u in the above equation is with respect to the angular position q in the orbit. The second term on the right side is what distinguishes this equation from the Newtonian equation of motion. An approximate solution of (1) for planetary orbits is given in Anomalous Precesion, but it’s also interesting to consider the power-series solution to this equation. Let u be represented by the infinite series



where c0, c1, … are real constants. Inserting into (1) and collecting terms by powers of q, and setting each coefficient to zero, we get the following conditions on the coefficients:



If the term of (1) involving u2 was absent, the convolutions in the above equations would also be absent, and we would have simply cn = -cn-2 / (n(n-1)). This is the Newtonian case, and with the initial conditions c0 = u0 and c1 = 0 we have



Thus the bound Newtonian orbits satisfy the equation



which is the equation of a conic section with the gravitating body at one focus. The ratio p = h2/m is the semi-latus rectum, and r0 is the perihelion radius. Letting a and b denote the major and minor radii and e the linear eccentricity of the ellipse, we have the elementary relations e2 = a2 – b2 , p = b2/a, and r0 = a – e, so the quantity in parentheses in the above equation can be written as



Thus the quantity in parentheses is the numerical eccentricity e = e/a. The orbit is an ellipse, a hyperbola, or a parabola respectively depending on whether |e| is less than, greater than, or equal to 1.


To account for the relativistic correction term involving u2 in equation (1) we must include the convolution terms in the series coefficients. We again find that, with the initial condition c1 = 0, all the subsequent coefficients of odd powers of q are also zero, so u(q) is an even function (meaning it is the same in both directions of q).


For the planet Mercury the orbital eccentricity is e = 0.026, the semi-major axis is a = (5.79)107 km, and the normalized mass of the gravitating body (the sun) is m = 1.475 km. From these we can compute the perihelion distance r0 = (1-e)a = (4.59)107 km, and we also have



Setting c0 = 1/r0, we can evaluate the power series coefficients of u(q) to give



and so on.  Based on the first 60 terms (which is 30 non-zero terms, since all the odd-ordered terms are zero), the function u(q) and its derivatives are plotted below.



The first derivative is zero at each perihelion and apihelion, and the first apihelion occurs at q = 3.14159290423…, which exceeds p by (2.506)10-7 radians. Since u(q) is an even function, we know the increment is the same for the apihelion near –p, so the total precession per revolution is twice this amount, which is (5.012)10-7 radians per revolution, equal to 0.1034 arc seconds per revolution. Mercury completes 414 revolutions per century, so it’s relativistic precession is 42.9 arc seconds per centuy, in close agreement with the observational value.


Algebraically the first several coefficients of the power series of u(q) are





The first order of approximation beyond the Newtonian solution is given by the terms in mc0 and (m/h)2. The coefficients of the former are 6, 30, 126, 510, …, which can be expressed as 2(4n – 1) for n = 1, 2, 3, …, and the coefficients of the latter are 18, 108, 486, …, which can be expressed as 8(4n – 1) – 6n for n = 1, 2, 3, …  Making these substitutions and consolidating the resulting series into trigonometric functions, we get the result



Naturally we still have u(0) = c0., since the contributions of the other terms cancel out at this special value. Now, since c0 and m/h2 are many orders of magnitude smaller than 1 (for typical astronomical situations), the constant term is closely approximated by simply m/h2. Also, the term involving cos(2q) represents just a very slight periodic oscillation affecting the shape of the orbit, with no secular effect. Thus to determine the secular effects we can examine the close approximation



The quantity 3(m/h)2q is very small for the first few revolutions, so for this part of the curve the quantity in the square brackets can be closely approximated by



Consequently the orbit precesses by 6p(m/h)2 per revolution. The term involving q sin(q) appears naturally from the conditions on the power series coefficients, but it is often introduced simply as a trial solution, and/or justified as a resonance term. Of course, for sufficiently large angles q, the approximation for u(q) increases without limit, but m/h is so small that the above expressions are valid over many revolutions.


Incidentally, for some special cases we can find an exact closed-form solution to equation (1). Notice that (1) can be written in the form



In the special case when h2 = 12m2. the right hand side vanishes, and we have the equation



This has the exact solution y = 2/(mq2), so we have




Hence the orbit spends infinitely much time arbitrarily close to a circular orbit at radius 6m for large negative values of q, but then as q approaches zero the orbit spirals down to the center of the field, as shown in the left-hand figure below.



Assuming it continues through the singularity as q passes through zero, the orbit then spirals back up and re-converges on a circular orbit of radius 6m, as illustrated in the right-hand figure above.


Incidentally, the relativistic contribution to a planet's orbital precession rate is often described as a "resonance" effect. Recall that the general solution of an ordinary linear differential equation contains a term proportional to elx for each root l of the characteristic polynomial, and a resonance occurs when the characteristic polynomial has a repeated root, in which case the solution has a term proportional to xelx. If there is another repetition of the root it is represented by a term proportional to x2elx, and so on.  As a means of approximating the solution of the non-linear equation (6), many authors introduce a trial solution of the form c0 + c1cos(f) + c2fsin(f), suggesting that the last term is to be regarded as a resonance, whose effect grows cumulatively over time because the factor f is not periodic, and therefore eventually has observable effects over a large number of orbits, such as the 414 revolutions of Mercury per century.  Now, provided (c2/c1)f is many orders of magnitude smaller than 1, we can indeed use the small-angle approximations sin(x) ~ x and cos(x) ~ 1 to write the solution as



where we’ve used the trigonometric identity cos(x)cos(y) + sin(x)sin(y) = cos(x-y). This yields the correct result, but the interpretation of it as a resonance effect is misleading, because the predominant cumulative effect of a resonant term proportional to fsin(f) in the long run is not a precession of the ellipse, but rather an increase in the magnitude of the radial excursions of a component that is at an angle of p/2 relative to the original major axis. It just so happens that on the initial cycles this effect causes the overall perihelion to precess slightly, simply because the phase of the sine component is beginning to assert itself over the phase of the cosine component.  In other words, the apparent "precession" resulting from the fsin(f) term on the initial cycles is really just a one-time phase shift corresponding to a secular increase in the radial amplitude, and does not actually represent a change in the frequency of the solution. It can be shown that a term involving fsin(f) appears in the second-order power series expansion of the solution to equation (6), which explains why it is a useful curve-fitting function for small f, but it does not represent a true resonance effect, as shown by the fact that the ultimate cumulative effect of this term is discarded when we apply the small-angle approximation to estimate the frequency shift.


Return to MathPages Main Menu