The perturbative approach for the weak deflection angle

Both null and timelike rays experience trajectory bending in a gravitational field. In this work, we systematically develop a perturbative method to compute the deflection angle of rays with general velocity v in arbitrary static and spherically symmetric spacetimes and in equatorial plane of arbitrary static and axisymmetric spacetimes. We show that the expansion in the large closest approach x0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0$$\end{document} limit depends on the asymptotic behavior of the metric functions only, and the generated integrand is always integrable, resulting in a deflection angle in a series form of either x0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0$$\end{document} or b, the impact parameter. Using this method, the deflection angles as series of both x0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0$$\end{document} and b are found in Schwarzschild, Reissner–Nordström and Kerr–Newman spacetimes to 17-th, 15-th and 6-th orders respectively, for both lightrays and particles with general velocity. The effects of the impact parameter, velocity and other parameters of the spacatimes are briefly analyzed. Moreover, we show that for spacetimes whose metric functions are only asymptotically known, the deflection angle in the weak field limit can also be calculated. Furthermore, it is shown that the deflection angle in general static and spherically symmetric spacetime and equatorial plane of static and axisymmetric spacetime to the lowest non-trivial order, depends only on the impact parameter, velocity of the particle, and the effective ADM mass of the spacetime but not on other parameters such as charge or angular momentum. These deflection angles are used in an exact gravitational lensing equation and the corresponding apparent angles of the images of the source are also solved perturbatively.


Introduction
One of the classical and important consequences of General Relativity (GR) is the deflection of lightlike geodesics in curved spacetimes. One hundred years ago, Eddington's a e-mail: junjijia@whu.edu.cn (corresponding author) observation of the star position shift played a major role in helping GR won its acceptance by physicists and the public. Nowadays, the deflection of lightlike rays is usually observed by very long baseline radio interferometer [1] with an accuracy of sub-microarcsecond [2] . In 1991, the deflection of light due to Jupiter was observed [3] and that due to Saturn might be observed by Gaia [4].
The deflection of light leads to one important theoretical and observational tool in modern astronomy and cosmology: the gravitational lensing (GL). After the observation of the first GL in 1979 [5], many features, including luminous arcs [6,7], Einstein cross [8], Einstein rings [9], GL of CMB [10][11][12], and supernovas [13,14] and even composition features such as the SN Refsdal which combines the Einstein cross with the GL of supernova [15] have been observed. The GL were then used to study properties of the lens [16], coevolution of supermassive black holes (BHs) and galaxies [17], cosmological parameters such as large scale structure (for a review see [18]), properties of the supernova [19], dark matter substructure [20,21], and to discriminate alternative gravitational theories. More recently, with the discovery of gravitational wave (GW) [22][23][24][25], observation of its GL effects has also been proposed and put to use [26,27].
Although traditionally lightrays have been the main messenger in the observation of GL, with the observation of extragalactic neutrinos from SN 1987A [28,29] and blazer TXS 0506 + 056 [30,31] and the GWs, it is clear that these two kinds of messengers in principle can also go through the bending process and be observed in lensing scenarios. One of them, neutrinos, are long known to have non-zero small mass [32]; while for GWs, its speed can also deviate from the speed of light in gravitational theories beyond GR [33,34]. Therefore in considering their lensing or GL of any other massive particles, in principle one should compute the deflection angle and time delay formulas using timelike geodesics. Currently, most theoretical works on the GL of these two messengers are still using corresponding values obtained using lightlike rays [27,[35][36][37][38]. Recently, some of us computed the deflection angles in Schwarzschild and Reissner-Nordström (RN) spacetimes and time delay in Schwarzschild spacetime for general velocity [39]. It was shown that the difference of the apparent angles of lightrays and neutrinos can be correlated to the neutrino mass and mass hierarchy [40,41]. Moreover, the dependence of the deflection angle on messenger velocity can also cause small change of the BH shadow size in these spacetimes [40,42]. These results show that in considering effects related to the trajectory bending of massive particles, timelike instead of lightlike geodesics should be used if a high accuracy is desired.
In GLs of both massless and massive particles, the GLs in the strong field limit is important for a few reasons, especially for their applications in strong field test of GR and alternative gravitational theories. In the strong field limit of arbitrary statics and spherically symmetric (SSS) spacetime, Bozza [43] developed a general method to calculate the deflection angle for lightrays and showed that it always diverges logarithmically as ln(x − x 0 ) where x 0 is the closest approach of the ray. Many authors successfully calculated the deflection angles in the strong field limit for particular spacetimes using this method [41,[44][45][46][47]. However, from an observational point of view, all the light bending and GLs seen up to now are in the weak field limit, i.e., with small deflection angles. GL in the strong field limit will not be directly observable in near future experiments because of its high resolution and magnification requirement. For example, to resolve the relativistic images produced by the Milky Way galaxy BH requires an angle resolution of 10 −2 microarcsecond, which is 2 orders smaller than the resolution currently reachable [41].
In this paper therefore, we still concentrate on the deflection angle and GL in the weak field limit, but for general velocity. We propose a general formalism for computing the deflection angle of particles with arbitrary velocity in SSS spacetimes and in the equatorial plane of stationary and axisymmetric (SAS) spacetimes. This procedure allows us to expand the deflection angle in the power of the reciprocal of the closest approach x 0 or impact parameter b, which are large quantities in the weak field limit. Moreover, it will also be shown that this procedure can be used to calculate the deflection angle of metrics that are only known in the asymptotic region. This is particularly useful because in most spacetimes that have a complex matter distribution, the metric functions are usually not analytically solvable while their asymptotical behavior can usually be obtained using series method. We will also show that for all equatorial geodesics in any SAS and asymptotically flat spacetimes (note that SAS spacetimes cover SSS spacetimes), the deflection angle to the lowest non-trivial order always takes the form where m is the ADM mass of the spacetime, b is the impact parameter and v is the speed of the test particles at infinity. In other words, all other spacetime parameters such as effective charges or angular momentum will not influence the deflection angle at this order. Again, we emphasis that all results in this paper apply not only to null rays but signals with general velocity. Previously, the weak deflection angle has been calculated mainly using two slightly different but yet connected approaches. The first and most traditional way is the direct integration method which tackles the integral for the deflection angle directly. The second, which is more recent and also very promising, is to utilize the Gauss-Bonnet theorem to find the deflection in a somewhat more indirect but elegant way [48][49][50][51][52][53]. From this point of view, our work leans more towards the first category.
We arrange the paper in the following way. In Sect. 2, we present the perturbative procedure for computing the deflection angle in SSS spacetimes to any desired order. The results are then used in Sect. 3 to find the deflection angle to the minus 17th order of b in Schwarzschild spacetime and minus 15th order in RN spacetime. We also give an example (the SU(2) Yang-Mills-Einstein solution) for computing the deflection angle in asymptotically known spacetimes. In Sect. 4, the deflection angle in the equatorial plane of Kerr and Kerr-Newman (KN) spacetimes are computed to the minus 6th order of b. We then in Sect. 5 show how all these deflection angles can be used to find the apparent angle in GL perturbatively. Lastly, a few possible applications and directions of extension are discussed in Sect. 6. Throughout the paper we use the geometric unit G = c = 1.

Deflection angle in SSS spacetimes
For general SSS spacetime, the metric can always be written as ds 2 = −A(x)dt 2 + B(x)dx 2 + C(x)(dθ 2 + sin 2 θ dφ 2 ). (2) For this metric, the geodesic equations after two first integrals take the forṁ where˙denotes the derivative with respect to the proper time (or affine parameter) λ. Here we have already set θ(λ) = π/2 without losing any generality and κ = 0, 1 for lightlike and timelike particles respectively. E is the energy of the lightlike ray or that of unit mass of the test particles at infinity, which can be related to their velocity v at infinity through L is the angular momentum of unit mass of the test particle at infinity, satisfying where b is the impact parameter of the ray. Here and henceforth, we assign the impact parameter to be a positive constant, while the angular momentum L in principle can change sign. Using Eqs. (3)-(5), it is easy to show that after some simple algebra, one can integrate dφ/dx to obtain the deflection angle of a ray both originating and propagating from infinite radius as [54] α(x 0 ) = I (x 0 ) − (2n + 1)π (6) where the change of the angular coordinate, I (x 0 ), is Here n is an integer such that −π ≤ α(x 0 ) < π. In the weak field limit, we should always choose n = 0. Setting κ = 0, Eq. (7) reduces to the change of the angular coordinate of lightlike ray given in Ref. [43,54] The closest distance of approach x 0 can be linked with L by settingẋ = 0 in Eq. (3c). Further using Eqs. (4) and (5), one has the following relation between the impact parameter b and x 0 As we stated in Sect. 1, we are interested in the weak field limit of the deflection angle. To calculate α(x 0 ) in weak deflection limit, then we should let x 0 approach infinity. The main content of this paper is then to show that the change of the angular coordinate I (x 0 ) and consequently the deflection angle α(x 0 ) can be expanded as a power series of 1 x 0 and the integral can always be carried out. Moreover, we find a practical procedure to calculate the coefficient of each power, for both null and timelike rays.
To carry out the expansion and integration of Eq. (7), we first change variable from x to u = x 0 /x so that I (x 0 ) becomes where y(x 0 , u) stands for the integrand. This change of variable effectively transforms a partial dependence of I (x 0 ) on x 0 through the lower limit of the integral to a full dependence on x 0 through the integrand y(x 0 , u). For this integrand, since in the weak field limit, x 0 is much larger than any characteristic length intrinsic to the spacetime, we can always do an asymptotic expansion in powers of 1 where denotes the derivative with respect to z, and (15) and y n (x 0 , u) (n ≥ 2) in Eqs. (12) and (13) denote the coefficient of 1 x 0 n . Setting κ = 0 in Eq. (12) yields for lightlike rays, Integrand (12) and (16) have a few remarkable properties.
It is seen that the order 1 term is not shown explicitly here because of its length, but we can show that it only contains terms up to second derivative of the metric functions. In general, it is found that this expansion can always be carried out to arbitrary desired order of 1 x 0 and the coefficient for the 1 x 0 n order involves up to the nth derivative of the metric functions. Note that although the expressions might appear long and tedious as the order increases, the involved mathematics are only taking derivatives and finding limits, and therefore still algebraically simple. In this sense, this expansion is straightforward and tractable, especially when the metric functions are explicitly known. As will be seen in Sect. 3, for explicit spacetime metrics the coefficients for all orders of 1 x 0 n are concise. Therefore, length of the results of I (x 0 ) obtained after integrating u in Eqs. (12) and (16) will also be much shorter comparing to the length of the integrand.
The second property, which is simple but fundamental, is that the expansions use only the asymptotic behavior of the metric functions to the desired order. In other words, from the point of view of the deflection angle α(x 0 ), the rays with asymptotically large x 0 do not experience how the spacetime is curved in the central region. Spacetime in the large x 0 region to various orders of 1 x 0 is enough to determine the deflection angle to the corresponding order. Although intuitively this seems simple, it is this property that allows us to determine the deflection angle in spacetimes whose metric function are not fully but only asymptotically known. We shall see in Sect. 3.3 how the deflection in colored BH spacetime can be calculated.
The third property of these expansions, is that they are always integrable if the metric functions A(x), B(x) and C(x)/x 2 are finite at large x and have asymptotical expansions of the form where a n , b n and c n are finite. These conditions are certainly satisfied by many familiar spacetimes such as Schwarzschild, RN, Hayward [55], Bardeen [56], Gibbons-Maeda-Garfinkle-Horowitz-Strominger (GMGHS) [57][58][59] and Janis et al. [60] metrics. With these conditions, clearly all the terms in various numerators in y n (x 0 , u) in Eq. (12) can be expanded as series of their arguments, which become powers of u. The only non-trivial parts are the denominators in y n (x 0 , u), i.e., the terms 1/ h(z, u) n . In the limit z → ∞, they can be simplified as A further transformation of variable u → cos θ will transform the expansion (13) into the following form Here y n,m are functions obtained after the transformation and they do not depend on the integration variable. Because of the integrability of functions of type cos m θ (1+cos θ) n (m ≤ n), the expansion (12) becomes always integrable under conditions Eq. (17). In short, with the integrand expanded as (12), one can then do the integral of u to get the change of the angular coordinate as a power series of 1 where I n is the n-th order coefficient in I (x 0 ). It is important to note that the integral here can always be carried out because of the integrability of Eq. (19). Because the impact parameter b is much easier to be linked to observable quantities in GLs, it is also desirable to write the expansion (20) in series of 1 b . In that case, we can attempt to directly solve Eq. (10) to obtain the inverse function and then substituting into Eq. (12) and expanding around large b to obtain the desired expansion Indeed, because what is needed is only the expansion of (21) by not necessarily the inverse function l −1 1 b itself, one can use the Lagrange inversion theorem to directly find this expansion without solving l −1 1 b explicitly This will be especially useful for those metric functions in Eq. (10) that are hard to find the inverse function l −1 1 b analytically. Moreover, for ultrarelativistic massive particles, we can also further expand the integrand (12) around v = 1 and integrate to find I (x 0 ) or I (b) as a series expansion of the velocity difference Apparently, the leading terms in the m summation n=0 I n,0 x n 0 or n=0 I n,0 b n will be just the change of the angular coordinate of lightlike rays.

Application to particular SSS spacetimes
In this section, we apply the procedure in Sect. 2 to some known SSS spacetimes, namely the Schwarzschild, RN and the colored BH spacetimes. We will find the deflection angles in the weak deflection limit for general particle velocity, i.e., for both lightlike and timelike particles. The state of art of the deflection angle for Schwarzschild metric is to the sixth order of 1 x 0 and 1 b for lightlike rays [61], and for RN metric to the third order of 1 x 0 and 1 b for lightlike rays [61], and to various orders for some other interesting metrics as well [62][63][64][65][66][67]. In the following, we present the change of the angular coordinate to much higher orders for signals with general velocity. As we explained before, although the idea of this formalism is powerful and clear, the calculations, especially taking high order derivatives, simplifications and integration are tedious. Therefore in most cases we will avoid showing the intermediate steps.

Deflection angle in the weak field limit in Schwarzschild spacetime
For the Schwarzschild spacetime, the metric functions take the form After substituting into Eq. (12), simplifying and integrating over u, we find I (x 0 ) to the order 1 with and the S 5 to S 17 terms are given in Eq. (A1) because of their excessive length.
For lightlike particles, setting v = 1 in Eq. (26) yields the result with with S 5,γ to S 17,γ given in Eq. (A2). For ultrarelativistic particles, their change of the angular coordinates can be obtained by expanding (26) around (28) as series of (1 − v). To the order (1 − v) 1 , this takes the form Higher order terms in m x 0 are given in (A3) and terms of order (1 − v) n (n ≥ 2) can be easily obtained too.
In order to write the expansion (26) in power series of 1 b , we first substitute metric functions (25) and κ = 1 into Eq.
Using Eq. (23), we can find m x 0 in power series of m b as m where the coefficients are and the C S,5 to C S,17 terms are given in Eq. (A4).
Putting Eq. (32) into Eq. (27), one finally obtain the change of the angular coordinate in the power series of 1 b for general velocity with and S 5 to S 17 are given in Eq. (A5). Comparing the powers of 1 b and 1 v in each summand of (34), we can see that where δ = 0, 1 for odd and even n respectively. Therefore in order for the entire expansion to converge, the velocity shall not be indefinitely small. Rather, the range of convergence for v is v ∈ (O m b , 1]. Note that when b is large, this range covers velocities of all interested rays or particles that are (potentially) timelike. These include supernova neutrinos and GWs in some generalized theories of GR.
After substituting v = 1 we obtain where and S 5,γ to S 17,γ are given in Eq. (A5). Again, for the ultrarelativistic particle, the expansion in powers of (1 − v) 1 is The higher order terms in m b are given in Eq. (A7). Corrections in high orders of (1 − v) n (n ≥ 2) can be easily obtained too.

Convergence of the series and effect of velocity
Although the main purpose of this work is not to examine the effect of parameters in each metric, numerical study of the deflection angles such as Eqs. (34), (43) and (88) can not only reveal their validity but also the regions of converge for b and moreover the effect of velocity. To examine the convergence of the deflection angle (34), we plot in Fig. 1a the partial sums for lightrays and the corresponding exact values obtained using numerical method to an accuracy of 10 −15 . We also plot in Fig. 1b the contribution from each order S n m b n for lightrays. It is seen from Fig. 1a that for all fixed b/m, the partial sums converge to the exact deflection angle as the number of terms in the partial sum increases. The minimal convergent impact parameter can reach a sub 10m level. Indeed, as one can see from the inset of Fig. 1a that when b 10m, the deflection angle is already larger than 0.2π . This value certainly exceeds the traditional weak deflection limits. Therefore this nice applicability of the perturbative deflection angle (34) shows clearly the power of the perturbative result, especially when high orders in the expansion can be known.
It is also noticeable from Fig. 1a that as b/m decreases from large values, the low order contributions to Eq. (34) becomes less dominant. As b/m approaches the strong field limit which is b c /m = 3 √ 3 for lightrays, the high order contributions become more and more important and eventually cause the divergence of the total deflection angle. In order to determine the accuracy of the deflection angle calculated to certain order, we plotted the contributions from each order in Fig. 1b and a benchmark line of 1 [as]. Clearly in this plot the slope of contribution S n m b n should be −n. What is important is that for any fixed b, as the order increases, the contribution from each order decreases by a factor smaller than 1. Therefore this decrease guarantees that the series converges at all b considered. In particular, one sees that even for b as small as ∼ 9.7 m, the deflection angle expanded to the 17th order is still accurate to the 1 [as] level, which is roughly the limit of GL observations of galaxies and galaxy clusters.
To study the effect of signal velocity on the deflection angle, we plot in Fig. 2a the contribution from each order to the deflection angle for v = 0.9c, and in Fig. 2b for fixed b = 10m and increasing v from 0.35c to c. Comparing Fig. 2a to Fig. 1b, one can see that as v deviates from c, the contribution to the deflection angle from each order also increases. This leads to an increase of the impact parameter (from 9.4 to 10.5m) at which the accuracy of the expansion can reach the 1 [as] level. For a fixed b = 10m, it is seen in Fig. 2b that as v decreases, the deflection angle increases as dictated by the series (34). However, it is known that for massive particle with velocity v at infinity, there exists a critical impact parameter b c [41] b below which the particle will spiral into the black hole. Using this relation, it can be worked out then when v = v c = 0.43c, b c = 10m. Particles with velocities at or below this value will all have larger b c 's and therefore for these particles if their impact parameter b = 10m, then they will experience a divergent deflection angle and eventually be captured by the BH. This is also seen from the sharp deviation of the perturbative deflection angle (34) from the exact value starting from v/c = 0.5. Below this velocity, the deflection angle, Eq. (34), becomes invalid.

Weak deflection angle in RN spacetime
The RN metric is given by where q is the total charge of the spacetime. Substituting this into Eq. (12), simplifying and integrating over u we get the change of the angular coordinate for RN metric to the fifteenth order of 1 x 0 for general velocity v where whereq ≡ q/m, S 1 -S 4 are the corresponding terms in Schwarzschild spacetime given in Eq. (27) and R 5 -R 15 are given in Eq. (B1). The expansion coefficients (44) show that the effect of the charge in the RN metric starts to appear in the deflection angle from the second order in R 2 . Rewriting Eq. (44c), we have This implies that the total deflection angle decreases monotonically as |q| increases. This agrees with the observation in Ref. [40]. Similar to the Schwarzschild case, comparison of R n+1 and R n in the large x 0 and small v limit also suggests that in order for the expansion (44) to converge, the velocity of the particles should be bounded within v ∈ ( √ m/x 0 , 1]. For lightlike particles, setting v = 1 we obtain its corresponding I RN,γ (x 0 ) where the coefficients are given by where S n,γ are given in Eq. (29) and R 5,γ to R 15,γ are given in Eq. (B2). The result (47) agrees with the deflection angle in Ref. [61]. For ultrarelativistic particles, the change of the angular coordinates can be expanded around I RN,γ (x 0 ). To the order (1 − v) 1 , this is where for simplicity in the curl bracket the orders higher than were not shown. In order to transform the change of angular coordinates (43) into a power series of 1 b , we shall use the relation (10) which in this case becomes for κ = 1 Using this, m x 0 can be expanded as power where the coefficients are with C S,1 to C S,4 present in Eq. (33) and C RN,5 to C RN,15 given in Eq. (B3). Substituting this into expansion (43), the change of the angular coordinate becomes where and the higher order terms are given in Eq. (B4). For lightlike rays, this reduces to where the first three orders were known in Ref. [61]. Again, for the ultrarelativistic particle, the change of the angular coordinate expanded to the (1 − v) 1 order is and terms of order m b 5 and (1 − v) 2 can be similarly computed but not shown for simplicity reason here.

Effect of charge on deflection angle
In n and therefore the total deflection will all decrease as q increases. Figure 3c further illustrates that this effect indeed is persistent to the region thatq > 1, i.e., from an RN BH spacetime to naked singularity case. A comparison of the deflection angle (52) to the exact value also reveals that this expansion approximate the exact deflection angle perfectly in both the RN BH and naked singularity cases. As explained in Sect. 2, the change of the angular coordinate in the large x 0 limit shall be determined by the spacetime in the large x 0 region. In literatures, there are quite some interesting SSS spacetimes whose metric functions are not known explicitly and/or analytically due to the complexity of the Einstein equations. However, very often their asymptotic behavior can be known from series analysis. We take the colored BH in the SU(2) Yang-Mills-Einstein theory as an example [68]. We will show that in this spacetime the change of the angular coordinate can be solved to the order 1 x 0 1 without having to know the explicit form of the metric functions.
The metric of the colored BH are given at large x by [68] A where −1 < δ k < 0, 0 < M k < 1 are constants and k = 1, 2, . . . are indices characterizing colored BHs of different masses. It is clear that asymptotically this metric is almost identical to the Schwarzschild spacetime at the first order except the factor e −2δ k in front of A(x). Substituting the metric into Eq. (20), one finds the change of angular coordinate for general velocity to the order 1 It is seen that for massive particles, I CB (x 0 ) is affected by not only the asymptotic mass M k but also δ k . For lightlike rays, the change of the angular coordinate becomes To the order 1 x 0 , this equals the changes of the angular coordinate (28) in the Schwarzschild spacetime and (52) in the RN spacetime.
The expansion of M k x 0 in term of M k b in this metric can be found using the same procedure as in Eq. (23) to be Substituting this into Eq. (59), one can obtain the change of the angular coordinate to the order of 1 b . Moreover, it is also easy to see that for lightlike rays, this becomes

Weak deflection angle in Equatorial plane of SAS spactimes
To find the deflection angle in the weak field limit in a non-SSS spacetimes described by metric g μν , in principle one could attempt to derive through the geodesic equations the differential equation that the angular coordinate satisfies where f ns denotes some function derivable from the metric functions g μν . Here the coordinate x should resemble the meaning of distance when it is large, E denotes energy per unit mass of the ray at infinity and p collectively stands for all other parameters that might appear in the metric, such as spacetime mass m, angular momentum per unit mass a etc.
Integrating Eq. (63) then yields the change of the angular coordinate for a ray coming and going to spatial infinity In order for the above process to make physical sense and the integral to be eventually carried out, usually a few requirements should be satisfied. The first is that there should exist a spatial infinity (x → ∞) at which the spacetime is approximated by a flat spacetime. Only in this circumstance, the interpretation of the integral (64) subtracting π as the deflection angle makes physical sense. The second is that the metric should still possess certain amount of symmetries. Indeed, the geodesic equations are originally second order differential equations while Eq. (63) is already first order. For this to happen, the spacetime will have to have at least the necessary symmetries to allow the first integrals. Moreover, only when the spacetime is static, the function f ns (x, E, p) in Eq. (63) and so that the change of the angular coordinate can be static. Finally, from the practical point of view, the function f ns (x, E, p) should still be simple enough for the integration to be carried out. This can require one or more of the following conditions or techniques, such as simple enough metric functions, specially picked trajectories in the spacetime, or series expansions of the integrand over certain parameters.
The above considerations naturally singles out the SAS spacetimes as next suitable candidates for computation. We start from the most general SAS metric given by the Weyl-Lewis-Papapetrou line element [69][70][71] where (t, ρ, φ, z) are the coordinates and U, ω and γ are arbitrary functions of ρ and z only. For latter easier reference, we can rewrite this metric into the form where we have set . Note there is a relation between these functions We further assume that this spacetime allows a plenary motion for particles in a plane with fixed z, which will be called the equatorial plane henceforth. Indeed, without losing any generality, the coordinates can be shifted along the z axis so that this plane becomes z = 0. The motion in the equatorial plane then becomes effectively 2+1 dimensional whose metric after suppressing the z coordinate in the metric functions becomes We then concentrate on the motion of particles in this plane. The geodesic equations corresponding to metric (68) can be readily obtained aṡ where we have substituted (67). Here E and L are still energy and angular momentum per unit mass of the particle at infinite ρ. Note that because the SAS spacetimes might carry a nonzero angular momentum, the direction of the angular momentum L is important in determining the shape and deflection angle of the geodesics. Using Eqs. (70) and (71) to find dφ/dρ, it is easy to show that the deflection angle for rays both coming and going to infinite ρ takes the form where we add an absolute sign because I (ρ 0 ), the change of the angular coordinate defined as can be close to π or −π when L is positive (s = −1) or negative (s = 1) respectively. The angular momentum L can be linked to ρ 0 by usingρ ρ=ρ 0 = 0 in Eq. (71) to find Here since in the small spacetime spin limit the first term in the numerator dominates the second term, it is clear then the s = +1 and −1 correspond to the cases that the L are negatively or positively oriented, respectively. Equation (73) usually does not permit integration into closed form in terms of elementary functions. Similar to the procedure in Sect. 2 which leads I (x 0 ) to power series of 1 x 0 and then further to series of 1 b , one can also expand I (ρ 0 ) in Eq. (73) in the weak field, i.e., the large ρ 0 limit, and integrate to find the change of the angular coordinate in power series of 1 ρ 0 and eventually in powers of 1 b . We will not carry out these formal steps as in Sect. 2 here because they are exactly the same way and too tedious to show here. Rather we will only illustrate these steps in subsection (4.1) for equatorial motion in Kerr and KN metrics respectively.

Weak deflection angle in equatorial plane of the Kerr and KN spacetimes
With the above consideration, the next spacetimes we will consider is naturally the Kerr and KN spacetimes. We will show that the formalism for expanding the integrand in series of 1 x 0 can also be used here to find the change of the angular coordinate of both massive and massless particles to high orders. We will only illustrate the procedure and results in the KN spacetime because setting its charge to zero will yield the corresponding results in the Kerr metric. To our knowledge, the deflection angle in Kerr spacetime has been known to fourth order in 1 x 0 and 1 b for lightlike rays [72,73]. The KN spacetime has a metric where functions and are (x, θ) = x 2 + a 2 cos 2 θ, and m, a and q are respectively the mass, angular momentum per unit mass and charge of the KN spacetime.
In the equatorial plane, θ = π 2 and the metric (75) reduces to where The coordinate ρ in the standard metric (68) is related to the coordinate x in metric (77) by [71,74] ρ(x) = x 2 − 2mx + a 2 + q 2 sin θ where θ = π 2 was used. The metric functions in (77) and (68) are then related by Substituting Eqs. (78)-(80) into (73), the change of angular coordinate can be expressed as We can get rid of the angular momentum L in this equation by using Eq. (74) which in this case becomes Again, s = +1 and s = −1 correspond to the cases that L is negative and positive respectively. The integral (81) can not be carried out to find a closed form and therefore an expansion of the integrand for large x 0 is needed. After integrating this expansion, we obtain a series approximation of the change of the angular coordinate.
Here we skip these tedious middle steps and present the result directly to the order m where where R 0 to R 6 are given in Eq. (44).
For lightlike rays, setting v = 1 we easily get their change of the angular coordinate as And for ultra-relativistic particles, their change of the angular coordinate deviates from the lightlike rays' amount. At order (1 − v) 1 this deviation is given by In order to express the change of the angular coordinate in terms of the expansion of the impact parameter, again we use the relation (9) with L given by Eq. (82) to obtain a series expansion of 1 It is seen that when b > 0 is fixed and q = 0, in the retrograde case (sa > 0) the x 0 will be smaller than in the prograde case (sa < 0). Substituting into Eq. (83) we obtain where where R 0 to R 6 are given in Eq. (53). From Eq. (89) it is seen that the effect of the angular momentum per unit mass a only starts to appear from order 1 b 2 . Writing Eq. (89c) out explicitly, we have It is seen that when the geodesic motion is prograde, i.e., sa < 0 the deflection angle decreases as |a| increases. On the contrary, if the geodesic motion is retrograde, i.e., sa > 0, then the deflection angle is larger than the prograde case with same |a|. Note that when Q = 0 for Kerr spacetime, the first three orders agree with Eq. (96) of Ref. [75]. For lightlike rays, setting v = 1, this becomes where The first four orders here agrees with Refs. [72,73].

Effect of angular momentum on deflection angle
In Fig. 4a, we plot each order in the deflection angle (88) for v = 1 andâ = 0,â = 0.3 for both the retrograde and prograde motions. It is seen for the first nontrivial order of the deflection angles, the lines completely overlap for differentâ. This is indeed a consequence of Eq. (88), where the effect ofâ only starts to appear from the third order (the second nontrivial order). For third or higher orders, it is seen that the retrograde motion has a larger deflection angle than in the Schwarzschild spacetime which is yet larger than the prograde motion case. This is indeed a manifestation of the frame-dragging effect of a rotating spacetime. In Fig. 4b, the impact parameter is fixed at b = 10m and the velocity varies from 0.7c to 1.0c with two spacetime rotation directions witĥ a = 0.3. It is seen that as v decreases, the deflection angles at all orders increase for both rotation directions. This can be understand from the observation that slower particles (keeping b fixed) tend to pass by the BH with a closer distance and therefore are more influenced by it, leading to a larger deflection angle. Moreover, the deflection angle is increased regardless whether the BH is rotating or not, or the direction of the rotation. In Fig. 4c, the dependence of the deflection angles on the angular momentum is shown. It is seen that asâ increases, i.e., the rotation of the spacetime becomes faster, the deflection angle increases for retrograde motion and decreases for prograde motion. Moreover, similar to the case of RN spacetime, it is seen that the total deflection angle computed using Eq. (88) works even when the angular parameter is beyond its critical valueâ = 1 −q 2 for the extreme KN spacetime. In other words, the deflection angle (88) is valid for both the KN BH and naked singularity spacetimes. show that this feature is universal to at least geodesics on the equatorial plane of general SAS spacetimes which are also asymptotically flat. These spacetime certainly include all the spacetimes in previous sections, in particular the SSS spacetimes specified by metric (2). This value although sounds intuitively trivial, was never proven rigorously before. Furthermore and more importantly, we will also show that for geodesics in the equatorial plane of such spacetimes, to the next (the lowest nontrivial) order, the change of the angular coordinate always takes the form where m is the ADM mass of the spacetime. Recently, Ref.
[76] derived a deflection angle to the first nontrivial order in terms of an energy-momentum distribution of the SSS spacetime. From that point of view, our result here further proves that to the first non-trivial order, ADM mass is the only parameter of an asymptotically flat spacetime that will affect the change of the angular coordinate. Other parameters in these spacetimes such as charge or angular momentum can only appear at higher orders. In addition, Eq. (93) also dictates how the particle velocity will affect this angle. Expanding Eq. (93) around v = 1 yield the correction of the deflection angle of ultrarelativistic particle with respect to lightlike rays to the lowest nontrivial order of both 1 b and  (1 − v).
We now show (93). In Eq. (73) we have calculated the change of the angular coordinate in terms of general metric functions in Eq. (68). In order for this spacetime to be asymptotically flat, these metric functions should asymptotically satisfy the following conditions [77] where m is the effective total ADM mass of the spacetime. The effective parameters of charge and angular momentum of the spacetime will only appear in the higher orders of this expansion. Substituting Eq. (95) into Eq. (73), changing the integration variables from ρ to u = ρ 0 /ρ, series expanding the integrand in the large ρ 0 limit and then carrying out the integration, the result of the change of angular coordinate is found to be To express Eq. (96) in terms of the impact parameter b, we can again using Eq. (9) and (74) to find to the lowest order Substituting this into Eq. (96), we immediately have to the first non-trivial order which is the desired Eq. (93).

Gravitational lensing in the weak field limit
One main application of the deflection angles is in the GL in the weak field limit. The configuration of the GL in the weak field limit in the SSS or equatorial plane of SAS spacetime is illustrated in Fig. 5. We will denote the distance between the lens and source as d ls and observer and lens as d ol , and the angular position of the source and its image against the observer-lens axis as β and θ respectively. These angles as well as the deflection angle α are all small in the weak field limit. A lens equation linking these lengths and angles can be established, from which the apparent angle θ of the image are usually solved.
In the following, we will start from an exact lens equation derived directly from the geometric relations in GL but without using small angle or large length approximation [78] The reason that we use the exact length equation but not the usual first order equation is that we would like to compute the apparent angle θ using perturbative method to higher orders and therefore this requires a lens equation accurate to high orders too. Equation (100) although can be derived using different approaches (see Ref. [78] for a review and comparison of different approaches) and were used commonly in most GL computations, more exact equation such as Eq. (99) has to be used when the error induced by the Eq. (100) exceeds observational sensitivity. Using Eq. (99), we will then illustrate how the various orders in the deflection angle α (e.g., Eqs. (34), (52) and (88) after subtracting π ) will perturbatively determine the apparent angle θ to the corresponding order. Using this perturbative method, the apparent angle θ to any desired accuracy can be achieved, provided that the deflection angle α can be calculated in prior to that order. Note that our method and result here will apply to not only null signal but signal with general velocity.
For the solution process, we first substitute the deflection angle α whose general form is a power series of 1  After expanding sin β, the first equation (104a) coincides with the leading order lens equation (100). Note that in this equation system, for each i, only θ i is the unknown in the ith equation. θ j with j < i in the ith equation can always be solved in the equations prior to the ith equation. Moreover, except Eq. (104a) which is not a linear equation of θ 1 but still solvable, all other equations are linear in its unknown θ i . Therefore the perturbative method guarantees that the equation system is iteratively solvable so that the final apparent angle can be obtained. This solvability is one of the advantages of the method.
Solving the system (104), we find the solution to θ i as The higher order equations and solutions can be similarly obtained. It is seen from Eq. (105) that the solution to θ m (m ≥ 2) depends on the lowest order apparent angle θ 1 and the deflection angle from the 1st to the mth orders, α n (1 ≤ n ≤ m), as a power law function, except a common denominator which only depend on α 1 and θ 1 but not higher order ones. Moreover, counting the orders of distances d ol and d ls and the small angles θ 1 , one finds that θ m+1 is smaller than θ m by an order of d ol θ 1 . This implies that if a certain calculation accuracy of the apparent angle θ is desired, one can work out from Eq. (105) to what order the deflection angle should be used. This tractability is another advantage of the perturbative method. With the series expansion of θ known, we can find the magnification of the images using and expand around small angles of θ . Again, the result will be of a series form. For simplicity reason, we will not show these results explicitly.

Discussions
We have studied a perturbative method to compute the deflection angle in general SSS and equatorial plane of SAS spacetimes for arbitrary velocity v. It was shown that the involved integral in this method can always be carried out and therefore a series in either the closest radial coordinate x 0 or the impact parameter b can always be found where we used p to collectively denote any parameter of the spacetime, and C n (v, p) and C n (v, p) are two sets of coefficient functions. Using this method, we were able to compute the deflection angle to the 17th order in the Schwarzschild spacetime, 15th order in the RN spacetime and 6th order in the KN and consequently Kerr spacetime. Using these results, we studied how v, b and various parameters of the spacetime affect the total deflection angle and the deflection angle at each order. Two general features are particularly worth mentioning. The first is that although the deflection angles are obtained perturbatively as a series of 1/b when its large, the valid range of the found deflection angle can extend to much smaller b, even when the deflection angle is not small anymore. The second is that the found deflection angles in the RN and KN metrics describe accurately the deflection angle not only for their BH spacetimes but also for their naked singularity cases. This last point here lays the foundation to apply these deflection angle results in the corresponding GLs to reveal the relevant features, if any, of these naked singularity spacetimes. Using this perturbative method, we have shown that the deflection angle in the weak field limit in the asymptotically flat spacetimes depends only on the asymptotical behavior of the metric functions but not their values at small b. In particular, it was shown that the deflection angle of particles with general velocity in an EYM spacetime whose metric is only asymptotically known, can be computed. Moreover, we also illustrated that for equatorial motion in general SAS (including SSS) spacetimes, the deflection angle to the first nontrivial order depends only on the ADM mass of the spacetime and the asymptotical velocity of the particle in the specific way given by Eq. (98).
Regarding the extension of this method, the first and most apparent is to apply this method in other interesting SSS and SAS spacetimes to compute their deflection angles and study the GL effects in the weak field limit [79]. Results found for these spacetimes are expected to not only reveal effect of spacetime parameters to signal deflection or GL, but also properties of the messenger itself, such as the neutrino mass/mass hierarchy and massive GWs. Secondly, we also expect that this perturbative method is applicable to other computations involving integration of the geodesic equations, e.g., in the computation of time delay in GLs. We are currently working along this direction.
Acknowledgements This work is supported in part by the NNSF China 11504276 and MOST China 2014GB109004. The author greatly appreciate discussion with Dr. Nan Yang, Haotian Liu and Ke Huang.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and therefore has no data to deposit.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: High order terms of the deflection angles in Schwarzschild spacetime
For Schwarzschild metric, the S 5 to S 17 in Eq. (26) are and for lightlike rays these become The full form of the expansion (30) is For Schwarzschild metric, when x 0 is expressed in terms of b in Eq. (32), the high order terms are (A4m) The S 5 to S 17 in Eq. (34) are given by