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 nonlinear 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 its also interesting to consider the powerseries solution to this equation. Let u be represented by the infinite series 

_{} 

where c_{0}, c_{1}, 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 u^{2} was absent, the convolutions in the above equations would also be absent, and we would have simply c_{n} = c_{n}_{2 }/ (n(n1)). This is the Newtonian case, and with the initial conditions c_{0} = u_{0} and c_{1} = 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 = h^{2}/m is the semilatus rectum, and r_{0} 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 e^{2} = a^{2} b^{2} , p = b^{2}/a, and r_{0} = 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 u^{2} in equation (1) we must include the convolution terms in the series coefficients. We again find that, with the initial condition c_{1} = 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 semimajor axis is a = (5.79)10^{7} km, and the normalized mass of the gravitating body (the sun) is m = 1.475 km. From these we can compute the perihelion distance r_{0} = (1e)a = (4.59)10^{7} km, and we also have 

_{} 

Setting c_{0} = 1/r_{0}, we can evaluate the power series coefficients of u(q) to give 

_{} 

and so on. Based on the first 60 terms (which is 30 nonzero terms, since all the oddordered 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 its 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 

_{} 

where 
_{} 

The first order of approximation beyond the Newtonian solution is given by the terms in mc_{0} and (m/h)^{2}. The coefficients of the former are 6, 30, 126, 510, , which can be expressed as 2(4^{n} 1) for n = 1, 2, 3, , and the coefficients of the latter are 18, 108, 486, , which can be expressed as 8(4^{n} 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) = c_{0}., since the contributions of the other terms cancel out at this special value. Now, since c_{0} and m/h^{2} are many orders of magnitude smaller than 1 (for typical astronomical situations), the constant term is closely approximated by simply m/h^{2}. 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)^{2}q 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 closedform solution to equation (1). Notice that (1) can be written in the form 

_{} 

In the special case when h^{2} = 12m^{2}. the right hand side vanishes, and we have the equation 

_{} 

This has the exact solution y = 2/(mq^{2}), 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 lefthand figure below. 


Assuming it continues through the singularity as q passes through zero, the orbit then spirals back up and reconverges on a circular orbit of radius 6m, as illustrated in the righthand 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 e^{lx} 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 xe^{lx}. If there is another repetition of the root it is represented by a term proportional to x^{2}e^{lx}, and so on. As a means of approximating the solution of the nonlinear equation (6), many authors introduce a trial solution of the form c_{0} + c_{1}cos(f) + c_{2}fsin(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 (c_{2}/c_{1})f is many orders of magnitude smaller than 1, we can indeed use the smallangle approximations sin(x) ~ x and cos(x) ~ 1 to write the solution as 

_{} 

where weve used the trigonometric identity cos(x)cos(y) + sin(x)sin(y) = cos(xy). 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 onetime 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 secondorder power series expansion of the solution to equation (6), which explains why it is a useful curvefitting 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 smallangle approximation to estimate the frequency shift. 
