Perihelion advance and trajectory of charged test particles in Reissner-Nordstrom field via the higher-order geodesic deviations

By using the higher-order geodesic deviation equations for charged particles, we apply the method described by Kerner et.al. to calculate the perihelion advance and trajectory of charged test particles in the Riessner-Nordstrom spacetime. The effect of charge on the perihelion advance is studied and compared the results with those obtained earlier via the perturbation method. The advantage of this approximation method is to provide a way to calculate the perihelion advance and orbit of planets in the vicinity of massive and compact objects without considering Newtonian and post-Newtonian approximation.


Introduction
The problem of planets motion in general relativity is the subject of many studies in which the planet has been considered as a test particle moving along its geodesic [2]. Einstein made the first calculations in this regard for the planet Mercury in the Schwarzschild space-time which resulted in the equation for the perihelion advance where G is the gravitational constant, M is the mass of the central body, a is the length of semi-major axis for planet's orbit and e is eccentricity. Derivation of perihelion advance by using this method leads to a quasi-elliptic integral whose calculation is very difficult, which is then evaluated after expanding the integrand in a power series of the small parameter GM/rc 2 . For the low-eccentricity trajectories of planets, one can obtain the following approximate formula for the perihelion advance ∆ϕ = 6πGM a(1 − e 2 ) ≃ 6πGM a (1 + e 2 + e 4 + e 6 + · · ·), even for the case of Mercury up to second-order of eccentricity, the perihelion advance differs only by 0.18% error from its actual value [1]. It should be noted again that Einstein's method is only valid for the small values of GM/rc 2 . . Also, u µ is the unit tangent vector to the central world-line, n µ is the separation vector to the curve s = const and N µ = b µ − Γ µ λρ n λ n ρ is the second-order geodesic deviation [3].
In what follows, we show that one can obtain the same results (without taking the complex integrals) only by considering the successive approximations around a circular orbit in the equatorial plane as the initial geodesic with constant angular velocity, which leads to an iterative process of the solving the geodesic deviation equations of first, second and higher-orders [4,5,6]. Here, instead of the GM/rc 2 parameter the eccentricity, e, plays the role of the small parameter which is controlling the maximal deviation from the initial circular orbit. In this method, we have no constraint on GM/rc 2 anymore. So, one can determine the value of perihelion advance for large mass objects and write it in the higher-order of GM/rc 2 .
The orbital motion of neutral test particles via the higher-order geodesic deviation equations for Schwarzschild and Kerr metrics are studied in [1] and [5], respectively. Also, for massive charged particles in Reissner-Nordstrom metric, geodesic deviations have been extracted up to first-order [7]. In this paper, by using the higher-order geodesic deviations for charged particles [8], we are going to obtain the orbital motion and trajectory of charged particles. We also expect that our calculations reduce to similar one in Schwarzschild metric [1] by elimination of charge. In fact, we generalize the novel method used in reference [1] for neutral particles in the Schwarzschild metric to the charged particles in the Reissner-Nordstrom metric. Recently, an analytical computation of the perihelion advance in general relativity via the Homotopy perturbation method has been proposed in [9]. Also, one can study the perihelion advance of planets in general relativity and modified theories of gravity by using different methods in [10]- [24].
The structure of the paper is as follows: In section 2, by using the approximation method introduced in reference [8], we derive the higher-order geodesic deviation for charged particles. By using the first-order geodesic deviation equations, the orbital motion of charged particles is found in section 3. In section 4, we obtain the second-order geodesic deviations and derive the semi-major axis, eccentricity, and trajectory using the Taylor expansion around a central geodesic. The obtained results are discussed in section 5.
2 The higher-order geodesic deviation method As it mentioned above, the higher-order geodesic deviation equations for charged particles have been derived in [8] for the first-time. In this section, we are going to derive the geometrical set-up used in our work. The geodesic deviation equation for charged particles is [7] where D Ds is the covariant derivative along the curve and n µ is the separation vector between two particular neighboring geodesics (see figure 1). Here, u µ is the tangent vector to the geodesic, R µ λνκ is the curvature tensor of space-time, q and m are charge and mass of particles (particles have the same charge-to-mass ratio, q/m), and F µ ν is the electromagnetic force acting on the charged particles. For neutral particles, the above equation reduces to the following geodesic deviation [25,26] which is the well-known equation (Jacobi equation) in general relativity. We introduce the fourvelocity u α (s, p) = ∂x α ∂s as the time-like tangent vector to the world-line and n α (s, p) = ∂x α ∂p as the deviation four-vector as well. Practically it is often convenient to work with the non-trivial covariant form. It can be obtained by replacement of the trivial expressions for the covariant derivatives, the Riemann curvature tensor and use of the equation of motion in the left-hand side of equation (3) [7] The geodesic deviation can be used to compose geodesics x µ (s) near a given reference geodesic x µ 0 (s), by an iterative method as follows. Considering this, one can write Taylor expansion of x µ (s, p) around the central geodesic and obtain the first-order and higher-order geodesic deviations for charged particles our aim is to obtain an expression in terms of the deviation vector. As shown in the above equation, the second term, ∂x µ ∂p , is the definition of deviation vector and shows the first-order geodesic deviation. But in the third term, ∂ 2 x µ ∂p 2 , is not vector anymore. Therefore, we define the vector b µ as follow to change ∂ 2 x µ ∂p 2 into the expression showing the second-order geodesic deviation. By substituting equation (7) into equation (6), one can obtain the expression in terms of the order of vector deviation In the above expression, one can make some changes for simplification. We consider δ n x µ (s) as n-th order of geodesic deviation and by assuming (p − p 0 ) as a small quantity, ǫ, we rewrite equation (8) as follows where δx µ (s) = ǫn µ (s) is the first-order geodesic deviation and δ 2 x µ (s) = ǫ 2 (b µ − Γ µ νλ n ν n λ ) is the second-order geodesic deviation. In order to obtain the second-order geodesic deviation equation, one can apply the definition of the covariant derivative on equation (7)(for more details see reference [8] and appendix therein) Similar to the first-order geodesic deviation equation (5), we can write equation (10) in the nonmanifest covariant form As it clears, the left-hand side of the second-order geodesic deviation equation (11) is same to the left-hand side of equation (5). As in the case of the second-order geodesic deviation, the higherorder geodesic deviation equations have the same left-hand side and different right-hand side. A non-manifest covariant form of the third-order geodesic deviation equation is given in appendix 1.
The successive approximations to the exact geodesic (b) have been shown in figure 1. Lines (c) and (d) represent the first-order approximation i.e.
In the next section, we are going to obtain the components of n µ from the first-order geodesic deviation, equation (5), for a circular orbit of charged particles. Then by substituting them into equation (11), we can solve the second-order geodesic deviation equations, b µ . Finally, by substituting n µ and b µ into equation (8), we will find the relativistic trajectory of charged particles in Reissner-Nordstrom space-time.
3 The first-order geodesic deviation 3

.1 Circular orbits in Reissner-Nordstrom metric
The Reissner-Nordstrom metric is a static exact solution of the Einstein-Maxwell equations which describes the space-time around a spherically non-rotating charged source with mass M and charge Q (in the natural coordinate with c = 1 and G = 1) where Also, the vector potential and the electromagnetic field of Maxwell's equations is [7] By assuming that M 2 > Q 2 , we are going to obtain the equation of motion for test particles which have mass m and charge q. Now, we consider a circular orbit with a constant radius R. On the other hand, we know that the angular momentum of particles which are bounded to the spherically symmetric condition is limited to the equatorial plane. For this purpose and for simplicity, we limit the space to the plane of θ = π 2 in which the angular momentum is in the z direction. By using of the Euler-Lagrange equation, one can lead to the following constants of motion where l = J m is the angular momentum per unit mass,φ = ω is the angular velocity and ε is the energy per unit mass.
Finally, from equations (12), (14) and (15) one obtains two constraints namely, the conservation of the absolute four-velocity and the radial acceleration. Now, due to the fact that the radius of the circular orbit is constant (r = R), two mentioned constraints vanish at all times and this creates two relations between R, l, and ε as follows and As we expect that by eliminating charge, all obtained equations reduce to the similar equations in the Schwarzschild metric. In summary, we obtain the following four-velocity vector for a circular orbit with radius R in an equatorial plane In the next subsection, we obtain the orbital motion by using the higher-order geodesic deviation method and compare the results with the perturbation method.

First-order geodesic deviation around the circular orbits
Now let us calculate the first-order geodesic deviation for the components n t , n r , n θ and n φ , by using of equation (5) where the matrix elements are given by which is similar to the Schwarzschild case. So in this case we can neglect this solution (n θ = 0), because the new plane of orbit is a new one inclined, or just a change of coordinate system [5]. Now, by eliminating the derivatives of n t and n φ in the differential equation of n r , we obtain the following oscillating equation d 2 n r ds 2 + w 2 n r = 0, with the characteristic frequency By considering n r 0 > 0, we choose the following solution for n r n r = −n r 0 cos(ωs).
Also, from the n t and n ϕ geodesic deviation equations, the solutions for n t and n ϕ are given by where the amplitudes depend on n r In this way, the trajectory and the law of motion are obtained by r = R − n r 0 cos(ωs), where the argument phase of the cosine function is taken by s = 0 for perihelion and s = π ω for aphelion. Now, (28) can be written as By direct solution of the Euler-Lagrange equations, the trajectory of motion for particles is obtained in terms of centrifugal inertia [27] Obviously, equations (31) and (32) show that we have the same results. It means that if we bring up the eccentricity e to n r 0 R and the semi-major axis a to R, the same results are extracted, but there is also a difference that the circular frequency, ω, is lower than the circular frequency of the unperturbed circular motion, ω 0 . So, if the circular frequency decreases, the period increases. Then we obtain an expression for the periastron shift per one revolution as It can be seen from above equation that the charge parameter, Q, decreases the perihelion advance.
In the perturbation method (Einstein's method), the orbital motion for charged particles moving in the equatorial plane of the Reissner-Nordstrom source is given by [23] comparing equation (33) with equation (34) shows that the presented method can be used in the vicinity of very massive and compact objects which is having a non-negligible ratio of M R . When the source is neutral and for the small values of M R , equation (33) reduces to the standard formula for Perihelion advance of planets [26]. If we also compare equation (33) to equation (2), it is clear that in the first-order deviation, we hold only the terms up to e 2 . In order to obtain the △ϕ for the higher values of the eccentricity, we must go beyond the first-order deviation equations. Therefore in the next section, we solve the second-order geodesic deviation equations in Reissner-Nordstrom space-time.

The second-order geodesic deviation
In this section, by using the first-order geodesic deviation equation and inserting equations (23), (24) and (25) into equation (10) and also doing a set of hard calculations, a linear equations system for the second-order geodesic deviation vector b µ is obtained where the constants C t , C t q , C r , C r q , and C ϕ , C ϕ q contain quantities depending on M , R, ω , ω 0 , q and Q and Here we do not have used any approximation in C i , (i = r, θ, φ) but in what follows we neglect terms of higher order of the small parameters M R , Q M and q m . Solving the matrix equation (10) for b µ is similar to the approach used in the previous section (for the first-order geodesic deviation vector n µ ) which contains the terms with characteristic frequency ω. Here we are only interested in a particular solution because of the oscillating general solution with the angular frequency ω already take into account for n µ (s). The particular solution of the above equation which is containing the oscillating terms with the angular frequency 2ω, the linear terms in the proper time s and constants. To obtain the trajectory x µ according to the equation (9), we need to calculate 1 2 δ 2 x µ . Also for x µ , the perihelion is extracted by ωs = 2kπ and the aphelion is derived by ωs = (1 + 2k)π where k ∈ Z.
In appendix 2, we have put the particular solution of the above equation, b µ , the second-order geodesic deviation δ 2 x µ , and the semi-major axis a and eccentricity e, respectively.
Finally, successive approximation brings us to trajectory by substituting s(ϕ) to ϕ(s) In the Schwarzschild limit, we have an elliptical orbit with [1] a = R + (n r 0 ) 2 R Also, for the Schwarzschild case the shape of the orbit is described up to second-order of ( which is in agreement with equation (62) of reference [1].

Third-order geodesic deviation and Poincaré-Lindstedt's method
In the previous section, we have calculated the trajectory of charged particles up to second-order. To find a more accurate trajectory, we need to obtain the higher-order terms of expansion (9). Using the first and second-order solutions and third-order equation (51) for δ 3 x µ , we have where m ij are defined in equation (19) and the coefficients As one can see the right-hand side of equations (46) have a frequency that is same as the eigenvalues of the differential matrix in the left-hand side (resonant terms). This makes a new problem i.e. infinite solution for δ 3 r which is called the secular term (growing without bound). For avoiding these unbounded deviations we use the Poincaré's method. In this method by replacing ω by infinite series in power of the infinitesimal parameter ǫ = n r 0 R as the correction frequencies ω 1 , ω 2 , ω 3 , · · · can be chosen such that the Poincaré's resonances vanish. By considering a differential equation for x µ as Now, by developing both of the sides in terms of a series of the parameter ǫ, for avoiding the secular terms, we find some algebraic relations on ω 1 , ω 2 , ω 3 , · · ·. In the Schwarzschild limit, we have [5] where ǫ 2 = (n r 0 ) 2 R 2 . The resonant terms will also appear at the fifth-order approximation, by terms cos 5 (ws), sin 3 (ws) cos 2 (ws), etc, this problem can be solved in a similar way.
Finally, we note that the electric charge of any celestial body being practically close to zero anyway. Therefore, it is worth to investigate the geodesic deviation and higher-order geodesic deviations in a more realistic background such as the Schwarzschild metric in a strong magnetic dipole field or magnetized black holes [28,29,30]. The study of them will be the subject of the future investigations.

Conclusion and discussion
Many of significant successes in general relativity are obtained by approximation methods. One of the most important approximation scheme in general relativity is the post-Newtonian approximation; an expansion with a small parameter which is the ratio of the velocity of matter to the speed of light. A novel approximation method was also proposed by Kerner et. al. which is based on the world-line deviations [1].
The calculation of the perihelion advance by means of the higher-order geodesic deviation method for neutral particles in different gravitational fields such as Schwarzschild and Kerr metric was first studied in several papers [1,5]. In the present paper by using of the higher-order geodesic deviation method for charged particles [8], we applied this approximation method to charged particles in the Reissner-Nordstrom space-time.
We first started with an orbital motion which is close to a circular orbit with constant angular velocity which is considered as zeroth-approximation (unperturbed circular orbital motion) with the orbital frequency ω 0 . In the next step, we solved the first and second-order deviation equations which reduced to a system of the second-order linear differential equations with constant coefficients. The solutions are harmonic oscillators with characteristic frequency. From equation (35), the first and second-order corrections are oscillating terms with angular frequency ω and 2ω, respectively.
Finally, we have obtained the new trajectory by adding the higher-order geodesic deviations (non-linear effects) to the circular one, equation (42). The advantage of this approach is to get the relativistic trajectories of planets without using Newtonian and post-Newtonian approximations for arbitrary values of quantity M/R.

Acknowledgements
We wish to express our thanks to the anonymous referees for their comments and suggestions that helped us significantly to improve this manuscript.

Appendix 1
For solving the third-order geodesic deviation equation, we should invoke to the Poincare's method. For this purpose, it is better to write the third-order geodesic deviation as δ 3 x µ . The third-order geodesic deviation equation δ 3 x µ is related to the third-order geodesic deviation vector h µ where h µ = Db µ Dp . We derive the third-order geodesic deviation equation as by substituting δ 3 x µ in term of h µ into above equation, we obtain equation (72) for case q = 0 [1].

Appendix 2
The second-order geodesic deviation vector b µ is As explained in section 2 the second-order geodesic deviation, (δ 2 x = b µ − Γ µ νλ n ν n λ ), are given by also, the semi-major axis a and eccentricity e are a = R + (n r 0 ) 2 R in which for massive central objects we have neglected all terms of order qQ 4πmM .