Computation of the general relativistic perihelion precession and of light deflection via the Laplace-Adomian Decomposition Method

We study the equations of motion of the massive and massless particles in the Schwarzschild geometry of general relativity by using the Laplace-Adomian Decomposition Method, which proved to be extremely successful in obtaining series solutions to a wide range of strongly nonlinear differential and integral equations. After introducing a general formalism for the derivation of the equations of motion in arbitrary spherically symmetric static geometries, and of the general mathematical formalism of the Laplace-Adomian Decomposition Method, we obtain the series solution of the geodesics equation in the Schwarzschild geometry. The truncated series solution, containing only five terms, can reproduce the exact numerical solution with a high precision. In the first order of approximation we reobtain the standard expression for the perihelion precession. We study in detail the bending angle of light by compact objects in several orders of approximation. The extension of this approach to more general geometries than the Schwarzschild one is also briefly discussed.

General relativity is a very successful theory of the gravitational field, whose predictions are in excellent agreement with a large number of astronomical observations and experiments performed at the scale of the Solar System. In particular, three fundamental tests of general relativity, the perihelion precession of planet Mercury [1,2], the bending of light by the Sun [3,4], and the radar echo delay experiment [5,6] have all fully confirmed, within the range of observational/experimental errors, the predictions of Einstein's theory of gravity.
But the importance of these effects goes far beyond the limits of the Solar System. A fast full general relativistic method to simultaneously constrain the mass of massive black holes, their spin, and the spin direction by considering both the motion of a star and the propagation of photons from the star to a distant observer was developed in [7]. The spininduced effects on the projected trajectory and redshift curve of a star depend on both the value and the direction of the spin. The maximum effects over a full orbit can differ by a factor up to more than one order of magnitude for cases with significantly different spin directions. In [8] it was shown that the spin of the massive black hole at the Galactic Center can be constrained with 1σ error <∼ 0.1 or even greater than 0.02 by monitoring the orbital motion of a star with semi major axis <∼ 300 AU and eccentricity >∼ 0.95 over a period shorter than a decade through future facilities. An improvement in astrometric precision would be more effective at improving the quality of constraining the spin than an improvement in velocity precision. Short-period stars orbiting around the supermassive black hole in our Galactic Center can successfully be used to probe the gravitational theory in a strong regime [9]. By using 19 years of observations of the two best measured shortperiod stars orbiting our Galactic Center constraints on a hypothetical fifth force that arises in some extended theories of gravity or in some models of dark matter and dark energy were obtained in [9]. No deviations from General Relativity were found, and the fifth force strength wss restricted to an upper 95% confidence limit of |α| < 0.016 at a length scale of λ = 150 astronomical units. Moreover, 95% confidence upper limit on a linear drift of the argument of periastron of the short-period star S0-2 was obtained, a result that opens the possibility of testing gravitational theories using orbital dynamic in the strong gravitational regime of a supermassive black hole. The S-star cluster in the Galactic center allows the study of physics close to a supermassive black hole, including distinctive dynamical tests of general relativity [10], where a new and practical method for the investigation of the relativistic orbits of stars in the gravitational field near Sgr A* was developed, by using a first-order post-Newtonian approximation to calculate the stellar orbits with a broad range of periapse distance r p . For in depth discussions of the experimental and Solar System tests of general relativity see [11] and [12], respectively.
Due to its importance in many applications, the study of the motion of massive or massless particles in different geometries obtained as solutions of Einstein's gravitational field equations, and of their extensions, is a fundamental field of general relativity. The first exact solution of the vacuum field equations was the Schwarzschild solution [13], which can be used efficiently to explain all astronomical observations at the scale of the Solar System.
The exact equation of motion in Schwarzschild geometry is highly nonlinear, and therefore to obtain the observable physical parameters approximate methods must be used. The first order approximation of the equation of motion already gives the correct approximation of the perihelion precession of Mercury, and of the deflection of light by the Sun [11]. However, due to the importance of the problem many mathematical techniques for the study of the astrometric properties of the planetary motions and of the light have been developed. A standard approach is based on the solution of the Hamilton-Jacobi equation [13], where g µν are the components of the metric tensor, and m is the mass of the particle.
By representing S as S = −Et + Mϕ + S r (r), where E and M are the constants of the energy and angular momentum, one can obtain the full solution of the equation of motion in Schwarzschild geometry in an integral form as [13] ct = E mc 2 ϕ = Mdr The parameters of the orbits can be obtained by solving the integrals by using some approximate methods. The geodesic equations obtained from the Schwarzschild gravitational metric in the presence of a cosmological constant were solved exactly, and expressed in a closed form in [14] as where α S = 2GM/c 2 , ǫ is an arbitrary integration constant, ℘ is the Weierstrass function that gives the inversion of the elliptic integral In this approach the exact expression of the perihelion precession is , and e 1 is a root of the cubic equation in the integral. The perihelion precession and deflection of light have been investigated in the four-dimensional general spherically symmetric spacetime in [15], where the master equation has also been obtained. As an application of this master equation, the Reissner-Nordstorm solution and Clifton-Barrow solution in f (R) gravity have been investigated. The homotopy perturbation method, which was introduced in [16], was applied for calculating the perihelion precession angle of planetary orbits in General Relativity in [17,18]. The basic ideas behind the homotopy perturbation method are as follows [16]. We start from the nonlinear differential equation A(u) = g(r), r ∈ Ω , where A is a general differential operator, and g(r) is a known analytic function, with the boundary conditions B(u, ∂u/∂n) = 0; r ∈ Γ , where B is a boundary operator, and Γ is the boundary of the domain Ω . We assume that the operator A can be divided into two parts M and N, and we reformulate our initial equation as M(u) + N(u) = g(r). Then the homotopy v(r, p) : Ω × [0, 1] → IR is constructed in the following way: r ∈ Ω and p ∈ [0, 1] is an imbedding parameter, and y 0 is an initial approximation of the equation. Since 0 ≤ p ≤ 1, it can be considered as a small parameter, and one can assume that the solution of the equation can be expressed as a power series in p as v = ∞ i=0 v i p i When p → 1, then this series becomes the approximate solution of the equation, that is This series is generally convergent. The study and the applications of Adomian's Decomposition Method (ADM) [19][20][21][22], which allows to investigate the solutions of many kinds of ordinary, partial, stochastic differential and integral equations that describe numerous physical and/or mathematical problems has attracted a lot of attention in recent years. Historically, the ADM was first proposed and applied in the 1980's [23][24][25][26]. An essential advantage of the ADM is that with its use one can obtain analytical approximations to the solutions of a rather wide class of nonlinear (and stochastic) differential and integral equations, without the need of linearization, perturbation, closure approximations, or discretization methods. Usually the application of these methods could lead to the necessity of intensive numerical computation. Moreover, to make solvable and to obtain closed-form analytical solutions of a nonlinear problem implies the necessity of introducing some simplifying and restrictive assumptions.
It is important to mention that ADM can generate the solution of a given equation in the form of a power series. The terms of the series are obtained by recursive relations using the Adomian polynomials. Another important advantage of the ADM is that usually the series solution of the differential equation converges fast, and therefore the use of this method saves a lot of computational time. Moreover, in the ADM there is no need to linearize or discretize the differential equation. Reviews of ADM and its applications in applied mathematics and physics can be found in [19,20], respectively. Many studies have been devoted to the modification and improvement of the ADM in an attempt to increase its accuracy, and/or to extend the applicability of the initial method [21,22,[27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. An important improvement of the ADM s represented by the Laplace-Adomian Decomposition method [42], in which the Adomian Decomposition Method is applied not to initial equation, but to the Laplace transformed one.
Even that the ADM has been extensively used in the study of many problems in different fields of physics and engineering, it has been applied very little in astronomy, astrophysics or general relativity, the only exception from this "rule" known to authors being the papers [43] and [44]. It is the purpose of the present paper to investigate the equations describing the perihelion precession and light bending in general relativity for static gravitational fields by using the Adomian Decomposition Method, representing a very powerful mathematical method for the investigation of the solutions of nonlinear differential equations. For the geometry outside a compact, stellar type object (the Sun) we adopt some specific static and spherically symmetric vacuum solutions of general relativity. In particular the power series solution of the equation describing the motion of massive and massless particles in Schwarzschild geometry is investigated in detail. As a first step in our analysis we derive the equations of motion for particles in arbitrary spherically symmetric spacetimes, and we develop a general formalism for obtaining the equations of motion that can be used for any given metric. As the next step in our study we adopt the Schwarzschild form of the metric, and we apply the Laplace-Adomian decomposition method to obtain its approximate analytical power series solution for both massive and massless particles. We compare our solutions with the exact numerical solutions of the equations of motion, and it turns out that by truncating our series to five terms only we obtain a very good description of the solution of the equation of motion. Moreover, in the first approximation we can reobtain easily the standard expressions of the perihelion precession and the bending angle of light.
The present paper is organized as follows. We derive the equations of motion of massive and massless particles in arbitrary static spherically symmetric spacetimes in Section II.
The application of the Laplace-Adomian Decomposition Method to the case of second order nonlinear differential equations is presented in Section III. The power series solution of the equation of motion of massive particles by using the Laplace-Adomian Decomposition Method is obtained in Section IV, where the comparison with the exact numerical solution is also performed. The motion of photons in Schwarzschild geometry is investigated in Section V. Finally, in Section VI, we discuss and conclude our results.

SPACE-TIMES
In the following we will restrict our analysis to the case of static and spherically symmetric metrics, given by where we have denoted dΩ 2 = dθ 2 + sin 2 θdϕ 2 , and θ and ϕ are the standard coordinates on the three-sphere. The time variable takes real values only, t ∈ R while the radial coordinate r ranges over a finite open interval r ∈ (r min , r max ), so that −∞ ≤ r min ≤ r max ≤ ∞.
Moreover, we also require that the functions ν(r) and λ(r) are strictly positive, and that on the interval (r min , r max ) they are (at least piecewise) differentiable. This form of the metric is relevant for the study of the dynamics of particles (both massive and massless) in the Solar System. Important observational evidence for the correctness of the theory of general relativity is provided, at the level of the Solar System, by three fundamental tests, which also allow the testing of its extensions and generalization, as well as of alternative theories of gravitation.
These three essential tests are the perihelion precession of Mercury, the deflection of photons by the Sun, and the radar echo delay observations. These three effects have been successfully used to test the Schwarzschild solution of general relativity, as well as other predictions of the theory. However, it is also important to study these physical phenomena in arbitrary static spherically symmetric space-times for any given metric. In the present Section, we develop a formalism that can be used for obtaining the equations of motion, and compute the perihelion precession and light bending angle in any static spherically symmetric metric.
This formalism was first introduced to study the Solar System tests for some modified gravity vacuum solutions in [45][46][47][48].

A. The equation of motion of massive test particles
The geodesic equations of motion of a massive test particle in the gravitational field of the metric given by Eq. (5) can be derived with the use of the variational principle where by a dot we have denoted d/ds. It can be easily checked that the orbit is planar, and therefore without any loss of generality we can take θ = π/2. Hence ϕ is the only the angular coordinate in this problem. Since t and ϕ do not appear explicitly in Eq. (6), their conjugate momenta gives two constants of motion, denoted E and L, so that The constant E gives the energy of the particle, while the constant L is related to its angular momentum.
From the line element Eq. (5) we obtain the following equation of motion for r, By substitutingṫ andφ from Eqs. (7) gives the following relation, We introduce now a new variable u, defined as r = 1/u, as well as the transformation d/ds = Lu 2 d/dϕ. Then Eq. (9) takes the form We represent formally e −λ as e −λ = 1 − f (u), thus obtaining We take the derivative of the above equation with respect to ϕ, which gives where Eq. (12) gives the equation of motion of a particle in an arbitrary spherically symmetric geometry.

The precession of the perihelion
The root of the equation u 0 = F (u 0 ) gives a circular orbit with u = u 0 . Any deviation δ = u − u 0 from it can be obtained from the substitution of u = u 0 + δ into Eq. (12), which gives the equation In the first order of approximation , and hence Therefore, to first order in δ, the trajectory of the massive particle can be obtained as where δ 0 and β are two arbitrary constants of integration. The angles for which r is minimum are the angles of the perihelia of the orbit. Therefore they are determined from the condition that u or δ is maximum. Hence, from one perihelion to the next the orbital angle varies by a quantity ∆ϕ, given by The parameter σ introduced in the previous equation is called the perihelion advance.
From a physical point of view it represents the rate of advance of the perihelion after one rotation. As the test particle advances through ϕ radians in its orbit, its perihelion precesses by σ∆ϕ radians. From Eq. (17), σ can be expressed as or, for small (dF/du) u=u 0 , as For a complete rotation of the planet we obtain ϕ ≈ 2π(1 + σ). Hence the advance of the perihelion is To be able to obtain effective estimations of the perihelion precession we must know the expression of the angular momentum of the particle L as a function of the geometric parameters of the orbit. We will obtain now the expression of L in the Newtonian limit [45][46][47][48].
Let's assume that the planet moves on a Keplerian ellipse, with semi-axisā andb, respectively, whereb =ā √ 1 − e 2 , and by e we have denoted the eccentricity of the orbit. The ellipse has a surface area πāb. The oriented surface area of the ellipse is d A = ( r × d r) /2, and consequently the areolar velocity of the planet is given by On the other hand T can be obtained from Kepler's third law as T 2 = 4π 2ā3 /GM [49]. In the Newtonian limit of small velocities ds ≈ cdt, and the conservation equation of the angular momentum reduces to r 2 dϕ/dt = cL. Hence we obtain first L = 2πā 2 √ 1 − e 2 /cT , and hence

The equation of motion of massive particles in Schwarzschild geometry
As a first astronomical application of the formalism introduced in the previous Section we obtain the precession of the perihelion of a planet in the Schwarzschild geometry, with Then we immediately obtain f (u) = (2GM/c 2 ) u. On the other hand since for this geometry ν + λ = 0, we easily find and respectively. Therefore the equation of motion of a massive test particle in Schwarzschild geometry is given by The radius u = u 0 of the circular orbit is found as the solution of the algebraic equation with the only physically acceptable solution given by Therefore which is the standard general relativistic result [13].

B. Equation of motions of photons and the deflection of light
In a gravitational field a photon follows a null geodesic, given by ds 2 = 0. In this case the affine parameter along the trajectory of the photon can be taken as an arbitrary quantitỹ λ. In the following we denote again by a dot the derivatives with respect to it. Similarly to the case of the motion of massive particles we have two constants of motion, the energy E and the angular momentum L, which can be obtained from Eqs. (7).
The equation of motion of the photon is given bẏ By using the constants of motion the above equation can be transformed intȯ We change now the independent variable r to u = 1/r. With the use of the conservation equations we eliminate the derivative with respect to the affine parameter, thus obtaining We take the derivative of Eq. (31) with respect to ϕ, and thus we find the basic equation of the photon in an arbitrary static spherically symmetric geometry, as given by where we have denoted In the particular case of the Schwarzschild geometry we have ν + λ = 0 and f (u) = (2GM/c 2 ) u, giving P (u) = (2GM/c 2 ) u 3 and V (u) = (3GM/c 2 ) u 2 , respectively. Hence the equation of motion of photons in the Schwarzschild metric is obtained as

The deflection angle of light
In the lowest order of approximation we can neglect the term on the right hand side of Eq. (32). Then the solution is given by a straight line, where by R we have denoted the distance of the closest approach to the central massive gravitating object. In the next order of approximation Eq. (35) is substituted into the right-hand side of Eq. (32). Hence the equation of the trajectory is given by a second order linear inhomogeneous differential equation, which has a general solution given by u = u (ϕ). The photons travel towards the star from infinity at the asymptotic angle ϕ = − (π/2 + ε), and are deflected to infinity at the asymptotic angle ϕ = π/2 + ε. The angle ε can be computed by solving the algebraic equation u (π/2 + ε) = 0. For the total deflection angle of the photon beam we find δ = 2ε.
a. The deflection of light in Schwarzschild geometry. We consider now the case of the Schwarzschild geometry. In the lowest order of approximation from Eqs. (35) and (36) we obtain for the photon trajectory the second order linear differential equation having the general solution given by By substituting ϕ = π/2 + ε, u = 0 into Eq. (38) we easily find where we have used the simple trigonometric relations cos (π/2 + ε) = − sin ε, cos (π + 2ε) = − cos 2ε, and the approximations sin ε ≈ ε and cos 2ε ≈ 1, respectively. The total deflection angle of light in the Schwarzschild geometry is thus δ = 2ε = 4GM/c 2 R, a well-known result in general relativity [13]. Let's consider a nonlinear differential equation of the form where ω and b are constants, and g is an arbitrary nonlinear function of dependent variable With the use of the properties of the Laplace transformation we easily find With the use of the initial conditions for our problem we obtain As a next step in our analysis we assume that the solution can be represented in the form of an infinite series, where the terms y n (x) are computed recursively. As for the nonlinear operator g(y), it is decomposed as where the A n 's are the so-called Adomian polynomials, defined generally as [20] A The first five Adomian polynomials can be obtained in the following form, Substituting Eqs. (44) and (45) into Eq. (43) we obtain Matching both sides of Eq. (52) yields the following iterative algorithm for the power series solution of Eq. (40), ...
By applying the inverse Laplace transformation to Eq. (53), we obtain the value of y 0 .
Substituting y 0 into Eq. (47) to find the first Adomian polynomial A 0 . Then we substitute The other terms y 2 , y 3 , . . .,y k+1 , can be computed recursively in a similar step by step approach.

A. Power series solution of the equation of motion
Assume that the solution of Eq. (78) can be obtained in power series form, Now taking Laplace transform to Eq. (78) yields Hence we obtain We write down a few Adomian polynomials for U 2 , Substituting Eq. (79) and U 2 = ∞ n=0 A n (ϕ) into Eq. (82) gives the relation or, equivalently, Next we rewrite Eq. (88) in the recursive forms ..., With the help of the explicit expressions of the Adomian polynomials, we obtain The power series solution of the equation of motion of the massive test particles in Schwarzschild geometry is thus given by The absolute difference ∆(ϕ), defined as where e is the eccentricity of the orbit. The angular momentum L 2 can be expressed in terms of the geometric parameters of the orbit by using Eq. (21). Thus in physical units we Therefore for this choice of parameters and initial conditions the successive terms in the power series solution of the equation of motion can be obtained by the Laplace Adomian decomposition method as follows: (116)

The perihelion precession
In the first order approximation the solution of the equation of motion for massive particles in Schwarzschild geometry is obtained as By neglecting the constant terms and those who oscillate through two cycles on each orbit we obtain A. Power series solution The truncated power series solution of Eq. (120) can be obtained immediately by using the Laplace-Adomian decomposition method as a 4 180ϕ sin ϕ + 60ϕ sin(2ϕ) + 119 cos ϕ + 96 cos(2ϕ) + 9 cos(3ϕ) + By taking ϕ = π/2 + ǫ, we have By performing a first order series expansion with respect to ǫ, and taking u = 0, gives the light deflection angle as which in the limit GM/c 2 R << 1 reduces to the well-known general relativistic result In the second order of approximation, with u ≈ u 0 + u 1 + u 2 , we obtain by using the same procedure Finally, in the fourth order of approximation we obtain In the case of the vacuum spacetimes outside charged spherically symmetric objects carrying a charge Q the geometry is given by the Reissner-Nordstrom metric, In the case of rotating star the first order monopole correction to the Schwarzschild metric is given by [50] e ν = e −λ = 1 − 2GM where J = IΩ, where I is the moment of inertia, Ω is the angular velocity, and J is the angular momentum of the star, respectively. The physics of the particle motion in these metrics can be easily investigated by using the methods developed in the present paper. For example, in the case of the Reissner-Nordstrom metric, f (u) = 2GMu/c 2 − Qu 2 , and respectively. Therefore the equation of motion in the Reissner-Nordstrom metric is given by which can be solved easily by using the Laplace-Adomian Decomposition Method. The radius of the circular orbit u 0 can be obtained as a solution of the algebraic equation In the first order of approximation, and with the assumption Q/L 2 ≪ 1, u 0 can be approximated as u 0 ≈ GM/c 2 L 2 . Thus, the perihelion precession ∆φ in the Reissner-Nordstrom geometry is given by .
The first term in Eq. (138) gives the general relativistic correction term for the perihelion precession, while the second term gives the correction due to the presence of the charge. In the case of the motion of massless particles we have and Hence the equation of motion for photons is given by This equation can also be investigated easily by using the Laplace Adomian Decomposition Method, and the methods developed in the present paper.