Perturbative deflection angle, gravitational lensing in the strong field limit and the black hole shadow

A perturbative method to compute the deflection angle of both timelike and null rays in arbitrary static and spherically symmetric spacetimes in the strong field limit is proposed. The result takes a quasi-series form of $(1-b_c/b)$ where $b$ is the impact parameter and $b_c$ is its critical value, with coefficients of the series explicitly given. This result also naturally takes into account the finite distance effect of both the source and detector, and allows to solve the apparent angles of the relativistic images in a more precise way. From this, the BH angular shadow size is expressed as a simple formula containing metric functions and particle/photon sphere radius. The magnification of the relativistic images were shown to diverge at different values of the source-detector angular coordinate difference, depending on the relation between the source and detector distance from the lens. To verify all these results, we then applied them to the Hayward BH spacetime, concentrating on the effects of its charge parameter $l$ and the asymptotic velocity $v$ of the signal. The BH shadow size were found to decrease slightly as $l$ increase to its critical value, and increase as $v$ decreases from light speed. For the deflection angle and the magnification of the images however, both the increase of $l$ and decrease of $v$ will increase their values.


I. INTRODUCTION
Geodesic behavior of both null and timelike rays in the strong field limit (SFL) near the ultra-compact objects, such as black holes (BHs) and neutron stars (NS), is important for many phenomena. The accretion of materials [1], the shape of BH shadow [2], the supernova explosion details [3] and the ring-down phase of BH/NS binary mergers [4] are among the many such phenomena. In order to understand these well, to test BH and NS physics and to test GR in the SFL, therefore it is necessary to understand the general geodesic behavior of null/timelike rays near the gravitational center.
Previously, the null geodesics in this limit has been well studied and understood in many spacetimes. The existence of photon sphere (surface) around some BHs [5] has been well known, and the location of the innermost stable circular orbits [6,7] has been intensively investigated. In particular, the deflection of null rays in the SFL has been systematically studied by Bozza et al. for static and spherically symmetric (SSS) spacetimes and equatorial motion in stationary and axisymmetric (SAS) spacetimes [8][9][10] and has since been used to study the deflection and GL in the SFL for many interesting spacetimes [11][12][13][14][15][16][17][18][19][20]. Note that most of the deflection angles found in these works assumed that the source and detectors are located at infinite distance. Through these studies, it is found that when the impact parameter b approaches the critical value b c , the deflection angle di- * Electronic address: junjijia@whu.edu.cn verges as whereā,b are some coefficients. With this deflection angle, then the single most remarkable consequence is the existence of a series of relativistic image with infinite magnification for rays near the critical photon sphere.
The messenger for these GL images and for GL in the weak field limit were traditionally light rays. However, with the discovery of extragalactic neutrinos from SN1987A [21,22] and blazar TXS 0506+056 [23,24] and the gravitational waves (GWs) [25][26][27][28][29], and especially the lensed supernovas [30,31] and simultaneous observation of GW+GRB events [28,29], we know that timelike rays might also be messengers for GL or other gravitational effects in the SFL. Therefore in this work we will investigate the influence of the asymptotic velocity on the deflection angle and GL in the SFL for general SSS spacetimes. Moreover, we adapt a perturbative method developed previously for the calculation of the deflection angle and time delay in the weak field limit [32][33][34][35][36] to the SFL case, and use it to find a series form of the deflection angle for source and detector at finite distances. The leading two orders of this series will produce the deflection angle (1) when the asymptotic velocity is set to light speed and source/detector distances set to infinity. Therefore results in this work will be more general in that it includes the effect of timelike velocity, the finite distance of the source and detector and contributions from orders higher than constant.
The paper is organized as the following. In Sec. II we describe in detail this perturbative method. The result of the deflection angle is given in Eq. (23)- (26) and its realizations in some simple spacetimes are compared with known results. In Sec. III, we study the GL, in-cluding the apparent angles and the magnification, and BH shadow size in the SFL for signals with arbitrary velocity and source/detector from finite distances. In Sec. IV, the perturbative method and the GL results then are applied to the Hayward spacetime and the effects of the charge parameter and signal velocity are analyzed. We summarize and discuss possible extensions in the end.

II. PERTURBATIVE METHOD
We start from the most general SSS metric, which is described by where (t, r, θ, ϕ) are the coordinates and A, B, C are metric functions depending on r only. The function C(r) can be further chosen to r 2 (at least locally). However, sometimes in special coordinate systems this is not explicitly set so we will keep it free here.
The corresponding geodesic equations arė where without losing any generality we have set θ = π/2. κ = 0, 1 for massless and massive particles respectively.
Here L and E are the angular momentum and energy of the signal (per unit mass) respectively. They can be related to the impact parameter b of the signal and its velocity v at infinity by The change of the angular coordinate from source located at radius r s to detector at radius r d then is We assume that the detector is located far away from the gravitational center while r s is not necessarily so. Here r 0 is the closest approach defined byṙ| r=r0 = 0 which after using Eq. (3a) becomes This r 0 can also be related to the impact parameter b by using Eq. (4) and (6) We also assume that the spacetime contains a particle (or photon) sphere below which the particle (or photon) will spiral into the BH and not escape. This particle sphere has a radius r c defined as the critical point of the denominator of Eq. (5) If the closest approach r 0 equals particle sphere radius r c , then we call the impact parameter corresponding to this r 0 the critical impact parameter b c . That is, using Eq. (7), For example, for Schwarzschild spacetime, it is well known that for light signal with κ = 0, the photon sphere has a radius r c = 3M as predicted by Eq. (8) and a critical impact parameter b c = 3 √ 3M which can be verified using Eq. (9). While for light signal in RN BH spacetime, it is known that [37] We now define a function p(x) that is inspired by Eq.
and denote its inverse function as q(x). Clearly, these functions satisfy p 1 r0 = 1 bc − 1 b and Using function q(x), for ∆φ in Eq. (5) we then do a change of integration variables from r to ξ linked by such that after using Eqs. (4), (12) and (13) the various terms in the integral (5) becomes respectively, where q = q(ξ/b c ) and q is its derivative. For quantity η s,d , it is through them that ∆φ depends on the finite distance r s and r d . Also note that arcsin(1 − η d ) is the apparent angle of the GL images observed by the detector (see Ref. [38,39]) and in the infinite r s,d limit, we will have η s,d = 1 because C(r s,d → ∞) → r 2 s,d . The change of variable (15) is the key step that enables the later integrability of ∆φ in subsection II A. Note that the dependence of ∆φ on the finite distance r s and r d are through the η s,d . Organizing the terms in Eq. (15) together, Eq. (5) is transformed to where a ≡ 1 − b c /b. Then our target for calculating ∆φ in the SFL becomes to calculate this integral in the limit b → b + c , i.e., a → 0 + .

A. Perturbative Computation of ∆φ
To proceed from Eq. (16), we will denote the integrand of Eq. (16), except the factor 1/ √ ξ − a, as Since we are interested in the small a limit of the integral, to which the contribution from y(ξ, a) when ξ is small is expected to be large, we can directly expand y(ξ, a) for small ξ and then attempt to carry out the integral. The square-root factor of y(ξ, a) has the following expansion Denoting the rest factors of y(ξ) as f (ξ), then one can show that it has an expansion where the initial index is -1 due to the q factor in y(ξ) and f n are the coefficients. We can then combine the above two expansions and collect according to the power of ξ n where y n,m denote coefficients of corresponding powers of ξ, a and (2 − a). Note that these y n,m are neither ξ nor a dependant and can be directly worked out once the metric functions are known, as will be shown in subsection II A. The ∆φ in Eq. (16) consequently becomes At this point, the integrability of ∆φ becomes apparent because the integrals of the form It is seen that the only divergence in the SFL is proportional to ln a and caused by the k = m = 0 term. In the small a limit, all functions of a here can be further expanded and then this ∆φ can be thought as a quasi-series of a, with only one term in the coefficient of each a n (n = 0, 1, · · · ) containing one ln a. That is, where −C n and D n are coefficients depending on the spacetime behaviour near the particle sphere through y k,m , on the finite distance r s,d through η s,d and on the particle velocity v through both. To the O(a) 0 order, only the m = 0, j = k terms in Eq. (23) contribute to ∆φ, i.e., where C 0 and D 0 (η s , η d ) can be read off directly from Eq. (25) to be To link y n,m to the behavior of the spacetime in a more transparent way, we will first express them using expansions of the metric functions around the particle sphere. It will be seen that for any particular n, the y n,m are linked to the metric expansion coefficients up to the n-th order. And then in Table I, we will apply these results to particular spacetimes and compare ∆φ found using our method to known analytic formulas and numerical values.
Assuming that the metric functions A(r), B(r) and C(r) have the following series expansion at the particle sphere r = r c then we see that the very fact that r c is the radius of the particle sphere demands through Eq. (8) that This condition provides a constraint from which we can express a 1 in terms of a 0 , c 0 and c 1 Moreover, the critical impact parameter in Eq. (9) can also be related to these coefficients It is also interesting to note that this value depends only on the metric functions A(r) and C(r) but not B(r). Furthermore, without losing any generality, C(r) can be set to r 2 . In this case, the expansion (29c) of C(r) only contains three terms If the metric function B(r) = 1/A(r) as in many SSS cases, then the entire set of B(r)'s coefficients b n (n = 0, 1, · · · ) can be correlated with that of A(r) and r c . In these cases, the number of independent coefficients will be much lesser.
Using the expansions (29), together with Eqs. (12), (17), (21) and the Lagrangian inverse theorem, we can compute the coefficients y n,m to any desired order. The first few y n,0 which are crucial for ∆φ to the O(a) 0 order are where constants T 2 , T 3 are It is seen that y −1,0 , which determines the divergence as a → 0 + , is fixed by the metric expansion coefficients a 0,2 , b 0 and c 0,1,2 but not the higher order coefficients. Higher order y n,0 (n = 1, 2) are given in Appendix A in Eq. (34c-d) because of their excessive length. From these, one finds that for a general n, y n,0 are determined by coefficients up to a n+3 , b n+1 and c n+3 . Once the spacetime metric is known and its r c radius is found, then all these y n,0 can be readily obtained from Eq. (34).

B. Comparison with known results
Now we can apply the above Eqs. (25)-(28) for ∆φ to particular spacetimes to verify the validity of our methodology by comparing our results to others in the SFL. In Table I, we list the spacetimes and the quantities we compared, the method for comparison (analytical or numerical) and the corresponding references. All these literatures only worked for null signal and most of the them only presented deflection angle numerically. In all literature compared, our Eq. (25) were able to reproduce the corresponding quantities in that spacetime.

III. GL AND BH SHADOW SIZE IN THE SFL
GL in the SFL (and also the weak field limit), including the image apparent angular positions and their magnifications, are always solved from the GL equation that is applicable to the case. Previously, GL equation was usually built using some approximate geometrical relations between the source and lens distances, the source's angular positions, and the deflection angle for source and detector located at infinity. However, in simple spacetimes such as SSS spacetimes, using ∆φ that works for source and detector located at finite distances, we show now that an easier and more exact GL equation system can be established. This system contains two equations and is applicable to GL in both strong and weak field limits. It also involves very little geometrical approximation and is usually easy to solve for the image apparent angles.
The first equation links ∆φ to the source coordinate (r s , φ s ) and the detector coordinates (r d , φ d ) (see Fig.  1). Note that φ s , φ d ∈ [0, 2π) and without losing any generality, later on we can set φ d = 0. This equation is nothing but the very definition of ∆φ (36) where integer n ≥ 1 is the number that the trajectory loops around the center and p is used to denote collectively all other parameters of the spacetime. The s = ±1 is a sign corresponding to the anti-clockwise and clockwise looping directions respectively. From this equation, one is able to solve for each n the two b's from the two directions that allow the trajectory starting from the source to be received at the detector. Using the ∆φ to the O(a) 0 order, which are given in Eqs. (26), Eq. (36) becomes where we denoted ∆φ sd = φ s − φ d . Immediately the impact parameter can be solved as where b c is given in Eq. (9). Clearly, as n varies, there will be a discrete series of b allowing the signal to reach the detector. Moreover, because the ∆φ(r s , r d , b, p) we substituted into Eq. (36) become more accurate when b approaches b c , i.e., a → 0 + , the solution (38) will be more precise when n is larger. Note that b n,s are also dependant on the finite distances of source and detector through the coefficient D 0 (η s , η d ) with η s,d given in Eq. (15b), although this dependence is weak especially when r s , r d and n are large. The purpose of the GL is usually to solve the image apparent angles θ against the detector-lens axis (see Fig.  1) and their magnification. For geodesics in SSS spacetimes, this can be achieved by using an exact apparent angle formula once b is known  (s = ±1, n = 1, 2, 3, · · · ).
We refer the reader to Ref. [38,39] for the derivation of this formula. Substituting Eq. (38) into the above equation, it is immediately clear that one discrete series of θ n,s (n = 1, 2, 3, · · · ) is allowed from each side (s = ±1) of the lens These are the apparent angular locations of the relativistic images in the SFL. Here in Eq. (40) for the first time we give an explicit formula for the apparent angle that works for all SSS spacetimes, particle velocity, source and detector distance. An approximate form of this result for null rays was obtained by Bozza in Ref. [8] and used/re-derived by many others in particular spacetimes [11][12][13][14][15][16][17][18][19][20]37]. When n → ∞, θ ∞,s becomes independent of C 0 , D 0 as well as ∆φ sd and reaches its minimal value, i.e., the size of the black hole shadow .
It is seen that the shadow size is determined completely by the location of the detector r d , the particle sphere size r c , the metric functions at these locations and the velocity of the signals forming the shadow. Although Eq. (41) is not a difficult result, to our best knowledge, such simple relation between the spacetime metric and shadow size was never reported before, certainly not for signals with general asymptotic velocities.
With the apparent angles of the images known, now we show that the magnifications of the relativistic images form a geometrical series with a convergence factor determined by C 0 only. The magnification of the n-th relativistic image is defined as [37,44] µ n,s = θ n,s β dθ n,s dβ , where β is the angle of the source against the lensdetector axis if there was no GL (see Fig. 1). To facilitate the computation of µ n,s , we have to link β to other known quantities. This can be done through the geometrical relation in triangle SAD that [r d − r s cos ∆φ sd ] tan β = r s sin ∆φ sd , which links β to ∆φ sd . Therefore using Eq. (43) in Eq.
This magnification depends on n and ∆φ sd in some in-teresting ways. Clearly, as n increases, the magnifica-tions forms a geometric series with the convergence factor e − 2π C 0 . This means that the relativistic images will become weaker and weaker as one goes to inner images with smaller θ n,s . On the other head, if one holds n and varies ∆φ sd , then the denominator of Eq. (44) indicates that the magnification diverges at (1) ∆φ sd = 0, π, arccos r s r d , 2π − arccos r s r d if r s < r d ; (2) ∆φ sd = π if r s = r d ; (3) ∆φ sd = 0, π if r s > r d .
(see Fig. 5). Moreover, because 0 ≤ ∆φ sd < 2π, one can also verify from Eq. (44) that the magnification of the n-th relativistic image with ∆φ sd that is not close to the detector-lens axis is always larger than that of the (n + 1)-th image with any ∆φ sd which might be close to 0 or π. This implies that when discussing series of relativstic images in GL in the SFL, the angular location of the source is not limited to ∆φ ≈ π as in the case of GL in the weak field limit or ∆φ ≈ 0 as in the case of retrolensing.
Since usually the relativistic images with n ≥ 2 are not resolvable from each other in reasonable future (see Sec. IV), the images with n ≥ 2 will only form one image. The total magnification of this image is To deduce properties of the BH spacetime, which are characterized by C 0 and D 0 , we can use the flux ratio of this images against the first relativistic image [8], which is Moreover, one can also use the positions of θ 1,s and θ ∞,s to find two ratios From (47) and (48), one is able to solve the coefficients C 0 and D 0 , and the angular coordinate difference ∆φ sd D 0 = π ln (r µ + 1) 3 r θ,+ r θ,− ln (r µ + 1) , ∆φ sd = π ln [(r µ + 1) r θ,+ /r θ,− ] ln (r µ + 1) . Again, as pointed out by Bozza [8], by measuring the flux and angular position ratios, one is able to reconstruct the full SFL expansion of the deflection angle. Further more, since C 0 and D 0 depend on the nature of the BH, this also allow us to deduce information about the central BH.

IV. THE HAYWARD BH SPACETIME CASE
In this section, we apply our perturbative method for solving ∆φ and the GL equations in the SFL to the Hayward BH spacetime, whose metric is given by [45] A(r) = 1 − 2mr 2 r 3 + 2l 2 m , B(r) = 1 A(r) , C(r) = r 2 , (52) where m is the spacetime mass and the charge parameter l is smaller than a critical value l < 4m/(3 √ 3) ≡ l c in order for the spacetime to contain a BH. The equation determining the location of the particle sphere is given by Eq. (8), which for the Hayward BH spacetime becomes After solving this equation numerically for a given E (or v), m and l then using Eq. (9), one can obtain the corresponding critical impact parameter. In Fig. 2, we plot the solved r c and b c as a function of l from 0 to l c for two different velocities. It is seen that the with the increase of l, both r c and b c decrease monotonically from their Schwarzschild values, i.e. 3m and 3 √ 3m. Then one would expect that the gravitational effect become stronger if the test particles come close to the same distance to the critical sphere. The decrease of the signal velocity increases b c more dramatically than r c , which is also expected because it is known that in Schwarszchild spacetime as v decreases to zero, b c approaches infinity while r c only increases to 4m. Similar increase of b c and r c also happens in the RN BH spacetime [37].
Once r c is known, then the expansions of the metric functions in (52) around r c can be worked out. To the O(r − r c ) 4 , they are Reading off the corresponding coefficients a n , b n and c n and substituting into Eqs. (34), one can obtain the coefficients in the deflection angle of the Hayward spacetime Higher order terms can be similarly obtained but they are too long to presented here. Substituting them further into Eqs. (25), the deflection angle in the Hayward BH spacetime in the SFL becomes 2y n,0 η n+1 2 i 2 [ n+1 2 ]+ 1 2 (n + 1) + higher order terms.
(56) In Fig.  3 (a), we plot the relative difference |∆φ H /∆φ H,num − 1| between ∆φ H obtained using Eq. (56) with the truncation of the summation to order 0, 2, 4, 6, 8, 10 and 12 (from bottom to top) and the corresponding numerical value ∆φ H,num obtained using numerical integration of the definition Eq. (5) as functions of b as it approaches b c from above. It is seen that as the order of truncation in the summation of Eq. (56) increases, the deviation of our analytical result from the true value obtained by the numerical integration of Eq. (5) decreases, by a factor of ∼ 3 for every two orders. At order 12, the relative deviation is only about 10 −5 when 1 − b c /b is about 10 −7 , which corresponds to the ∆φ H ≈ 3π, i.e., the trajectory loops 1 circle around the BH. Moreover, for each fixed order, the deviation decreases monotonically as b approaches b c . Both the above features match the fact that our result (56) is essentially a perturbative series approximating the true ∆φ.
In Fig. 3 (b), the coefficients C 0 , D 0 and ∆φ H as functions of the parameter l for different values of v are plotted. The increase of C 0 and decrease of D 0 as the charge parameter l increases are very similar to the effect of electric charge Q in the case of deflection in the SFL in RN BH spacetime [37]. This plot also shows that the ∆φ H increases as l increases, which is in accordance to the decrease of b c in Fig. 2. On the other hand, although the decrease of the velocity from 1 to 0.9 in this case still increases ∆φ H , its effect is much weaker than the increase of l, unlike the case in Fig. 2. This is a reflection that the decrease of v indeed increases the particle sphere size and therefore the looping of the trajectory happens at a larger radius, reducing ∆φ H 's amount of increase.
Substituting the ∆φ H in Eq. (56) into Eq. (40), we are able to solve the apparent angles θ n,s . In Fig. 4 (a) and (b), we plot the Einstein rings in the SFL, which are located at θ n,s when ∆φ sd = π. Here we assume that m = 4.12 × 10 6 M is the SgrA* BH mass, r s = r d = 8.12 [kpc] is the distance to SgrA* BH [46]. However, note that in the zoom-out large plot, different rings correspond to the relativistic Einstein rings of different l (from outer to inner circle in Fig. 4 (a), l/m = 0, 1/2, 3/4) or different v (from inner to outer circle in Fig. 4 (b), v = 1, 0.9, 0.8 respectively). Only in the zoom-in insets the Einstein rings for a fixed l or v are resolved between θ 1,s and θ 2···∞,s . In general, the first relativistic ring θ 1,s is always larger than higher order ones θ 2···∞,s , separated from them by roughly O(0.1) [µas]. Therefore, these higher order rings are not even observationally resolvable in the near future. Their locations are practically the same as the location of the BH shadow of the Hayward BH. In Fig. 4 (c) top and bottom plots respectively, the dependence of the outmost ring location θ 1,s and the inner overlapped ring location θ 2···∞,s on the spacetime parameter l and signal variable v are plotted. It is clear that as l increases or v decreases, the Einstein ring and shadow sizes decrease or increase respectively. Moreover, change of v has a larger effect than that of l on θ n,s , in agreement with their effects on b c as indicated in Fig.  2. These features in the Hayward spacetime are similar to the effect of charge and velocity in the RN spacetime [37]. The magnification of the relativistic images in the Hayward BH spacetime can be obtained using Eq. (44). In Fig. 5 (a), we plot the magnification of the images for different n, ∆φ sd and r s /r d , and in Fig. 5 (b) the magnification as a function of l/m and v. It is seen from Fig.  5 (a) that as expected, for all ∆φ sd and r s /r d , the magnification µ 1,s of the outer-most image is always much larger than that of the inner images. Because the ratio of magnifications of the two successive images is given by exp(−2π/C 0 ) and C 0 < 2, µ 1,s is indeed larger than µ ,s , as predicted by Eq. (47). The three subplots also verify the three possible cases found in Eq. (45). The existence of these three cases is robust and purely due to the geometric variables ∆φ sd , r s and r d . Fig. 5 (b) shows the effect of l and v to the magnification. Given that both the decrease of v and increase of l will increase the value of C 0 , as can be seen from the Fig. 3 (b), using Eq. (44) then it is immediately clear that these changes will increase the magnification too. This implies that the the Hayward BH with larger charge or signal with slower velocity tends to have larger magnification of the relativistic images, which is confirmed in Fig. 5 (b).

V. DISCUSSIONS
In this work, we developed a perturbative method to compute the deflection angle of both null and timelike signals in the SFL in the SSS spacetimes. In doing so, the effects of finite source/detector distances are also taken into account. The result takes a quasi-power series form of Eq. (24), i.e., where a = 1 − bc b is the infinitesimal expansion parameter as b → b + c and the series coefficients of a n contain a weak logarithmic divergence. The a 0 order for light ray has been well known and widely applied in many spacetimes. The form of the first few orders in Schwarzschild spacetime was also known for light signal by expanding the elliptical functions representing the deflection angle [40]. However, the general form of ∆φ for higher order terms for general SSS spacetime, especially for timelike signal is first derived here. Using this result, we connected the behavior of the spacetime near the particle (or photon) sphere, the particle velocity and the finite distances to the deflection angle in the SFL.
Using a new set of GL equations, we were able to solve the apparent angles and their magnifications of the relativistic images in the SFL for signals with arbitrary velocity. A simple formula for the shadow size was found in terms of the metric functions, particle sphere size and detector radius.
Applying the perturbative method and results to the Hayward spacetime, we verified that the series approaches the true deflection angle as the truncation order increases and b → b + c . As the charge parameter l increases, the particle sphere size decreases, the deflection angle increases, the apparent angle decrease and magnification increases. On the other hand, as the decrease of the signal velocity from that of light, the particle sphere size, the deflection angle, the apparent angle and magnification all increase.
where T n are given in Eqs. 35 They can be proven by using successive changes of variables ξ = a + t 2 and t = √ a cot(s) and then use integral formula in Ref. [47]. Their correctness can also be directly checked numerically.