Higher order corrections to deflection angle of massive particles and light rays in plasma media for stationary spacetimes using the Gauss-Bonnet theorem

The purpose of this article is twofold. First, we extend the results presented in [Gabriel Crisnejo and Emanuel Gallo, Phys.Rev.D 97, 124016 (2018)] to stationary spacetimes. Specifically, we show that the Gauss-Bonnet theorem can be applied to describe the deflection angle of light rays in plasma media in stationary spacetimes. Second, by using a correspondence between the motion of light rays in a cold non magnetized plasma and relativistic test massive particles we show that this technique is not only powerful to obtain the leading order behavior of the deflection angle of massive/massless particles in the weak field regime but also to obtain higher order corrections. We particularize it to a Kerr background where we compute the deflection angle for test massive particles and light rays propagating in a non homogeneous cold plasma by including third order corrections in the mass and spin parameters of the black hole.

The purpose of this article is twofold. First, we extend the results presented in [1] to stationary spacetimes. Specifically, we show that the Gauss-Bonnet theorem can be applied to describe the deflection angle of light rays in plasma media in stationary spacetimes. Second, by using a correspondence between the motion of light rays in a cold non magnetized plasma and relativistic test massive particles we show that this technique is not only powerful to obtain the leading order behavior of the deflection angle of massive/massless particles in the weak field regime but also to obtain higher order corrections. We particularize it to a Kerr background where we compute the deflection angle for test massive particles and light rays propagating in a non homogeneous cold plasma by including third order corrections in the mass and spin parameters of the black hole.

I. INTRODUCTION
General relativity is a geometric theory of gravitation which predicts the existence of black holes (BHs). BHs are extremely important astrophysical objects characterised by many astrophysical processes. For example, they are important to explain the formation of extragalactic jets from a black hole accretion disk and tidal disruption of stars by supermassive black holes, to study the deflection of light and particles, but also to explore various predictions of quantum field theory and quantum gravity through Hawking radiation. The recent announcement of the first image concerning the detection of an event horizon of a supermassive black hole at the center of the neighboring elliptical M87 galaxy by the Event Horizon Telescope (EHT) collaboration is a great triumph of general relativity [2][3][4][5][6][7]. In general, the image of a black hole with a surrounding accretion disk appears distorted due to the strong gravitational lensing effect. In this way, black holes are expected to cast shadows on the bright background which is related to the existence of an event horizon and therefore an unstable photon region [8]. The shadow image is of great scientific significance, it can help us to probe the geometrical structure of the event horizon and perhaps to measure the angular momentum parameter of the black hole.
The gravitational deflection of light and particles is not only an important tool in astrophysics to test the existence of BHs, or to distinguish BHs from other compact objects but also to characterize the matter distribution at galactic or extragalactic scales. Quite amazingly, it turns out that the deflection angle can be linked with the topology of the spacetime by means of the Gauss-Bonnet theorem (GBT) applied to an associated optical geometry. Namely, Gibbons and Werner [9] showed that, in asymptotically flat and static spacetimes the deflection angle can be obtained by integrating the Gaussian optical curvature over an infinite domain of integration outside the massive body. This result was generalized to stationary spacetimes by Werner applying the Finsler-Randers geometry [10]. On the other hand, another method was applied by Ishihara et. al. [11] in order to calculate the deflection angle using the Gauss-Bonnet theorem. This method was generalized to stationary spacetimes by Ono et.al. [12,13].
Recently a new step forward was put by Crisnejo and Gallo [1] extending the use of the Gibbons-Werner method to more general situations where light rays propagate in a plasma environment. In that work the deflection angle of test massive particles is also studied using the Gauss-Bonnet theorem by taking into account a known correspondence between the dynamic of light rays in an homogeneous cold nonmagnetized plasma and test massive particles following geodesics at the same spacetime [14]. This correspondence was extensively used by other authors in the past. In [15] Bisnovatyi-Kogan and Tsupko used this correspondence to calculate the deflection angle of massive particles in static spacetimes in the weak-field approximation while in [16] and [17] the correspondence was used to study the strong deflection limit.
In [18] it was established for the first time a correspondence between the motion of light rays in a particular non-homogeneous plasma and the one of relativistic test charged massive particles in vacuum. In the same work, the Gauss-Bonnet method was applied to obtain an expression for the deflection angle in terms of the components of the energy-momentum tensor in a plasma environment generalizing in this way previous works restricted to the pure gravity case [19][20][21][22].
In many astrophysical situations the effect of the plasma medium on the light rays can be neglected but in others, as for example in the case of light rays propagating in the radio frequency band this is not suitable. A well know example is the influence of the solar corona on light rays propagating near to the Sun. The first works in this regard were made by Muhleman et. al. [23,24] who studied the influence of the solar corona on the time delay and in the light deflection in the weak field approximation. This is a case where a non-magnetized plasma is a good model of the medium but in other scenarios could be necessary to consider the presence of a magnetic field, that is, to consider a magnetized plasma. A rigorous derivation of the relation dispersion in this kind of plasma was presented by Breuer and Ehlers [25][26][27].
In this work we will only consider cold non-magnetized plasma for which a very detailed treatment was done by Perlick in [28]. In the same reference, Perlick derived an exact closed formula for the deflection angle of light rays in a Kerr metric surrounded by an arbitrary nonhomogeneous plasma. That expression is found in integral form in terms of the closest approach distance and is not limited to weak field regime.
The purpose of this paper is to extend the results presented in [1] to stationary spacetimes. To do so we will make use of a Finsler-Randers type metric. Although many papers in the literature applies the GBT to different stationary metrics, the problem of higher order correction terms in the deflection angle in stationary spacetimes is not in general discussed. Here we fill this gap by computing in a successfully way the deflection angle in the weak field regime including up to third order correction.
The plan of the paper is as follows: In Sec. II, we shall use the dynamic of light rays in a plasma medium to introduce the Finsler-Randers metric. In Sec. III, we review the GBT in stationary spacetimes. In Sec. IV, we obtain the deflection angle in stationary spacetimes surrounded by a plasma medium at linear order using the approach proposed by Ono, Ishihara and Asada in [12]. In Sec. V we extend the Werner's method to calculate the deflection angle of massive particles. In Sec. VI we study the problem of higher order corrections to the deflection angle by studying the deflection angle of massive particles and also light rays propagating in a nonhomogeneous plasma medium in a Kerr background. Finally in Sec. VII we comment on our results. Two appendices with relevant calculations are also included.
In the following we shall use natural units G = c = = 1 and we adopt the index conventions that Greek indices run from 0 to 3, early Latin indices (a,b,c,...) run from 1 to 3 and middle Latin indices (i,j,k,...) run from 1 to 2.

II. OPTICAL METRIC FOR STATIONARY AXISYMMETRIC SPACETIMES IN A PLASMA ENVIRONMENT
Let us consider a stationary axisymmetric spacetime (M, g αβ ) filled with a cold non-magnetized plasma described by the refractive index n, where ω(x) is the photon frequency measured by an observer following the integral curves of a timelike Killing vector field, and ω e (x) is the plasma frequency given by, where e and m e are the charge and the mass of the electron, respectively while N (x) is the number density of electrons in the plasma. Following the notation introduced in [29] and since then used in other works by other authors we have used the definition, in Eqn. (2). Note that only light rays with ω 2 (x) > ω 2 e (x) propagate through the plasma.
The dynamic of the light rays in this kind of medium is usually described through the Hamilton's equations associated with the Hamiltonian, with the dispersion relation, As explained by Perlick in [28], this Hamiltonian define a ray-optical structure on M and even more, under appropriated conditions it is possible to build a reduced ray-optical structure on a 3-dimensional manifoldM in order to discuss the spatial paths of light rays. In this section we will carry out the reduction process described in [28] (see Theorem 6.5.1 of that reference for a complete description of this process) to build this reduced ray-optical structure for an arbitrary stationary and axisymmetric spacetime given by the line element, To do that, we have to restrict our attention to the region on the spacetime where the Killing vector field ∂ ∂t is timelike and the propagation of light rays is allowed, i.e ω 2 (x) > ω 2 e (x). First of all we want to remark that the Hamiltonian (4) can be rewritten as follows, where Ω 2 := (g 03 ) 2 − g 00 g 33 As the metric is stationary, p 0 is a constant of motion which can be identified with the frequency measured by a stationary observer at infinity, that is, p 0 := −ω ∞ . Because of the basic equation from which the dynamic of light rays can be deduced is the dispersion relation (5), we can multiply the Hamiltonian (7) by a non-zero function without altering the dynamic of the light rays 1 . Then, we have wherẽ and as we will see Ω is a non-zero function within the region of interest. In particular, by making use of the following identities, (g 03 ) 2 − g 00 g 33 we can express the Hamiltonian (10) as follows, with As discussed in [28] a reduced ray-optical structure given by the HamiltonianH can be simply obtained by replacing the conserved momentum coordinate p 0 by the constant −ω ∞ , i.e. the photon frequency measured by an asymptotic observer. Then, the 3-dimensional reduced ray-optical structure is defined by the HamiltonianĤ , where being n the refractive index of the medium given by Eqn. (1) and where the gravitational redshift has already been taken into account, In this way we can re-write the Hamiltonian (15) as follows,Ĥ wherê and with inverseĝ ab (defined asĝ abĝ bc = δ c a ) given by, To obtain the above expression we have used the identity, On the other hand,β a are the components of the oneformβ given by,β It is important to note thatĝ ab is a 3-dimensional Riemannian metric and as follows from the Hamilton equations associated to (18) the motion of light rays are determined by [28], x a +Γ a bcẋ bẋc =ĝ ad ∂ eβd − ∂ dβe ẋ e whereΓ a bc are the Christoffel symbols of the metricĝ ab . From (24) we can see that the ifβ = 0 light rays of the reduced optical structure in general do not follow geodesics with respect to the metricĝ ab . For static spacetimes, they are indeed geodesics ofĝ ab , and it is known as Jacobi metric.
Finally, the optical path length takes the form, where is a Finsler-Randers type metric. According to the Fermat principle light rays between two points of the reduced ray-optical structure are the extremes of the functional (25) and in particular they follow geodesics with respect to the Finsler-Randers metric F (see [42] for a detailed treatment on Finsler geometry). We want to point out that an explicit expression for the Finsler-Randers metric (26) was obtained for the particular case of photons propagating in a cold non-magnetized plasma over a Kerr spacetime by Perlick in [28]. Expressions (20) and (22) are their generalization to any stationary spacetime, and in fact even when to our knowledge they were not explicitly presented before in the literature, they are already implicit in [28]. On the other hand, it was recently reported in [43] an expression of the Finsler-Randers metric for massive particles propagating in a Kerr spacetime. Expressions (20) and (22) also contain this particular case through the known correspondence between the motion of test massive particles in a given stationary axisymmetric spacetime and photons moving in a homogeneous plasma environment in the same background. We will use this correspondence later to discuss the motion of test massive particles.

III. DEFLECTION ANGLE IN STATIONARY SPACETIMES USING THE GAUSS-BONNET THEOREM
As it is well known, the Gauss-Bonnet theorem connects the intrinsic geometry of a two-dimensional surface with its topology.
Precisely, it can be enunciated as follows. Let D ⊂ S be a regular domain of an oriented two-dimensional surface S with a Riemannian metricg ij , the boundary of which is formed by a closed, simple, piece-wise, regular, and positive oriented curve ∂D : R ⊃ I → D. Then, (27) where χ(D) and K are the Euler characteristic and Gaussian curvature of D, respectively; k g is the geodesic curvature of ∂D; and ǫ i is the exterior angle defined in the ith vertex, in the positive sense.
In order to apply the Gauss-Bonnet theorem to calculate the deflection angle of light rays propagating in the equatorial plane in stationary spacetimes we will follow the procedure described first by Gibbons and Werner in [9] for the pure gravity case and extended by two of us to the plasma case [1]. That is, we will apply the Gauss-Bonnet theorem to a domain D R as in the Fig. 1 but at difference with the aforementioned works, we will consider in this case that light rays do not necessarily follow geodesics with respect to the optical metric g opt ij 2 . In all the cases we are going to consider it is assumedthat light rays are propagating in the equatorial plane of an stationary axisymmetric spacetime. This implies in particular 2 Later, we will identify this optical metric with the metricĝ ij .
The points S and R represent the source and receiver position while γp indicates a light ray emitted by the source and that reaches the observer. b is identified with the impact parameter and CR is a semi-circle of radius R connecting S and R.
that the geodesic curvature of the light rays projected to the optical manifold is not zero. We obtain the following expression for the bending angle α, where the limit R → ∞ in both sides of the equation is understood. As was mentioned in the past, for asymptotically flat spacetimes one has that [k g dσ dϕ ] CR → 1 as the radius of C R goes to infinity which allows to write the deflection angle in an even more simple way, In this work we will restrict our attention to asymptotically flat spacetimes and by then the expression (29) will be enough to calculate the deflection angle. We mention that the expression (29) is similar to the one obtained by Ono et. al. (see Eqn. (30) in [12]) where the difference is in the domain of integration because they are computing finite distance corrections to the deflection angle and in the present paper we are not interested in that case.
It has been reported in the literature two different ways to implement the equation (29) to calculate the deflection angle in a stationary and axysimmetric spacetime for the vacuum case. Here, we briefly describe both of them. One method (developed in [12]) consists in using the Riemmanian metricĝ ab restricted to the equatorial plane to compute the Gaussian and geodesic curvatures in (29) to calculate the bending angle. In that case due to light rays do not follow geodesics with respect to this metric, the geodesic curvature associated with γ p will not be zero. In particular, for light rays moving in the equatorial plane characterized by θ = π/2, its geodesic curvature can be calculated as follows, whereĝ is the determinant ofĝ ab . The other method which was reported by Werner in [10] for stationary spacetimes in the pure gravity case consists in building an osculating Riemannian metric from the Finsler-Randers one by using the Nazim's method [44]. The advantage of this method is that the light rays follow geodesics with respect to the resulting Riemannian metric and by them the second term in (29) vanishes. The disadvantage is that the computations are more cumbersome even at linear order.
In the following two sections we will compute the deflection angle of light rays propagating in a plasma environment from the equation (29) in the weak field approximation.

IV. DEFLECTION ANGLE IN WEAK FIELD APPROXIMATION USING ONO ET.AL. APPROACH
We are interested in the computation of the deflection angle in the weak field regime for an stationary and axisymmetric spacetime surrounded by a cold-plasma medium using the method proposed in [12]. In particular, in this Section we will obtain a general expression for the deflection angle taking into account first order corrections. That is, we will consider an stationary metric written as and we will keep only those terms which are linear in the bookkeeping parameters ǫ 1 , ǫ 2 and ǫ 3 . In general is assumed that h tϕ is proportional to a spin parameter a.
We will restrict all the discussion to the equatorial plane θ = π/2. At the considered order the Finsler-Randers metric is determined bŷ For a cold non magnetized plasma with a charge density profile of the form N (r) = N 0 + N 1 (r), with N 0 = constant, it follows that the refraction index can be written as where we have used that ω 2 e = K e N (r) = ω 2 e0 + K e N 1 (r). The deflection angle is determined by: Let us start with the computation of the Gaussian curvature at the considered order. It is given by: Let us assume thatK ǫ2 ,K ǫ1 andK plasma are different from zero, i.e., neither h rr , rh ′ tt or rN ′ 1 take constant values. In the following and the rest of the paper, we will only consider prograde orbits, that is orbits whose orbital angular momentum are in the same direction as the spin of the considered metric.
As we are interested in the first order contribution to the deflection angle it is enough to integrate along a domain D r where the curve γ p is approximated by the path followed by light rays in the flat background, that is, straight lines. This is called the Born approximation. Therefore the equation for γ p in the integration domain D r can be replaced by the flat geodesic straight line given by r ϕ = b sin ϕ , while the polar angular coordinate ϕ runs from 0 to π.
In that situation, the integration of (35) gives (39) As we will see in Sec. VI A, for higher order corrections this approximation should be improved.
Let us see α kg , The geodesic curvature k g obtained from (30) is given by, On the other hand for prograde orbits, Then, Finally the following expression for the deflection angle results, For static spacetimes (ǫ 3 = 0) this expression can be compared to the one obtained by Bisnovatyi-Kogan and Tsupko in [15] (see Eqn. (30) on that reference). The expression (44) is given in radial coordinates instead of the Cartesian ones used in [15] but after a change of coordinates and integration by part it can be shown that they are completely equivalent.
The case of retrograde orbits can be analyzed in a similar way under a correspondent domain D ′ r , resulting in a similar expression to (44) but with an opposite sign in the last term. Alternatively, it can be understood as a result that the motion of a retrograde photon (with orbital angular momentum p ϕ < 0) in the original coordinate expression for the metric (with positive spin a and starting in the asymptotic region in ϕ S = 0 and with the receiver position at ϕ R = −π) is equivalent to the motion of a photon with positive angular momentum p ϕ > 0 in a metric of the same coordinate form as (31) but obtained by replacing a by −a and studying orbits starting in the asymptotic region in ϕ S = 0 and ending in ϕ R = π.
It is worthwhile to note that unlike the exact expression for the deflection angle in terms of the closest approach distance in the Kerr spacetime for nonhomogeneous plasma derived by Perlick in [28], in the practice it is usually more convenient to work with an approximate expression in terms of the impact parameter instead of the closest approach distance because it is a parameter at infinity. Therefore, for applications in the weak-field approximation the expression (44) is more convenient to the one obtained by Perlick in [28]. Of course, the expression for the deflection angle described in [28] is not limited to the weak field regime, and it must be full considered in the case of strong deflection. In Sec. VI we will derive an expression for the deflection angle with the same characteristics of (44) in the Kerr spacetime and for different electronic density profiles taking into account higher order corrections.

A. Correspondence between homogeneous plasma and massive particles
As is well known there exist a correspondence between the motion of a photon in a cold non magnetized plasma and the motion of a test massive particle in the same background [14]. This correspondence is obtained by identifying the electronic frequency of the plasma ω e with the mass m of the test massive particle and the total energy E = ω ∞ of a photon with the total energy E ∞ = m/(1 − v 2 ) 1/2 . With this identification at hand, and by setting the bookkeeping parameters as ǫ 2 = ǫ 1 = ǫ 3 = 1, we can express the deflection angle for test massive particles as follows, with s = +1 for prograde orbits and s = −1 for retrograde ones. As applications of this expression in the next subsection we will consider the propagation of a test massive particle in a Kerr-Newmann background and in a rotating Teo wormhole.

Kerr-Newman spacetime
Let us consider a Kerr-Newman spacetime given by the following line element at first order in a, M and Q 2 , By restricting our attention to the equatorial plane θ = π/2 we have, Using the expression (45) we obtain the deflection angle for prograde/retrograde orbits of massive particles, which coincides with the one obtained in [45] for uncharged particles.

A rotating Teo wormhole
As a second application, let us consider a rotating Teo wormhole described by the line element, where In order to apply the expression (45) to calculate the deflection angle of massive particles at the equatorial plane in this spacetime, it is enough to consider the line element (51) at first order in the parameters b 0 and a: By using the expression (45) for the deflection angle of prograde/retrograde orbits we finally obtain, which coincides with the one obtained in [46].

V. DEFLECTION ANGLE IN WEAK FIELD APPROXIMATION USING WERNER'S APPROACH: A CASE OF MASSIVE PARTICLES
We show here that an alternative method proposed by Werner to compute the deflection angle of light rays in stationary spacetimes, can also be used to describe the motion of test massive particles. Due to this method needs the non trivial construction of an osculating Riemannian metric even in the situation of a linearized and slowly rotating source is a cumbersome task to derive an equation like (45). For this reason, using this approach we will show its effectiveness by recovering the results previously presented for the Kerr spacetime.

A. Kerr spacetime
Let us consider the Kerr metric in the equatorial plane θ = π/2 given by where In what follows we will use the metric (57) to find out the deflection angle of of massive particles. Towards this purpose let us find the corresponding Finsler-Randers metric for the Kerr metric. It is well known that a Finsler-Randers metric F with manifold M and x ∈ M, X ∈ T x M is characterised by the Hessian given as [10] g ij (x, X) = 1 2 Using the correspondence between the motion of photons in a homogeneous plasma and the motion of test massive particles as discussed in the subsection (IV A) we can study the propagation of test massive particles over the Kerr spacetime.
In particular the particle of mass m is assumed leaving the asymptotic region with a speed v as measured by an asymptotic observer and therefore with an energy In the same way let us assume that the particle has an angular momentum As before, our lensing setup is adapted for prograde orbits. Then, we can consider the massive particle as propagating through an effective medium with refraction index given by, The associated Finsler-Randers type metric reads, As we have said before, according to the Fermat's principle in general relativity test massive particles (due to the correspondence with photons in a homogeneous plasma) follow geodesics curves with respect to the Finsler-Randers type metric F . That is, test massive particles follow geodesics given by the condition, It is worth noting that γ F represent a geodesic of the Kerr-Randers optical metric F . Following Werner's method [10], one can now construct a Riemannian manifold (M,ḡ) which osculates the Randers manifold (M, F ) using the so-called Nazım's method [44]. In particular, one can choose a vector fieldX tangent to the geodesic γ F , withX(γ F ) =ẋ, and the new metric characterised by the Hessianḡ At this point, we recall the striking property according to which the geodesic γ F of the Finsler-Randers manifold is also a geodesic γḡ of (M,ḡ), in other words γ F = γḡ ( for details see [10]). In what follows we are going to use the osculating Riemannian manifold (M,ḡ) to compute the deflection angle of a test massive particle. At leading order, it suffices to use the Born approximation for the path of light ray (r(ϕ) = b/ sin ϕ). Near the light ray, one can choose the vector field as follows As shown in the next subsection for asymptotically flat spacetimes the deflection angle will be given by, where K and dS are the Gaussian and surface element associated with the osculating optical metricḡ.

Gaussian optical curvature and deflection angle
Let us now proceed to use of the osculating Riemannian manifold (D R ,ḡ) in a domain (M,ḡ) with the boundary ∂D R = γḡ ∪ C R . In that case, one can apply the GBT stated as with K being the optical Gaussian curvature and κ = |∇γγ| the geodesic curvature. Furthermore we have used the fact that the sum of exterior jump angles at S and R satisfies θ R + θ S → π and the Euler characteristic χ(D R ) = 1. In this way, letting R → ∞ the GBT reduces to the following form integrated over the optical domain D ∞ outside the light ray γḡ. In addition we have used the fact that since γḡ is geodesic one must have κ(γḡ) = 0. Thus in the limit R → ∞ it is not difficult to show that the geodesic curvature yields (71) While for a constant radial coordinate the optical metric yields Combining these expressions we find In other words, the deflection angle (66) is immediately obtained. Finally we can compute the corresponding metric componentsḡ of the osculating Riemannian geometry as well as the Gaussian optical curvature K. Applying the Hessian (59), along with the vector field Eqn. (66), we find the following expressions g ϕϕ = r 2 v 2 + 2M r − 2M av sin 2 ϕr(2 sin 4 ϕr 2 + 3 cos 2 ϕb 2 ) Indeed one can easily observe that in the limit v = 1, these equations reduces to light deflection case [10]. The Gaussian optical curvature is defined by the following relation Making use of the above expression we find the following result in leading order terms where f (r, ϕ) = 1 r 2 sin 4 ϕ + cos 2 ϕ b 2 7/2 30 cos 4 ϕ sin 8 ϕb 2 r 3 − 6 sin 14 ϕr 5 + 12 cos 2 ϕ sin 10 ϕb 2 r 3 − 48 cos 4 ϕ sin 7 ϕb 3 r 2 − 24 cos 2 ϕ sin 9 ϕb 3 r 2 − 30 cos 6 ϕ sin 4 ϕb 4 r − 27 cos 4 ϕ sin 6 ϕb 4 r − 12 cos 2 ϕ sin 8 ϕb 4 r + 12 cos 6 ϕ sin 3 ϕb 5 + 6 cos 4 ϕ sin 5 ϕb 5 .
Using the equations (79) and (80) and following (68) the deflection angle can be expressed as follows, (81) Finally by evaluating the above integral we find which coincides with (50) by setting Q = 0. In the limit v = 1, we do recover the deflection angle of light by a Kerr black hole [47][48][49]. In this way we show that it is also possible to calculate the deflection angle of massive particles using the Werner's method, or equivalently to study the motion of photons in a homogeneous plasma. In order to study the propagation of photons in a nonhomogeneous plasma we prefer to apply the method discussed in section IV because it is more direct. We are going to study this case in Section VI A. For completeness, and in order to show how Werner's method is a bit more cumbersome to obtain the deflection angle, in Appendix (A) we apply it again to calculate the deflection angle in the rotating Teo spacetime.

VI. HIGHER ORDER CORRECTIONS TO THE DEFLECTION ANGLE
A. Deflection angle of massive particles up to third order in a Kerr background As we have shown, at linear order the expression (45) is a general formula that follows from the use of the Gauss-Bonnet method, and therefore at this order it is not necessary to repeat the computations of the Gaussian and geodesic curvatures for each one of the linearized metrics under study.
However, in order to compute higher order corrections a general formula is not so convenient. Moreover, at higher order, the Born approximation will not be enough, and the integration domain will explicitly depend on the details of the orbit. To our best knowledge, all the existing computations for the deflection angle in stationary spacetimes based in the Gauss-Bonnet theorem have been limited to find linear order corrections due to the intrinsic spin of these metrics (See [13,50,51] to cite some examples). Now we will procedure to fill this gap by computing the deflection angle for test relativistic massive particles in a Kerr background including third order corrections, that is containing terms of higher order as M 3 b 3 , M 2 a b 3 and Ma 2 b 3 . In particular, we will see that for the case of massless particles our general expression reduces to a known formula obtained in the past using other techniques [49,52].
As mentioned above, at higher order we need to go beyond the Born approximation. In particular, in order to compute third order corrections, we need to know the orbit of the massive/massless particle at second order.
Before to do that let us write the general expression for the orbit in the equatorial plane of a general stationary spacetime given by the metric: From the Hamiltonian H = 1 2 (g αβ p α p β + m 2 ), it follows that the orbit equation for the prograde orbit of a particle which is assumed leaving the asymptotic region with a speed v as measured by an asymptotic observer and therefore with an energy and angular momentum given by Eqn.(60) and (61) is given by with ∆ = AD + H 2 and where u γ = 1/r γ . For massless particles (v = 1) it reduces to the expression given in [12]. Let us concentrate on the Kerr spacetime, for other stationary metrics the procedure will be similar. In that case the metric components at the equatorial plane are given by: In this case, the orbit equation reduces to: As we are working in the weak regime, let us consider the bookkeeping parameters As we are interested in an expression for the deflection angle which be correct up to third order in these parameters, that is an expression containing terms of the form γ 3 , γ 2 δ and γδ 2 , we need to know u γ at second order in these parameters, that is, we require an expression for u of the form: u γ = u 0 +u 1 γ+u 2 δ+u 3 γ 2 +u 4 δγ+u 5 δ 2 +O(γ 3 , γδ 2 , δγ 2 , δ 3 ).
(91) By putting this ansatz into Eqn.(89), and by imposing the asymptotic condition that lim ϕ→0 u = 0, results in the following iterative solution for u γ : The fact that u 5 (proportional to a 2 b 2 ) is different from zero even for a flat spacetime (M = 0) should be expected because the coordinates r, ϕ are spheroidal instead or spherical, and therefore at order a 2 in this coordinates the straight geodesic line is not given by r = b/ sin ϕ. For the same reason, the fact that u 2 is zero (proportional to a b ), should also be expected because when M = 0, the correspondent flat metric depends on a 2 , and therefore at linear order in a the coordinates behave as the standard polar coordinates.
From (20) we obtain the associated optical metric: with the refractive index given by Eqn.(62).
Now we proceed to compute the deflection angle using Eqn. (34). Before to do that, we must to define the integration domain. This domain will be bounded from below by the orbit under analysis determined by r γ ≡ 1/u γ and r = ∞ from above. For the angular coordinate ϕ, we have set the source position at ϕ S = 0, and for the position of the observer we must to set ϕ R = π + α (1) , with α (1) the first order correction in the mass to the deflection angle, that is, in our case given by the first term of Eqn. (50). The reason that second order contributions to α are not needed to add to ϕ R can be found in [53]. In fact, in that reference we show that to obtain second order constructions to the deflection angle the approximation ϕ R = π was enough. We will back to this issue later.
At the consider order KdS reads; and therefore the contribution to the deflection angle in terms of u γ reads Let us remark that in the integration of the angular variable, appear some trigonometric functions that must be valued in ϕ R , and therefore they must be re-expanded in terms of γ and δ around π. As anticipated above, it can be checked by direct computation that if instead of taking ϕ R = π + α (1) we take ϕ R = π + α (1) + α (2) , (with α (2) been formed by terms of order δ 2 = M 2 /b 2 and δγ = M a/b 2 )) the contribution to α K of the resulting new terms are of higher than third order. It is easily seen from the fact that the only terms that could potentially contribute at third order to the deflection angle (by introducing second order corrections in ϕ R ) arise from the integration in the angular variable of (97) However, as the integration of (97) in the angular variable will produce a term proportional to with a 1 , a 2 , a 3 constants, it follows that after reexpanding this term in powers of M , only the first term in M of ϕ R (proportional to a 1 ) will contribute at the desired third order. Hence, in order to compute the deflection angle up to third order, we need to know the orbit a second order and the deflection angle at first order. Using the same argument it can be shown that in order to know the deflection angle at n order using (34), we need to know the orbit at n − 1 order and the deflection angle at n − 2.
To finalize, we need also to compute the contribution α κg to the deflection angle. In our casê The associated geodesic curvature is given by On the other hand, dl is obtained from At the considered order we obtain that κ g dl is given by Hence, the contribution α κg to α is obtained from Taking cognizance of (96) and (103); and by using the adimensional parameterâ = a M (which for a black holes must satisfy |â| ≤ 1) the resulting deflection angle reads For massless particles it reduces to: which agrees with known expressions [49,52]. As in the previous cases, for retrograde orbits we must to change the sign of the terms which are linear inâ. Therefore, we have shown that the method proposed by Ono et.al. [12] can be applied in a successfully way to compute the de- flection angle to higher order for stationary spacetimes, and moreover, it can be applied to massive particles. We are not aware of a previous presentation of (104) in the literature.
Note also that using the correspondence between the motion of test massive particles and photons in an homogeneous plasma, (104) also gives the deflection angle of light rays with frequency ω ∞ in a homogeneous plasma characterized by ω e = constant, by simply replacing v in (104) by the group velocity v gr of the light ray in plasma given by v gr = n 0 = 1 − ω 2 e /ω 2 ∞ in units where c = 1. In Fig.(2) and Fig.(3) we consider a black hole with parameters similar to the super-massive black hole in the center of our galaxy Sgr A*, namely we take M = 4.1 × 10 6 M ⊙ ,â = 0.6 and we assume that b = 100M . In both figures, we plot the different contributions α (1) , α (2) and α (3) to the total deflection angle. In (2) we plot the α (i) 's components in terms of the veloc-ity of the massive particle and in Fig.(3) in terms of the quotient ω e /ω ∞ . For both situation we consider prograde and retrograde orbits (s = ±1) and also the nonspinning caseâ = 0.
In Appendix (B) we will present and alternative derivation of (104) using a method originally proposed by Aazami, Keeton and Petters [49].
B. Deflection angle of light rays in a non-homogeneous plasma up to third order in a Kerr background Let us consider a stationary and axysymmetric spacetime surrounded by a cold nonmagnetized plasma and whose restriction to the equatorial plane has the form (83). In such a situation, the orbit equation can be derived to the Hamiltonian (4) resulting in In general for a prograde photon propagating in a plasma environment p t = − ω ∞ and p ϕ = −p t bn 0 , with b the impact parameter and n 0 the asymptotic value of the refraction index. In the next discussion we assume that the electronic number density profile is given by, withÑ 0 a given constant and k a positive integer number. In the case of spherically symmetric spacetimes a detailed discussion of the leading order contribution to the deflection angle and also to associated optical quantities produced by this kind of density profiles was carried out by Bisnovatyi-Kogan and Tsupko in [15]. In particular, this family of profiles can be useful in the study of plasma environments in galaxies and galaxy-clusters. Note also, that for the case of the supermassive black hole in the center of our galaxy Sgr A*, a plasma density profile with a radial dependence of the form r −1.1 has also been considered by different authors [54,55]. Even when it is not exactly in the family described by Eqn.(107), for k = 1 we can have an estimation of the behaviour of the deflection angle for this kind of system. Higher powers of k (including k = 6 and k = 16) also appear in the discussion of light ray propagation in an extended solar corona model [53,[56][57][58][59], In [15] analytical closed expressions for the deflection angle were derived for any k. To higher order, a general expression is more difficult to obtain, and the intermediate steps are not very illustrative. For this reason, in the following we only consider the cases k = {1, 2, 3}, showing in detail the way to proceed for the case k = 2, and only given the final expressions for the other two cases. For any other value of k not considered in this work, the expression for the deflection angle can be found in exactly the same way as in the case k = 2. Note also that our final expressions can not only be used when the rotation of the spacetime is taken into account, they can also be useful for spherically symmetric spacetimes if higher order corrections due to the presence of plasma need to be considered. In the following, we also assume that the frequency ω ∞ of the photon at the asymptotic region is much bigger than the value of plasma electronic frequencyω e ≡ ω e (b) at r = b. That is, we assume that ǫ =ω 2 e ω 2 ∞ ≪ 1. Note that for the present situation, n 0 = 1. Hence, the orbit equation can be rewritten as Let us restrict our attention to the Kerr metric with the functions A, H, B and D given by (85)-(88). Using the ansatz, and by imposing the asymptotic condition that lim ϕ→0 u = 0, result in the following iterative solution for u γ : At the considerer order KdS reads; In order to compute the deflection angle to higher order, we also need to know the expression for the deflection angle at first order. It can be directly obtained from Eqn. (44) with the result: Therefore the contribution α K to the deflection angle including third order terms is determined by 3 The geodesic curvature reads, while dl is given by Hence, Therefore, the contribution α κg to α is Finally, collecting the results of eqs. (121) and (125), and writing them in terms ofâ = a/M we obtain: The first three terms in Eqn.(127) represent the contribution to the deflection angle due to pure plasma and the rest of the terms give the contribution to the deflection angle due to the coupling between plasma and gravitational effects. As before, the deflection angle for retrograde orbits is obtained from (127) by changing the sign in all the terms linear inâ.
The contribution of each component to the deflection angle will depend on the parameters M, a, b and the quotient ǫ between the electronic frequency of the plasma and the photon frequency ω ∞ . Now, we write without derivation expressions for the deflection angle in the cases k = 1 and k = 3 obtained using similar steps. For k = 1, i.e. N =Ñ 0 /r, the final expression for the deflection angle reads: with α p given by Note that for this case, there is not a pure plasma contribution at second order in ǫ. Similarly, for k = 3, i.e. N =Ñ 0 /r 3 , the contribution α p to the deflection angle reads: In all the cases, the linear term in ǫ can be considered as particular cases of the known expression obtained in the past by Giampieri [57] and reobtained in a more general context by Bisnovatyi-Kogan and Tsupko in [15]: with Γ(x) the Gamma function. In Fig.(4) we consider a black hole with the same parameters as in Fig.(2) and Fig.(3) assuming again that b = 100M . We plot the leading order contribution to the deflection angle given by α p for the considered cases in terms of ǫ ∈ [0.01, 0.018]. For this range of values, the deflection angle produced by the plasma is of the same order as the Kerr vacuum angle (but with opposite direction). That is, as is well known, and can also be directly seen from Eqn.(44) the effect of a non uniform plasma that decays with r is to produce a divergent lens effect. Note that for the same ǫ value, the faster the plasma profile decays, the greater the divergence effect on the light rays. Of course, it is a chromatic effect, and different values of ǫ will contribute at different ways. Similarly, in Fig.(5) we plot the contribution to the deflection angle due to the rest of the terms in α vac and α p for the three considered plasma density profiles. The case s = ±1 corresponding to prograde/retrograde orbits respectively. For completeness, we also plot the contributions to the deflection angle for a nonrotating BH (â = 0).

VII. CONCLUSION
In this paper we have used the Gauss-Bonnet theorem to study the gravitational deflection of massive particles and light rays in stationary spacetimes. For that purpose we have obtained a Finsler-Randers type metric which was utilized to calculate the deflection angle of light rays and massive particles. More specifically, using the Finsler-Randers metric for stationary spacetimes along with the correspondence between the motion of light rays in a plasma medium and relativistic test massive particles we have extended the recent method proposed by Werner. As particular examples we have obtained the deflection angle in a rotating black hole and a rotating wormhole geometry. However, the most important result arising from our study, is the fact that we have been able to obtain higher order corrections in the expression for the deflection angle in the Kerr spacetime background applying the Ono et.al. approach. We subsequently analysed also the deflection angle for test mas-sive particles and light rays propagating in a non homogeneous cold plasma by including third order correction terms.
Appendix B: Deflection of light rays up to the third order in a Kerr spacetime using direct integration As an independent check of our expression for the deflection angle given by (104) we will use a direct inte-gration of the orbit equation following a straightforward extension of the method presented by Aazami, Keeton and Petters and Hansen in [49] orignally introduced to study null geodesics.
The orbit equation for both prograde (s = +1) and retrograde orbits (s = −1) in the equatorial plane can be reexpressed as: where Expression (B1) reduces to Eqn.(B14) of [49] in the case of null rays (Γ = 0). The deflection angle is obtained as in the usual way by where r 0 the value of the radial coordinate at the closest approach to the lens. The relation between the impact parameter b and r 0 follows from require dr dϕ r=r0 = 0, which from (B1) implies that The positive real solution of the last equation is By substituting Eqn.(B7) into (B1), and changing the integration variable r in Eqn.(B5) by x = r0 r , the deflection angle is determined by the following expression: with h = M/r 0 . By assuming that h ≪ 1, we can do a Taylor expansion in h of the integrand of (B8), which after integration and conserving terms up to third order in h gives, α = b 0 π+a 1 h+(a 2 +b 2 π)h 2 +(a 3 +b 3 π)h 3 +O(h 4 ); (B9) These expressions generalize the results found in [49] to the case of test massive particles. We can also write the deflection angle in terms of an expansion in γ = M/b using that h can be perturbatively obtained from (B7), In terms of γ the deflection angle reads: A 3 = 1 2G 7/2 6a 1 Γ F 2 + 4F 2 (a 2 + b 2 π) G + Γ(Γ a 1 + 4a 2 + 4b 2 π) + 2(a 3 + b 3 π) G 2 + 5a 1 F 4 (B21) Of course, as F and G are already functions of γ, the A i coefficients must be expanded in γ, and after that, by replacing in Eqn.(B17) we recover the deflection angle as given by Eqn.(104) with the linear terms inâ multiplied by s.