An elementary approach to simulating the perihelion of Mercury

The relativistic correction to the precession of the perihelion of Mercury provided key evidence for the accuracy of general relativity as a theory of gravity. This example still has a large amount of potential to introduce students to the power of numerical simulations in theoretical physics, but existing approaches may be too detailed for many students and involve them beginning to learn a programming language at the same time. In this article, we take a simpler approach which uses as little coding as possible. The equation for the orbit of a planet is solved with and without relativistic corrections. It is shown that there is precession of the perihelion in the relativistic case, whereas in the Newtonian case, the orbit does not rotate about the origin. Quantitative information is extracted on the precession of the perihelion of Mercury and shown to match with observations.


Introduction
It is well-known that general relativity introduces some subtle corrections to Newtonian mechanics.One of the most famous of these is the precession of the perihelion of Mercury [1,2].The observed rate of precession of Mercury's perihelion disagrees with Newtonian mechanics, even after all other factors have been accounted for (the perihelion of an orbit being the point where the body comes closest to the Sun) [3].Most of the observed secular precession of the perihelion of Mercury can be explained via precession of the equinoxes and perturbations of the orbit by other planets, but after these effects are dealt with, there is still an observed precession of approximately 43 arcseconds per century.One takes the Kepler problem for a bound orbit perturbed by a central force and includes a perturbation Hamiltonian due to the Schwarzchild solution for weak gravitational fields, one obtains a value for the average rate of precession which is 42.98 arcseconds per century.Observation currently tells us that the precession due to relativistic effects is 43.1 ± 0.5 arcseconds per century, an insignificant deviation from theory [4].The perihelion of Mercury is still a topic of active research and new relativistic corrections are currently under consideration by some authors [5][6][7].A good history of work on the advance of the perihelion of Mercury may be found at [8].For more background, readers can also find an interesting theoretical analysis of these historical attempts to explain the perihelion in terms of Lagranian mechanics at [9].
In this article, we discuss the pedagogical utility of the perihelion of Mercury as an introduction to general relativity and computer solutions of differential equations.This has been discussed in the literature before by Koerber et al, who carried out a very thorough study of the perihelion [10].However, although the work of Koerber et al claims to be accessible to advanced high school students, this seems questionable considering the amount of mathematics that they use, as well as a large amount of sophisticated Python code which might be confusing even to undergraduates.We should emphasise that students in certain European countries are required to learn steps for programs to be implemented in Python, in which case [10] would fit well into the curriculum.For this reason, we view our approach as complementary to that of Koerber et al, with a set of aims which are different and simpler.One could also view our work as a primer for the much more detailed approach of Koerber et al, which might seem intimidating to a student considering this topic for the first time.In their work, it was stated that only advanced students would be able to extract any quantitative conclusions from their approach using input parameters from general relativity, whereas in our approach we believe that all students should be able to follow both the mathematical argument and the code to quantify the precession of the orbit.Koerber et al use the programming language Python to study the motion and a separate module VPython to visualise trajectories.In our opinion, this is not ideal for high school students and undergraduates as it may lead to them becoming wrapped up in the complexity of starting to learn a particular progamming language and losing sight of the physics which is involved.Again, we should emphasise that in certain countries it is a required part of the curriculum to learn Python, in which case this criticism does not hold.
We instead suggest using Maple which is much easier to get started with and does not require a clean programming style or any knowledge of programming [11].Mathematica could also be used if preferred.In section 2, we write down the orbit equation which one obtains by considering the classical Kepler two-body problem and solve it numerically to show that there is no precession in this case, with an eternal ellipse that does not rotate about the origin.In section 3, we write down the corresponding relativistic correction and solve this new equation to show that there is a precession of the perihelion when general relativity is brought into play.In section 4, we show that our code can reproduce the value for the precession of the perihelion of Mercury when suitable input parameters are introduced.An advantage of our approach is that this does not involve any additional complexity, whereas extracting quantitative information from the code of Koerber et al was considered to be for more advanced students [10].We should note that the simulation we discuss in this article is most suitable for general undergraduate students, whereas the topic of general relativity is often designated for advanced undergraduates or even beginner graduate students.However, even if this might be the case, it is still desirable to at least start introducing the concept of general relativity as soon as possible in the undergraduate curriculum so that students are aware of it, in which case our simulation seems ideal.In section 5, we finish by discussing the pedagogical utility of our work in comparison to similar studies in the literature.

Solving the orbit equation
The Kepler problem considers two bodies interact by a central force F whose strength varies depending on the square of the distance.The classic example is that of a light body like Mercury orbiting a heavy body like the Sun.In vector form, the central force vector can be written as k r F r.

= Note
that this force can be attractive or repulsive depending on the sign of the constant k.
It can be shown by solving the Kepler problem that the trajectory of the light body around the heavy body restricts to a two-dimensional plane and is given by a differential equation r r ⎞ ⎠ where G is the gravitational constant, M is the mass of the heavy body, and h = Mr 2 ω.In orbital analysis, one usually introduces a new variable u = 1/r and re-writes this as As we are keeping things elementary, the derivation is too involved and uses too much trickery to be reproduced here, but it can be found in classical mechanics textbooks [4].If we take this equation as a black box which gives the orbit of a body, we can now start with the Maple code which will solve this equation numerically.To best illustrate the effect of the precession, we begin by setting all the constants equal to unity: restart; with(plots): M:= 1; G:= 1; h:= 1; We next need to specify the differential equation and two initial conditions since it is a second-order ODE: we will use r r 0 2 3, 0 0 = ¢ = ( ) ( ) .

Solving the orbit equation with the relativistic correction
In general relativity, the orbit equation receives an additional correction term and one now has ⎞ ⎠ where c is the speed of light in a vacuum.Similarly to the previous section, the derivation of this equation is too advanced to reproduce here but can be found in relativity textbooks [3].Students with knowledge of the Schwarzschild metric should be able to follow the derivation.The basic idea is that in general relativity spacetime has curvature.This curvature corresponds to the gravitational field and the heavy body contains almost all of the mass in the system, so the light body follows a trajectory which is slightly deflected by the essentially static heavy body.In terms of forces, this is observed as a new attractive force.This force means that the light body is slightly more attracted to the heavy body than it would be in Newtonian gravity, causing the orbit to precess very slightly in the direction of rotation.
As before, we choose values for the constants which will help to make the precession more obvious (we choose c = 8).We then solve the ODE numerically and create a plot of the orbit (r(f), f) for 0 f < 10π.The result is shown in figure 1 (right).It can be seen that the elliptical orbit has rotated five times about the origin.

Quantifying the precession for Mercury
In common with [10], we have shown that Maple can be used to visualise the precession of the orbit of a light body due to relativistic corrections.However, we can now go further and make a quantitative connection with the precession of Mercury (in [10] this part of the analysis was considered to be for advanced students only).To begin, we need to find the angular positions f 1 , f 2 , f 3 and f 4 of the first four perihelia in figure 1 (right).We search for the minimum distance of the orbit from the origin for each revolution (ie.whenever r (f) has a minimum).We include the second derivative condition for a local minimum in the definition of the ODE and the initial conditions, specifying that a revolution stops when this condition is met: It is found that the advance of the perihelion δf at each revolution is approximately constant and is given by 2 0.336 25 radians.
This is a substantial precession but since the relativistic correction is inversely proportional to c 2 , in reality, this correction will be negligible unless the orbital velocity is comparable to c. Relativistic effects are mostly negligible when considering orbits of the planets in the Solar System, with the very important exception of Mercury.We plug in the appropriate input parameters and initial conditions for the orbit of Mercury and then calculate the precession as before: In this case, the advance of the perihelion at each revolution is approximately 0.506 × 10 −6 radians.Mercury orbits the Sun every 88 days and there are 36 525 days in 100 Earth years, so dividing through we have 415.056818 2 revolutions.Note that we keep a large number of decimal places in our work as accuracy is very important here.To get the total precession in arcseconds, we multiply by the number of radians of precession in one revolution, convert to degrees, and then multiply by 3600 (since 1 arcsecond is equal to 1/ 3600 of a degree).This gives us 43.319 476 76 arcseconds per terrestrial century.

Conclusions
In this work, we have revisited the idea of using the advance of the perihelion of Mercury as a topic to motivate the use of numerical simulations in theoretical physics.Although the work of Koerber et al is ambitious, we do not think that it achieves its goal of being accessible to high school students (or even most undergraduates) and that the treatment in terms of the coding is unnecessarily complicated.Obviously, there is the important exception of many European countries where the use of Python and writing of computer programs is an obligatory part of the physics curriculum [10].In their work, it was stated that they wished to communicate the beauty of theoretical physics to beginner students with an important historical example, but the amount of detail provided was perhaps overwhelming for students approaching this topic for the first time.For this reason, we were motivated to provide a simpler approach with slightly different aims which can be followed by all students and used as a primer for those who wish to move on to the more involved discussion in Koerber et al [10].In our work, we retain the concept of using numerics to solve a differential equation with a subtle relativistic correction and provide visualisations and quantitative information about the precession of the perihelion, but the coding is now much easier to follow and does not ask students to start learning a programming language in order to carry it out.

Figure 1 .
Figure 1.(Left) Plot of numerical solution of the orbit equation and (right) plot of numerical solution of the orbit equation with the relativistic correction.