Analytical Computation of the Perihelion Precession in General Relativity via the Homotopy Perturbation Method

We propose a new approach in studying the planetary orbits and the perihelion precession in General Relativity by means of the Homotopy Perturbation Method (HPM).For this purpose, we give a brief review of the nonlinear geodesic equations in the spherical symmetry spacetime which are to be studied in our work. On the basis of the main idea of HPM, we construct the appropriate homotopy what leads to the problem of solving the set of linear equations. First of all, we consider the simple example of the Schwarzschild metric for which the approximate geodesics solutions are known, in order to compare the HPM solution for orbits with those obtained earlier. Moreover, we obtain an approximate HPM solution for the Reissner-Nordstorm spacetime of a charged star.


Introduction
The advance of the perihelion in the orbit of Mercury is a relativistic effect [1]. Together with the observation of the deflection of light, this result offers unique possibilities for testing General Relativity (GR) and exploring the limits of alternative theories of gravitation.
The calculations of the perihelion precession in General Relativity and some modified theories of gravity, were recently considered in [2]- [6]. For example, a way to obtain information about higher dimensions from observations by studying a brane based spherically symmetric solution is considered for the classic tests of General Relativity in [7]. The analytical computation of the Mercury perihelion precession in the frame of relativistic gravitational law and comparison with general relativity is presented in [8].
In [9], the perturbations determined by a generic alternative theory of gravity to the GR solution describing the gravitational field around a central mass are worked out. In [10], the authors used recent observations from solar system orbital motions in order to constrain f (T ) gravity. In particular, the spherical solutions of the theory are used to describe the Sun's gravitational field and advances of planetary perihelia in order to obtain upper bounds on the allowed f (T ) corrections.
It is well known that the geodesics equations in RG are nonlinear, and therefore cannot in general be solved exactly. For instance, the geodesic equations resulting from the Schwarzschild gravitational metric element are solved exactly by the Weierstraß Jacobi modular form [11]. Mostly, the perihelion precession of planetary orbits based on Einstein's equations had been calculated in different approximations for a general spherically symmetric line element.
The idea of the Homotopy Perturbation Method which is a semi-analytical method was first proposed by Dr. Ji-Huan He [12]- [14] for solving differential and integral equations. Later, the method is applied to solve the non-linear and non-homogeneous partial differential equations [15].
The HPM has a significant advantage providing an analytical approximate solution to a wide range of nonlinear problems of the fundamental and applied sciences [16], [17]. Sometimes, this method allows to find even an exact solution with the help of a few iteration [18]. This method and a wide spectrum of its application have been extensively developed for several years by numerous authors (see [19] and references therein).
Recently there were studies in which this method was used for analytical calculations in the field of cosmology and astrophysics (see, e.g. [20], [21]). Our aim is to give one more application of the method to the problem of planetary motion in the spherically symmetric gravitational field in General Relativity.

Geodesics in the Spherical Symmetry Spacetime
As well known [1,3], the line element of the 4dimensional general spherically symmetric stationary spacetime can be written as and the perihelion precession is usually treated as the time-like geodesic in spacetime. Let us consider the geodesics γ(τ ) in the above spherically symmetric spacetime. We set the geodesic γ(τ ) expressed in the spherical coordinates x µ = (t, r, θ, ϕ) as x µ (τ ), which are satisfied The geodesic γ(τ ) can be obtained by solving the above equation. However, taking into account the symmetry of spacetime (1), one could use the following simple way to obtain the geodesic γ(τ ). First, we can find that one component of the geodesic γ(τ ) can always be chosen as θ(τ ) = π/2, which means that the geodesic can always be chosen to lay in the equatorial plane of the spherically symmetric spacetime. Thus, t = t(τ ), r = r(τ ), θ = π/2, ϕ = ϕ(τ ). Let us denote the tangent vector of geodesic γ(τ ) as U µ ≡ dx µ /dτ . For the time-like geodesic, we chose τ to be the proper time . Hence, from (1) we can obtain where we have used θ = π/2. Second, it could be noted that ξ a = (∂/∂t) a and ψ a = (∂/∂ϕ) a are two Killing vectors in the spherically symmetric spacetime (1). Therefore, there are two conserved quantities along the geodesic γ(τ ), the total energy and the angular momentum per unit mass After inserting (3) and (4) into (2), one could obtain This equation contains only one function r(τ ), and it could be solved in principle. Then, after inserting the solved r(τ ) into (3) and (4), the rest components t(τ ) and ϕ(τ ) of geodesic could be finally obtained. However, it should be pointed that perihelion precession is usually related to the orbit of geodesic, i.e. r(ϕ). Therefore, it is convenient to rewrite the equation (5) with the help of (4) as It has been found that the coordinate u ≡ 1/r is more convenient than r to derive the perihelion precession. Thus, the main equation investigated in our paper could be simply obtained from equation (6) by converting r into u du dϕ Finally, differentiating the equation (7) with respect to ϕ, we get the second-order geodesic equation in the following form

Homotopy Perturbation Method
To illustrate the basic ideas of HPM [12] for solving nonlinear differential equations, let us consider the following nonlinear differential equation: with the boundary conditions B(u, ∂u/∂n) = 0; r ∈ Γ , where A is a general differential operator, B is a boundary operator, g(r) is a known analytic function, Γ is the boundary of the domain Ω . Suppose the operator A can be divided into two parts: M and N . Therefore, (9) can be rewritten as follows: The homotopy v(r, p) : Ω × [0, 1] → IR constructed as follows [13] H where r ∈ Ω and p ∈ [0, 1] is an imbedding parameter, and y 0 is an initial approximation of (9). Hence, one can see that and changing the variation of p from 0 to 1 is the same as changing H(v, p) from M (v) − M (y 0 ) to A(v) − g(r), which are called homotopic. In topology, this is called deformation. Due to the fact that 0 ≤ p ≤ 1 can be considered as a small parameter, by applying the perturbation procedure, one can assume that the solution of (11) can be expressed as a series in p, as follows: When we put p → 1, then equation (11) corresponds to (10), and (13) becomes the approximate solution of (10), that is u( This series is convergent for most cases. However, the convergent rate depends upon the nonlinear operator A(v). Sometimes, even the first approximation is sufficient to obtain the exact solution [12]. As it is emphasized in [13] and [17], the second derivative of N (v) with respect to v must be small, because the parameter p may be relatively large, i.e. p → 1, and the norm of L −1 ∂N/∂v must be smaller than one, in order that the series converges.

Orbits and Perihelion Precession via HPM
First, we consider the simplest case of the metric (1), namely, the Schwarzschild metric describing the gravitational field of an uncharged non-rotating star. For the Schwarzschild solution, the two functions are f (r) = h(r) = 1 − 2M/r, or where M is a mass of the star. Therefore, Eq.(8) for the time-like geodesic can be Compared with the case in Newton's gravity the term 3u 2 M comes from the correction of general relativity. Note that, the analytical solution of (15) is absent like the Schwarzschild case. However, there is an approximation solution of (15) is the analytical elliptical solution which has already been found in Newton's gravity, and e is the orbital eccentricity which has been considered as a small constant. Now on, we consider the HPM of solving equation (15). For this end, we suppose the following homotopy where the prime denotes derivative with respect to ϕ, and assume that the solution of (15) can be expressed as a series in p by u(ϕ) = u 0 (ϕ) + p u 1 (ϕ) + p 2 u 2 (ϕ) + ... .
According to (17), the initial conditions for u 0 (0) and u i (0) can be chosen as follows where i ≥ 1. The substitution of (19) into equation (18) yields It is noteworthy that we obtain the set of linear equations. Their solutions with the initial conditions (20), (21) can be readily found as u 0 (ϕ) = M L 2 (1 + e cos ϕ), where we have deliberately limited our calculation by the minimum degree of approximation. All subsequent approximations can also be obtained easily.
In accordance with the HPM, it follows from (19) and (25) that the solution of equation (15)  Comparing our result (26) with the approximate formula (17), one can conclude that these solutions are different in the term −(M 3 /L 4 )(3 + e 2 ) cos ϕ, and therefore they do not give the same value of the perihelion shift. In order to demonstrate the difference in the numerical and approximate solutions of the equation (15) more clearly, we take some hypothetical parameters of the system, providing significant perihelion shift compared to that of Mercury. So we set M = 0.04, L = 0.2 in conventional units, and e = 0.8. Typically, the perihelion shift is found from the appropriate approximate solutions, neglecting e 2 compared to e. Thus, for small eccentricity, the equation (17) can be rewritten as where ǫ = 3M 2 /L 2 . For the perihelion of orbit, it satisfies cos(ϕ − ǫϕ) = 1. Therefore, the precession angle of perihelion equals to ∆ϕ = 2πǫ = 6πM 2 /L 2 . A more precise formula for the shift can be obtained from the condition u ′ (ϕ) = 0 in the perihelion, and ϕ = 2π + ∆ϕ, assuming that ∆ϕ ≪ 1. As applied to the approximate solution (17), this method yields At the same time, the use of this approach to our solution (26) gives the following result Note that even for a small value of eccentricity, we have ∆ϕ HP M = ∆ϕ, but both of them lead to the same expression for ∆ϕ ≈ 6πM 2 /L 2 when L 2 /M 2 ≪ 1.
Finally, let us obtain the HPM orbits and shift in the Reissner-Nordstorm spacetime of a charged star. In this case, we have [1] where Q is the charge. According to (8) and (28), the main equation (15) is replaced by the following one Assuming that the unperturbed equation should have solution (8), consider the following homotopy where p ∈ [0, 1]. Substituting (19) into equation (30), we get where the simplest approximation is taken. The set of linear equations (31), (32) with the initial conditions (20), (21) can be easily solved, giving for the approximate solution u(ϕ) = u 0 (ϕ) + u 1 (ϕ).
With the help of our solution (33), one can easily obtain the value of shift angle. We do not provide it here due to a cumbersome nature of its expression.

Conclusions
Thus, in this work we have considered a simple analytical computation of the perihelion precession in General Relativity with the help of the Homotopy Perturbation Method. First, we have studied the example of geodesic motion in the Schwarzschild metric, in order to approbate HPM in the problem of planetary motion in GR, and present the main steps in solving by this method.
The comparison of our solution (26) with the approximate solution (17), obtained earlier by the perturbation method, demonstrates the better degree of accuracy. This result has been shown in Fig 1. A direct result of the differences in solutions consists in that the perihelion precession (27), derived from the HPM solution (26), yields the better accuracy compared to the previously known result. It is important that formula (27) for ∆ϕ HP M differs from the standard ∆ϕ, even for the small values of eccentricity. Only in the case of small parameter L 2 /M 2 ≪ 1, the HPM approximation could give the commonly used value of precession 6πM 2 /L 2 . It is worthy to note, that all these results were obtained by the single iteration. This may give us the hope that the next iteration could provide us the result with much greater accuracy and a minimum size of computations.
Moreover, in our work we have obtained HPM solution for the Reissner-Nordstorm spacetime of a charged star. For this case, the additional terms in the HPM approximation were obtained as compared to the similar solution in the standard approximation (see, e.g. [3]).
According to the results of this work and our work published earlier [20], we could make the following remarks. Foremost, the undoubted advantage of this method consists of that there is no need to establish a small parameter for solving a problem in some approximation, because such a small parameter sometimes could destroy the main feature of the exact solution. On the other hand, the approximate solution, and the rate of its convergence in this method greatly depend on the construction of homotopies. Nevertheless, we argue that HPM is able to provide a rather high degree of accuracy in the problems of astrophysics and cosmology associated with nonlinearity of their main equations.