We revisit the gravitational lensing of light or gravitational waves by Schwarzschild black hole in geometric optics. Instead of a single massless particle we investigate the collective behavior of a congruence of lightgravitational rays described by the geodesic deviation equation GDE. By projecting on the NewmanPenrose tetrad GDE is decoupled and we find an analytical Dysonlike series solution in the weak deflection and thin lens limits. Based on such a solution we study the evolution of crosssectional area and axis ratio. Finally we reproduce the magnification and axis ratio of the lensing images up to the second order of weak deflection approximation and improve some missing corrections in previous works.
['General Relativity and Quantum Cosmology']
Schwarzschild Lensing From Geodesic Deviation
{'No Header': '2', 'Schwarzschild Lensing From Geodesic Deviation': 'Zhao Li, 1, 2, 3, ∗ Xiao Guo, 4 Tan Liu, 4, 5 Tao Zhu, 6 and Wen Zhao 1, 2, †  \n- 1 Department of Astronomy, University of Science and Technology of China, Hefei, Anhui 230026, China  \nSchool of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, China  \n3  \nDepartment of Physics, Kyoto University, Kyoto 606-8502, Japan  \n4 School of Fundamental Physics and Mathematical Sciences,  \nHangzhou Institute for Advanced Study, University of Chinese Academy of Sciences,  \nNo.1 Xiangshan Branch, Hangzhou 310024, China  \n5 University of Chinese Academy of Sciences, 100049/100190 Beijing, China  \n6 Institute for Theoretical Physics and Cosmology,  \nZhejiang University of Technology, Hangzhou, 310032, China,  \nUnited Center for Gravitational Wave Physics (UCGWP),  \nZhejiang University of Technology, Hangzhou, 310032, China  \nWe revisit the gravitational lensing of light or gravitational waves by Schwarzschild black hole in geometric optics. Instead of a single massless particle, we investigate the collective behavior of a congruence of light/gravitational rays, described by the geodesic deviation equation (GDE). By projecting on the Newman-Penrose tetrad, GDE is decoupled, and we find an analytical Dyson-like series solution in the weak deflection and thin lens limits. Based on such a solution, we study the evolution of cross-sectional area and axis ratio. Finally, we reproduce the magnification and axis ratio of the lensing images up to the second order of weak deflection approximation and improve some missing corrections in previous works.', 'I. INTRODUCTION': "Gravitational lensing is one of the most important phenomena predicted by general relativity (GR), which occurs when a massive object lies between a distant source and Earth [1-3]. The gravitational field bends the light rays from the source, creating multiple images, or producing arcs and rings [4]. A direct consequence is magnification, allowing us to observe objects that would otherwise be too far away and too faint to be seen at high redshift [5]. Gravitational lensing offers a powerful tool to probe the Universe and test gravity at the astrophysical scale [6-9]. Because the lensing effect directly depends on the lens object's mass distribution, including both visible and dark matter, it allows researchers to measure the dynamical mass of the galaxy and galaxy cluster, and study the distribution of dark matter [10].  \nLike light, gravitational waves (GWs) are lensed when passing through massive objects [7, 11-17]. Several studies have begun to search the lensed GW signal in the current catalog by LIGO/Virgo/KAGRA network, however, finding no strong evidence for lensing imprint [18-24]. As predicted, the future third-generation ground-based GW observatory is expected to detect hundreds of lensing GW events, providing rich information on the gravitational theory, cosmic structure, dark matter, and so on [25-27].  \nAn ideal situation is that the typical wavelength of the lensed signal is much shorter than the background curvature scale. Through eikonal expansion, one finds the behavior of light and GWs are the same as the massless particle in the lowest order. The photon and graviton paths are nothing but the null geodesics [28-30]. In this framework, the lensing process is described by the lens equation [31, 32], for a given surface mass density model of the lens object. The main observational quantities of the lensing events consist of magnification and shear, which can be calculated by the Jacobian matrix [31, 32].  \nIn this work, we revisit the lensing process by a Schwarzschild black hole through a different viewpoint. The lensed signal is not made up of a single photon/graviton, but a congruence of light/gravitational rays [33]. During the propagation, the congruence cross-section area and axis ratio are significantly affected by the tidal of the Schwarzschild lens [34, 35]. We study the evolution of geodesic congruence by solving the geodesic deviation equation (GDE) [2]. Under the weak deflection limit, where the impact distance is much larger than the gravitational radius of the lens [31, 32], we find the analytical solution through a Dyson-like series expansion [36]. Based on that, we analytically  \n∗  \n†  \nreproduce the magnification and axis ratio of the images. In addition, we improve the previous lens equation and the corresponding results, in which some higher-order corrections are missed.  \nThis paper is organized as follows. In Sec II, we review the GDE and project it onto the Newman-Penrose (NP) tetrad. We present the Dyson-like series solution to the GDE and the physical meaning of the optical scalars in Sec III. In Sec IV, we investigate the Schwarzschild lensing, especially, the evolution of the deviation vector and the optical scalars, and most importantly, we reproduce the magnification and shear of the Schwarzschild lens. We take a brief conclusion and discussion in Sec V. We do not consider the cosmological background and redshift in this work. Throughout the paper, we work in geometric units in which c = G = 1, where c is the speed of light in the vacuum and G is the gravitational constant.", 'II. GEODESIC CONGRUENCE': 'In this section, we investigate the behavior of the geodesic congruence [33]. The arm length of GW detectors and the aperture of optical telescopes are much smaller than the scale of lens objects. Therefore, the detected light and gravitational rays are very close to each other and approximately located inside a common geodesic congruence.  \nFor simplicity, we only study two neighboring points, whose trajectories are denoted by  \nx α ( λ, µ a ) , and x α ( λ + δλ, µ a + δµ a ) , (1)  \nwhere λ is the affine parameter and µ a is the a -th motion constant of massless particle in the given background. The difference δλ and δµ a are kept sufficiently small. The geodesic deviation vector is defined as the difference between such two trajectories, and can be expanded in terms of small parameters δµ a and δλ ,  \nξ α ≡ x α ( λ + δλ ; µ a + δµ a ) -x α ( λ ; µ a ) = ∂x α ∂λ δλ + ∑ a ∂x α ∂µ a δµ a + O (2) . (2)  \nwhere O (2) means the second and higher order of δλ , ξ α , and δµ a . We denote the tangent vector of the geodesics as k µ ≡ ∂x µ /∂λ , which satisfies the geodesic equation, k α ( ∇ α k ν ) = 0. Then the difference in wavevectors associated with each trajectory is expanded by the deviation vector,  \nδk ν = ( ∂ α k ν ) ξ α + O (2) . (3)  \nMeanwhile, this difference can also be expanded in terms of small quantities δλ and δµ a ,  \nδk ν = ∑ a ∂k ν ∂µ a δµ a + ∂k ν ∂λ δλ + O (2) = ∑ a ∂ 2 x ν ∂µ a ∂λ δµ a + ∂ 2 x ν ∂λ 2 δλ + O (2) = ∂ ∂λ [ ∑ a ∂x ν ∂µ a δµ a + ∂x ν ∂λ δλ ] + O (2) = ∂ξ ν ∂λ + O (2) = k α ∂ α ξ ν + O (2) , (4)  \nCombining Eqs. (3) and (4) up to the linear order of δλ and δµ α , we get [34, 35]  \nξ α ∇ α k ν = k α ∇ α ξ ν , (5)  \nwith ∇ α being the covariant derivative operator compatible with the background metric. An equivalent explanation to Eq. (5) is the Lie derivatives of deviation vector along geodesic vanish, i.e., L k ξ ν = k α ( ∇ α ξ ν ) -( ∇ α k ν ) ξ α = k α ( ∂ α ξ ν ) -( ∂ α k ν ) ξ α = 0. Multiplying k ν on both sides of Eq. (5), and combining the geodesic equation, one gets conservation law that  \nk α ∇ α ( k ν ξ ν ) = 0 , (6)  \nSetting the inner product k ν ξ ν as zero means that the points x α and x α + ξ α belong to the same wavefront, with equal phases, defined by k µ ≡ -∇ µ Φ [2]. This is proved by  \nΦ( x α + ξ α ) ≈ Φ( x α ) + ( ∂ α Φ) δx α = Φ( x α ) -k α ξ α = Φ( x α ) . (7)  \nThe GDE dominates the evolution of the congruence, that is  \nD 2 ξ ν Dλ 2 = -R ν ραβ k ρ k β ξ α , (8)  \nwith D/Dλ ≡ k α ∇ α and R ν ραβ being the background Riemann tensor. To solve Eq. (8) more conveniently, we construct the NP tetrad along the null geodesics, whose four legs are denoted as [3, 37]  \ne ( a ) = { k , n , m , ¯ m } . (9)  \nThe subscript ( a ) is the tetrad indices. The first leg is the wavevector of the lensed waves, defined as the gradient of a scalar field denoted by Φ. The second one is along the spatially opposite direction. The third and fourth legs are orthogonal to the propagation direction of the waves, they are complex conjugate to each other. The plane spanned by m and ¯ m is generally called the polarization plane, representing the oscillation direction of light or GWs. The orthogonality of the NP tetrad is presented as -k · n = m · ¯ m = 1, and other possible inner products vanish. Most importantly, all of the tetrad legs should be required to be parallel-transported along the geodesics, k α ∇ α e µ ( a ) = 0. We project the deviation vector onto the NP tetrad,  \nξ α = ξ ( a ) e α ( a ) , with ξ ( a ) = ξ α e ( a ) α . (10)  \nSetting ξ α to satisfy the equal-phase condition, we have ξ n = k α ξ α = 0 [53]. Substituting Eq. (10) into Eqs. (5, 8), we obtain [34, 35, 37-41]  \nD ξ m = ρξ m + ¯ σξ ¯ m , (11)  \nand  \nD 2 ξ m = ¯ Ψ 0 ξ ¯ m , (12)  \nfor the m -components. The derivative operator is D ≡ k µ ∂ µ = d/dλ . One obtains similar equations for ¯ m -component by taking complex conjugation of these two equations. In Eqs. (11, 12), the Greek letters ρ and σ represent the spin coefficients, defined as ρ ≡ m µ ( ∇ µ k ν ) ¯ m ν and σ ≡ m µ ( ∇ µ k ν ) m ν . These two scalars are also known as optical scalars, describing the geometry and evolution of the geodesic congruence. The Weyl scalar Ψ 0 is defined as Ψ 0 ≡ -R µναβ k µ m ν k α m β . We have assumed the spacetime to be a vacuum throughout this work, and then the Weyl tensor is identical to the Riemann tensor.', 'III. SOLUTION TO DEVIATION VECTOR': "In this section, we present the solution to the geodesic vector. The Weyl scalar Ψ 0 is a real number for Schwarzschild lensing, which simplifies the GDE (12). We consider Eq. (12) and its complex conjunction. By adding and subtracting this pair of equations, and defining ξ ± ≡ ξ m ± ξ ¯ m , we get  \nD 2 ξ ± = ± Ψ 0 ξ ± . (13)  \nThus, Eq. (13) for ξ ± are decoupled. Transforming these two equations into a set of first-order differential equations, we get the equivalent matrix form  \nDY = A · Y , with Y ≡ ( ξ ± D ξ ± ) , and A ≡ ( 0 1 ± Ψ 0 0 ) . (14)  \nIn a typical lensing system, the impact parameter of the photon/graviton paths is much larger than the gravitational radius of the lens object. Meanwhile, the source and observer are far from the lens, and the background curvature is sufficiently approaching zero at these two points. These two facts are usually called weak deflection and thin lens limits. We assume that the lensed wave is emitted from a point source, and propagates freely near the source. Therefore, the initial conditions of the deviation vector can be set as ξ ± = 0 and D ξ ± = k ± . The initial derivative is proportional to the open angle of the congruence near the wave source. In these two limits, the right-hand side of Eq. (14) is always perturbation, and its 0-th solution can be set as  \nY 0 = ( k ± λ k ± ) . (15)  \nThen the Dyson-like series solution is written as [36, 40, 42]  \nY 1 = ∫ λ 0 dλ ' A ( λ ' ) · Y 0 ( λ ' ) , (16)  \nY 2 = ∫ λ 0 dλ ' A ( λ ' ) · Y 1 ( λ ' ) = ∫ λ 0 dλ ' ∫ λ ' 0 dλ '' A ( λ ' ) · A ( λ '' ) · Y 0 ( λ '' ) , (17)  \n· · ·  \n.  \nY n +1 = ∫ λ 0 dλ ' A ( λ ' ) · Y n ( λ ' ) (18)  \nThe solution to ξ ± is finally iteratively given by,  \nξ ± = k ± · ˜ ξ ± , ˜ ξ ± ≡ ∞ ∑ n =0 ( ± 1) n · I n , (19)  \nwhere the integrations I n are  \nI n = ∫ λ 0 dλ ' ∫ λ ' 0 dλ '' Ψ 0 ( λ '' ) I n -1 ( λ '' ) = ∫ λ 0 ( λ -λ ' )Ψ 0 ( λ ' ) I n -1 ( λ ' ) dλ ' , with I 0 = λ. (20)  \nand ˜ ξ ± are real. Because ξ + is real and ξ -is purely imaginary, such that the integration constants k + and k -are real and purely imaginary, respectively.  \nThe GDE (12) is a second-order equation, there should be two independent solutions. Supposing ξ m 1 and ξ ¯ m 2 to be two independent solutions of deviation vector, with different k ± . Therefore, from Eq. (11), the optical scalars are determined by [34]  \n( ρ σ ) = ( ξ ¯ m 1 ξ m 1 ξ ¯ m 2 ξ m 2 ) -1 ( D ( ξ ¯ m 1 ) D ( ξ ¯ m 2 ) ) , (21)  \nand explicitly [43],  \nρ = D ln √ ∣ ∣ ∣ ˜ ξ + ˜ ξ -∣ ∣ ∣ , and σ = D ln √ ∣ ∣ ∣ ˜ ξ + / ˜ ξ -∣ ∣ ∣ . (22)  \nwhich is independent of the integration constant k ± .  \nIt is noted that the null tetrad leg m relates a spacelike tetrad by m = ( e x + e y ) / √ 2. Therefore, ˜ ξ ± is identified to be the scale of geodesic congruence on the e x and e y directions. Such that, the product ˜ ξ + ˜ ξ -is proportional to the cross-sectional area, and ˜ ξ + / ˜ ξ -represents the axis ratio of the congruence. We can define such area and axis ratio as  \nA ≡ ˜ ξ + ˜ ξ -, and R≡ ˜ ξ -/ ˜ ξ + , (23)  \nThe corresponding evolution equations are derived as [34, 43, 44]  \nD A = 2 ρ A , and D R = -2 σ R . (24)  \nThe physical meanings of the optical scalars are the relative variation rate of the cross-sectional area and the axis ratio of the null congruences.", 'A. Null tetrad and Weyl scalars': 'In coordinate { t, r, θ, φ } , the line element of Schwarzschild black hole with mass M is [1-3]  \nds 2 = -( 1 -2 M r ) dt 2 + ( 1 -2 M r ) -1 dr 2 + r 2 ( dθ 2 +sin 2 θdφ 2 ) , (25)  \nand the tangent vector of the null geodesics equations are  \nk µ = { ( 1 -2 M r ) -1 , U r , 0 , L r 2 } , (26)  \nDue to the spherical symmetry of the Schwarzschild background, the orbits of the massless test particle are always on the equatorial plane of the black hole, with θ = π/ 2. L is the angular momentum of the photon or graviton. The radial potential U r is defined as  \nU r ≡ ± √ 1 -L 2 r 2 ( 1 -2 M r ) , (27)  \nwhere the overall sign is negative and positive when the particle moves toward or backward to the pericenter, respectively. Following the steps presented in Ref. [45] and Appendix A, the null tetrad is constructed as  \nn µ = { r 2 2 L 2 [ -2 λ r U r + ( 1 + λ 2 r 2 )]( 1 -2 M r ) -1 , r 2 L 2 [ U r 2 ( 1 + λ 2 r 2 ) + λ r ] , 0 , -1 2 L ( 1 -λ 2 r 2 ) } , (28)  \nand  \nm µ = 1 √ 2 { r L ( λ r + U r )( 1 -2 M r ) -1 , -r L ( 1 + λ r U r ) , i r , -λ r 2 } . (29)  \nThe Weyl scalar Ψ 0 can be computed as  \nΨ 0 = 3 ML 2 r 5 , (30)  \nwhich is a real number, and ensures that we can apply the framework in Sec III to the Schwarzschild lensing.', 'B. Analytical solution under the weak deflection limit': "We use the Dyson-like series to re-express the solution to the geodesic deviation vector and optical scalars. This work analytically integrates these integrations under weak deflection and thin-lens assumptions. In these limits, the angular momentum L (with mass dimension) is much larger than the Schwarzschild radius, therefore, the impact parameter b ≡ L/M ≫ 1. The distances of the emitter and observer are much larger than the impact distance, i.e. r s , r o ≫ L . Ref. [3] and Appendix B discuss the solution to the null geodesic equation.  \nTo perform the integration (20), it is convenient to express the affine parameter in terms of a new parameter x , defined as x ≡ u/u m , with u ≡ M/r and u m ≡ M/r m , where r m is the pericenter radius of the photon/graviton path. The solution to λ = λ ( x ) is presented in Appendix C, with two different forms in intervals [ λ s , λ m ] and [ λ m , λ o ], denoted by λ I ( x ) and λ II ( x ). λ s , λ m , λ o , are the affine parameter λ evaluated at the source, pericenter, and observer positions. Such that the integration becomes different forms in these two intervals. When λ s < λ < λ m ,  \nI n ( x ) = ∫ x x s [ λ I ( x ) -λ I ( x ' )]Ψ 0 [ λ I ( x ' )] I n -1 [ λ I ( x ' )] dλ I ( x ) dx ' dx ' , ( n = 1 , 2 , 3 , · · · ) . (31)  \nWhen λ m < λ < λ o ,  \nI n ( x ) = ∫ 1 x s [ λ II ( x ) -λ I ( x ' )]Ψ 0 [ λ I ( x ' )] I n -1 [ λ I ( x ' )] dλ I ( x ) dx ' dx ' + ∫ x 1 [ λ II ( x ) -λ II ( x ' )]Ψ 0 [ λ II ( x ' )] I n -1 [ λ II ( x ' )] dλ II ( x ) dx ' dx ' , ( n = 1 , 2 , 3 , · · · ) . (32)  \nThe analytical expressions of integrations I n are shown in Appendix D. Due to the complex form of these expressions, we performed numerical calculations to more clearly demonstrate their behavior.  \nThe typical configuration of gravitational lensing under the weak deflection and thin lens limits is depicted in Fig. 1. The two panels correspond to the two possible photon/graviton paths, forming two images in the geometric optics.  \nFor convenience, we denote the left case as image 1 and the right one as image 2. In Fig. 1, O, L, and S represent the observer, lens object, and wave source. The light blue line approximately characterizes the photon/graviton path, and I is the lensing image on the source plane. Correspondingly, OS is the free-propagation path. And α , β , and θ are the deflection angle, the angular coordinate of the source and image, that is related by lens equation, discussed in Appendix B. D L , D S , and D LS are the distances between observer, lens plane, and source plane. In our numerical calculation, we set D L = 780 M and D LS = 750 M . The well-known Einstein angle is defined by [4]  \nθ E ≡ √ 4 M D LS D S D L , (33)  \nthe value in our setup is about 0 . 0501381 rad ≈ 2 . 8727 · , which is a small angle, consistent with the weak deflection assumption. After giving the distances D L and D LS (in the unit of lens mass), the unique free parameter is the angular coordinate of the wave source, β . For convenience, the normalized coordinate is β 0 ≡ β/θ E . We set β 0 = { 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 } and the corresponding values of the image position, deflection angle, and impact parameter are listed in Table I. In the weak deflection limit, such quantities satisfy β, θ, α, u m ≪ 1 and b ≫ 1.  \n)  \nFIG. 1: Geometrical configuration of gravitational lensing system (left: for image 1, right: for image 2). The points O, L, and S represent the observer, lens object, and wave source. The straight line ON is the optical axis. The planes passing through the lens and source, and being vertical to the optical axis are the lens plane and source plane. The distances from the observer to the lens and source plane are denoted by D L and D S , and then D LS = D S -D L . α is the deflection angle. β and θ are the angular coordinates of the source and image.  \n<!-- image -->  \nTABLE I: We list all of the example values of the lensing parameters involved in our discussion. We have fixed D L = 780 M and D LS = 750 M , and the system is completely determined by source position. The first two columns are the (normalized) source coordinates. The rest of the columns are the normalized image coordinates, deflection angle, impact parameter, magnification, and axis ratio, for the two images produced by the Schwarzschild lensing.  \nIn Fig. 2, we display the evolution of the first four terms in the Dyson-like series in Eq. (20). The first two terms I 0 = λ , shown as a straight line with slope 1, and I 1 are the dominant terms of the geodesic deviation vector. There  \nis an obvious turning point at λ/M ∼ 750, because the external gravitational field deflects the photon/graviton path by slightly different angles. The subsequent two terms, I 2 and I 3 , are almost always zero because they represent highorder weak-deflection correction. This also indicates that the Dyson-like series converges under the weak-deflection limit.  \nFIG. 2: The evolution of each term in Dyson-like series (19) along the geodesic. These curves are plotted based on the analytical solution shown in Eqs. (C2, C4, D1, D2, D3).  \n<!-- image -->  \nFig. 3 shows the evolution of the transversal components of the geodesic deviation vector along the null geodesics. Before entering the external gravitational field, the evolutions of ˜ ξ + and ˜ ξ -are almost equal. However, the slopes of ˜ ξ + and ˜ ξ -are obviously different after leaving the gravitational field. From the expressions of tetrad leg m µ , we conclude that ˜ ξ + lies on the equatorial plane of the Schwarzschild black hole, and ˜ ξ -is vertical to that. Due to the tidal force, ˜ ξ + always increases, and the increasing rate gets larger after leaving the lens object than before entering. Distinguishingly, the growth of ˜ ξ -is always suppressed. For relatively large b , the tidal force is strong enough to make ˜ ξ -decrease rapidly enough, as shown in the second line of Fig. 3, ˜ ξ -arrives at zero at the moment, where the cross-sectional area and axis ratio tend to 0. This point is usually called a caustic point, denoted by λ c and determined by ˜ ξ -( λ c ) = 0. In this case, the detected lensing image is inverted (negative parity). Conversely, when the impact parameter is relatively small, there is no caustic point, and ˜ ξ -remains positive on the entire path, ultimately resulting in an upright (positive parity) image. Besides the approximate analytical solution, we implement the numerical solution to Eq. (13), with initial values being ξ ± ( λ = 0) = 0 and D ξ ± ( λ = 0) = 1. These results are also plotted in Fig. 3, verifying that the Dyson-like series solution is valid under the weak deflection and thin lens limits.  \nAs shown in Eq. (24) and around it, we can define the cross-sectional area ( A ) and axis ratio ( R ) from the deviation vector. Fig. 4 depicts their evolution along the path. The area and axis ratio are always positive for image 1 (the left panels), but meet zero at the caustic point for image 2 (the right panels). The relative changing rates of A and R are described by optical scalars, ρ and σ , respectively, which are shown in Fig. 5. When closing to the wave source (point like as we have assumed), all of the geodesics in congruence focus on a point, and subsequently, ρ is infinity. At the starting point, ˜ ξ + and ˜ ξ -are set to the same initial conditions, such that σ is keeping to zero before approaching the lens object. More importantly, in image 1, these two optical scalars are continuous, but there is a discontinuity in image 2, caused by the caustic point.  \nThe weak deflection and thin-lens approximation bring errors in the analytical solutions. For example, one can find slight deviations between the slopes of the numerical and analytical solutions to ˜ ξ ± (see Fig. 3). Such error is amplified when calculating σ , especially when approaching the observer (see the left-bottom panel in Fig. 5). It is noted that the example parameters in Table I significantly deviate from the realistic situation. A typical lens object has mass M ∼ 10 8 M ⊙ and redshift z ∼ 1. The angular separation of the image and lens is about θ ∼ 1 '' . We can estimate the distance and impact parameter for such a system as D L ≈ 4 . 3 × 10 3 Mpc ≈ 9 . 0 × 10 14 M and b ∼ 4 . 8 × 10 -6 D L . This  \nestimation indicates that the impact distance is much greater than the Schwarzschild radius, and D L is much greater than the impact distance, ensuring the validity of the weak deflection and thin lens approximations.  \nFIG. 3: The evolution of the transversal components of the geodesic deviation vector, ˜ ξ ± from different source positions and images. The red curves show the evaluated results based on the analytical expressions from the Dyson-like series (20). The blue curves are obtained by numerically solving the GDE (13), together with the radial geodesic equation in Eq. (26).  \n<!-- image -->", 'C. The magnification and axis ratio of the lensing images': 'In the last subsection, we investigate the behaviors of the geodesic deviation vector, cross-sectional area/axis ratio, and optical scalars along the lensing path. In this subsection, we study the observational imprints by such evolution. The magnification of the lensing image is defined as the ratio between the lensed and unlensed energy flux at the observer. In geometric optics, the energy flux of light and GWs is covariantly conserved within the congruence, and then proportional to 1 / A . Therefore, the magnification is rewritten as  \nµ = ̂ A ( λ o ) A ( λ o ) . (34)  \n̂ A means the cross-sectional area of a congruence freely propagating along the unlensed path (OS, see Fig. 1). By setting the initial conditions as ̂ A ( λ = 0) = 0 and D ̂ A ( λ = 0) = 1, its solution is given by  \n̂ A ( λ ) = 1 OS 2 , with OS = D S √ 1 + tan 2 β = r o √ 1 + tan 2 β [ 1 + √ 1 -(1 + tan 2 β ) ( 1 -r 2 s r 2 o ) ] . (35)  \nWe recall that r s and r o are the radial Schwarzschild coordinates of the source and observer, relating D L and D S by  \nr o = D L , and , r s = D 2 S (1 + tan 2 β ) -2 D L D S + D 2 L (36)  \nThe area of the lensed congruence is calculated from the Dyson-like series (20). At the observer position, we get the expressions of the first four terms in this series,  \nI 0 ( λ o ) = λ o ≈ M u m ( 1 x s + 1 x o )[ 1 -1 2 x s x o ( 1 -4 u m x s + x o )] , (37)  \nand  \nI 3 ( λ o ) ≈ M x s x o · 16 5 u 2 m . (40)  \nUnder the weak deflection and thin-lens assumptions, u m ≪ 1, and x s,o ≡ r m /r s,o ≪ 1 with the same order of u m , we have expanded the expressions in terms of small parameters u m , x s , and x o , up to the zeroth order, in Eqs. (37 40). Inserting Eqs. (37 - 40) into Eq. (19), and then into the definition Eq. (23), the area A ( λ o ) is obtained. Then, combining this result with the unlensed area in Eq. (35), the magnification is given from its definition Eq. (34).  \nHowever, the parameters used here are different with those in the standard textbook [4], and we transform them into the commonly used lensing parameter as follows. (1) Transform u m into b by Eq. (B3), and then, from the geometry shown in Fig. 1, b = ± D L sin θ/M for image 1 and 2, respectively. (2) Transform x s and x o into β , D L , and D LS by their definition and Eq. (36). (3) Substitute the relation between β and θ from lens equation (B12, B12). (4) In the weak deflection and thin lens limit, we expand the image position as  \nθ = θ E [ θ 0 + θ 1 ϵ + θ 2 ϵ 2 + O ( ϵ 3 ) ] , (41)  \nwhere the bookkeeper is of order,  \nFIG. 4: The evolution of the area and axis ratio of the geodesic congruence. Same with Fig. 3, the blue and red curves are obtained from the analytical solutions and corresponding numerical solutions to Eq. (13).  \n<!-- image -->  \nI 1 ( λ o ) ≈ 2 M x s x o { 2 + ( 45 π 16 -4 ) u m + [( 192 5 -135 π 16 ) u 2 m +2 u m ( x s + x o ) -( x 2 s -x s x o -x 2 o ) ]} , (38)  \nI 2 ( λ o ) ≈ M x s x o { 15 π 8 u m + u m [( 48 -45 π 8 ) u m -4( x s + x o ) ]} , (39)  \nϵ ∼ θ E 4 D , with D ≡ D LS D S . (42)  \nFIG. 5: The evolution of the spin coefficients ρ and σ . Same with Fig. 3, the blue and red curves are obtained from the analytical solutions and corresponding numerical solutions to Eq. (13). The vertical dotted lines in the right panels represent the caustic points. The corresponding values are listed in Table I.  \n<!-- image -->  \nAfter the above operations, we obtain the final result of the magnification, which is expanded in powers of ϵ as  \nµ = θ 4 0 θ 4 0 -1 ∓ 15 πθ 3 0 16(1 + θ 2 0 ) 3 · ϵ + θ 2 0 384(1 -θ 2 0 ) 2 (1 + θ 2 0 ) 5 { -1024(3 -2 D 2 ) +1024(9 + 12 D14 D 2 ) θ 2 0 +[4096(3 -D 2 ) -2025 π 2 ] θ 4 0 -[ 2048(6 + 15 D17 D 2 ) -2025 π 2 ] θ 6 0 -1024(9 + 6 D8 D 2 ) θ 8 0 +1024 D (3 + 18 D20 D 2 ) θ 10 0 +6144 D (1 -D ) θ 12 0 } ϵ 2 + O ( ϵ 3 ) , (43)  \nThe leading-order result can be found in several standard textbooks, e.g., [4]. For the case where θ = θ E , and θ 0 = 1, the magnification tends to infinity. The positions with infinite magnification is defined as the critical curve. When the image is located at the outside/inside of the critical curve, the magnification is positive/negative, as shown in Table I, representing whether the image is upright (positive parity) or inverted (negative parity). In the second-order terms, the ± corresponds to image 1 and image 2. The higher-order corrections are derived by Ref. [31], through the Jacobian matrix. In Appendix E, we compare our result (43) with the previous work [31].  \nAnother direct observational quantity is the axis ratio of the images, describing the shape distortion of the background wave source. Similar to the deduction of magnification, the axis ratio R ( λ o ) is given by,  \nR = θ 2 0 -1 θ 2 0 +1 ∓ 15 π ( θ 2 0 -1)(1 + 3 θ 2 0 ) 16 θ 0 (1 + θ 2 0 ) 3 · ϵ + ( θ 2 0 -1) 768 θ 2 0 (1 + θ 2 0 ) 5 { -9(2048 -225 π 2 ) -[4096(18 + 3 D4 D 2 ) -8100 π 2 ] θ 2 0 -3[4096(6 + D-D 2 ) -3375 π 2 ] θ 4 0 -8193 D 2 θ 6 0 +6144(3 -2 D +2 D 2 ) θ 8 0 -4096 D (3 -4 D ) θ 10 0 } ϵ 2 + O ( ϵ 3 ) . (44)  \nWhen the lensing image is located at the critical curve, with θ 0 = 1, the axis ratio is zero. This means the extended  \nsource is distorted to be a line segment, with zero area and infinite amplitude. Similar to the magnification, the axis ratio is positive/negative for the images located outside/inside of the critical curve, representing whether the image is upright (positive parity) or inverted (negative parity). The numerical results of R are listed in Table I.  \nThe magnification and axis ratio are related to the convergence ( κ ) and shear ( γ 1 and γ 2 ) by  \nµ = 1 (1 -κ ) 2 -γ 2 , and R = 1 -κ -γ 1 -κ + γ , (45)  \nwith γ = √ γ 2 1 + γ 2 2 . The convergence depends on the lens mass density within the geodesic congruence. The shear γ 1 characterizes the stretching of the lensing images and γ 2 characterizes its orientation. For the Schwarzschild lens, that is vacuum and axisymmetric, κ = 0, γ 2 = 0, and γ = γ 1 . Through the Jacobian matrix [4], the shear of the point-mass lens is calculated as γ = 1 /θ 2 0 at the leading order of weak deflection approximation, resulting in µ = θ 4 0 / ( θ 4 0 -1) and R = ( θ 2 0 -1) / ( θ 2 0 +1), as shown in the leading terms of Eqs. (43) and (44).', 'V. CONCLUSION': 'In this work, we revisit the Schwarzschild lensing of light or GWs from geodesic deviation. The detected lensed signal consists of a congruence of light/gravitational rays, and their collective behavior is described by the GDE. To solve it, we project the GDE onto the NP tetrad constructed along the null geodesics and get a pair of decoupled equations for ξ + and ξ -(13) for real Weyls scalar Ψ 0 . Similar to the Dyson series, we present the solution to the transversal components of the deviation vector in Eqs. (19) and (20). In addition, the relationship between the deviation vector and the optical scalars, ρ and σ , is presented in Eq. (11). The physical meaning of these two scalars is the relative changing rate of the congruence cross section and the axis ratio [see Eq. (22)].  \nSubsequently, we consider the lensing process in the weak deflection limit, where the impact parameter is sufficiently large. The general signature is shown in Fig. 1, two images are formed by Schwarzschild lensing. In this assumption, the Dyson-like series is calculated from Eq. (31, 32) analytically. We depict such results in Fig. 2, and compare them with the numerical solution for a group of parameters listed in Table I. And we find that I 0 and I 1 are dominant and I 2 , I 3 are close to zero, implying the series (19) is convergent in weak deflection limit. Fig. 3 shows the evolution of the deviation vector. ˜ ξ + increases more rapidly after lensing and keeps positive. And ˜ ξ -increases more slowly or even decreases after lensing, meeting the caustic point, and finally forms an inverted image at the observer. Correspondingly, Figs. 4 and 5 show the cross-sectional area, axis ratio, and optical scalars, defined in Eqs. (23, 24).  \nIn subsection IV C, we reproduce the magnification from the geodesic deviation [see Eq. (43)], consistent with that from the Jacobian matrix in Ref. [31]. Meanwhile, we find some missing corrections in the previous calculation, e.g., the deflection point deviates from the lens plane and the lens equations are slightly different for the two images. Additionally, we present the result of the axis ratio, equivalently, the shear, of lensing images up to the second weak deflection approximation [see Eq. (44)].  \nThis work is potentially worthwhile for future studies. As discussed in Ref. [35, 44], the calculation of spin coefficients is indispensable when investigating the lensing effects beyond geometric optics, in which richer information on light and GW polarization is presented [46-51]. In this work, we calculate ρ and σ following the procedure in Ref. [34] in the simplest case, the Schwarzschild lensing. Furthermore, the astrophysical lens objects are more complex than a Schwarzschild black hole [52]. Extending the computation to a more general lens model is helpful to the lensing observation and gravitational testing in the future.', 'Acknowledgments': 'We would like to thank Takahiro Tanaka, Emanuel Gallo, Shaoqi Hou, Sam. R. Dolan, Haofu Zhu, Xinyue Jiang, and Donglin Gao for their helpful discussions and comments. This work is supported by Strategic Priority Research Program of the Chinese Academy of Science (Grant No. XDB0550300), the National Key R&D Program of China Grant No.2022YFC2200100 and 2021YFC2203102, NSFC No.12273035 and 12325301 the Fundamental Research Funds for the Central Universities under Grant No. WK3440000004, and the science research grants from the China Manned Space Project with No.CMS-CSST-2021-B01, and China Scholarship Council, No. 202306340128. X. G. acknowledges the fellowship of China National Postdoctoral Program for Innovative Talents (Grant No. BX20230104). T.L. is supported by NSFC No. 12003008. T.Z. is supported in part by the National Key Research and Development Program of China under Grant No. 2020YFC2201503, the Zhejiang Provincial Natural Science Foundation of China under Grant Nos. LR21A050001 and LY20A050002, the National Natural Science Foundation of China under Grant No. 12275238.', 'Appendix A: Null tetrad in Kerr spacetime': 'In this section, we review the main procedure to construct the null tetrad in Kerr spacetime, shown in Ref. [45]. Let us denote the Boyer-Lindquist coordinate as { t, r, θ, φ } , the background metric as g µν , the black hole mass as M , and spin as a in this section. Before construction, it is necessary to summarize several important tensors. The first one is the principal tensor,  \nγ µν =       0 r a 2 sin θ cos θ 0 -r 0 0 ar sin 2 θ -a 2 sin θ cos θ 0 0 a ( r 2 + a 2 ) sin θ cos θ 0 -ar sin 2 θ -a ( r 2 + a 2 ) sin θ cos θ 0       , (A1)  \nsatisfying ∇ α γ µν = g αµ ξ ν -g αν ξ µ , where ξ µ = (1 , 0 , 0 , 0) is the Killing vector, indicating the energy conservation of the particle moving along the geodesics. Its Hoge dual defines the Killing-Yano tensor by  \nf µν ≡ -ˆ γ µν = -1 2 √ -gϵ αβµν γ αβ =       0 -a cos θ ar sin θ 0 a cos θ 0 0 -a 2 sin 2 θ cos θ -ar sin θ 0 0 r ( r 2 + a 2 ) sin θ 0 a 2 sin 2 θ cos θ -r ( r 2 + a 2 ) sin θ 0       , (A2)  \nsatisfying ∇ ( α f µν ) = 0. ϵ αβµν is Levi-Civita symbol.  \nFirstly, we need to construct a set of basis vectors, denoted by { λ (1) , λ (2) , λ (3) , λ (4) } , as follows,  \nλ α (1) = k α , (A3)  \nλ α (3) = γ αβ λ (1) β -A ( r, θ ) λ α (1) , (A4)  \nλ α (2) = γ αβ λ (3) β -B ( r, θ ) λ α (1) , (A5)  \nλ α (4) = f α β λ β (1) , (A6)  \nThe wavevector in Kerr spacetime is presented by  \nλ α (1) = k α = { ( r 2 + a 2 ) 2 -a 2 ∆sin 2 θ -2 Mξar ∆Σ , R Σ , Θ Σ , 2 Mar +(Σ -2 Mr ) ξ csc 2 θ ∆Σ } , (A7)  \nwith ∆ ≡ r 2 (1 -2 M/r + a 2 /r 2 ) and Σ ≡ r 2 + a 2 cos 2 θ . A ( r, θ ) and B ( r, θ ) are two underdetermined functions. ξ, η being the conserved angular momentum, and Carter constant. The first leg is parallel-transported. Letting the third leg be parallel-transported, we get the constraint on the unknown functions,  \ndA dλ = ξ β k β , and dB dλ = h γβ ξ γ k β -1 2 dA 2 dλ . (A8)  \nFrom the above equation, one can find that the function B is determined by directly integrating on the right-hand side once A is given. Based on the above discussion, we can list all of the expressions of the tetrad. The third one is  \nλ α (3) =       { a 2 sin 2 θ/ Σ -[( r 2 + a 2 ) 2 -2 Maξr ] / ∆Σ } A -r ( r 2 + a 2 ) R / ∆Σ -(Θ / Σ) a 2 sin θ cos θ -r ( r 2 + a 2 -aξ ) / Σ -( R / Σ) A a cot θ ( ξ -a sin 2 θ ) / Σ -(Θ / Σ) A -[2 Mar +(Σ -2 Mr ) ξ csc 2 θ ] A/ (∆Σ) + a [(Θ / Σ)cot θ +( r/ ∆)( R / Σ)] ,       (A9)  \nUsing the known Killing vector, the constraint on A ( r, θ ) becomes dA/dλ = -1. We shall suppose that the function A is variable-separable, which gives the form A ( r, θ ) = F ( r ) + G ( θ ). An obvious set of solutions is  \nF ( r ) = -∫ r 2 R dr, G ( θ ) = -∫ a 2 cos 2 θ Θ dθ. (A10)  \nUsing the principal tensor (A1) and Killing vector, the constraint on B becomes  \nB ( r, θ ) = 1 2 [ r 2 -a 2 cos 2 θ -A 2 ( r, θ )] . (A11)  \nAnd then the second leg is  \nλ t (2) = -(1 / 2)( r 2 -a 2 cos 2 θ -A 2 ) { ( r 2 + a 2 ) 2 -a 2 ∆sin 2 θ -2 aMξr } / (∆Σ) + r ( r 2 + a 2 )( R A + a 2 r -aξr + r 3 ) / (Σ∆) + a 2 cos θ [Θsin θA + a cos θ ( a sin 2 θ -ξ )] / Σ , λ r (2) = { 2 rA [( r 2 + a 2 ) -aξ ] + R (Σ + A 2 ) } / (2Σ) , λ θ (2) = [ -2 a cot θA ( ξ -a sin 2 θ ) + Θ A 2 -ΣΘ ] / (2Σ) , λ φ (2) = ( A 2 + a 2 cos 2 θ -r 2 )[ a 2 ξ cot 2 θ +2 aMr + ξr csc 2 θ ( r -2 M )] / (2∆Σ) +2 ar ( R A + a 2 r -aξr + r 3 ) / (∆Σ) -a cot θ ( a cot θ ( ξ -a sin 2 θ ) -A Θ) / Σ , (A12)  \nand the fourth leg is  \nλ α (4) =       a [ R ( r 2 + a 2 ) cos θ -r Θsin θ ∆] / (Σ∆) a cos θ [( r 2 + a 2 ) -aξ ] / Σ r csc θ ( ξ -a sin 2 θ ) / Σ csc θ ( a 2 R sin θ cos θ -r ∆Θ) / (∆Σ)       . (A13)  \nAfter constructing { λ (1) , λ (2) , λ (3) , λ (4) } , the NP tetrad is defined by  \ne α ( a ) = { λ α (1) , 1 C 2 λ α (2) , 1 C 1 √ 2 [ λ α (3) + iλ α (4) ] , 1 C 1 √ 2 [ λ α (3) -iλ α (4) ] } , (A14)  \nwhich satisfies the orthogonality and parallel-transported condition. The normalization constant is  \nC ≡ √ ( a -ξ ) 2 + η. (A15)', 'Appendix B: Schwarzschild lensing': 'One defines u ≡ M/r , and then the photon/graviton trajectory is dominated by [3]  \n( du dφ ) 2 = 1 b 2 -u 2 (1 -2 u ) ≡ f ( u ) , (B1)  \nThe impact parameter is b = | L | /M . When b > 3 √ 3, the function f ( u ) possesses two different zero points, denoted by 0 < u 2 < u 3 , the allowed orbits should be located at interval [0 , u 2 ] and [ u 3 , + ∞ ], where f ( u ) > 0. These first kinds of orbits are open, in which the massless particle moves between an infinite region and a pericenter, corresponding to the gravitational lensing situation. The radial range of allowed motion is r 0 < r < max { r o , r s } , or equivalently, min { u o , u s } < u < u m ≡ M/r m . r o and r s denote the radial coordinates of the wave source and observer. correspondingly, u o ≡ M/r o and u s ≡ M/r s .  \nIn this work, we only focus on the weak deflection limit and thin lens approximation, in which the impact distance ( Mb ) is much larger than the gravitational radius of the lens black hole, and the distances from emitter/observer to lens object are very larger than the impact distance. Using a bookkeeper, ϵ , the order of magnitude impact parameter, source/observer distance are  \nu m ∼ 1 b ∼ ϵ, r m r s , r m r o ∼ ϵ. (B2)  \nThe following calculations are based on the Taylor expansion by the bookkeeper ϵ . Additionally, the explicit relationship between u m and b is given by f ( u m ) = 0,  \nb = 1 u m 1 √ 1 -2 u m = 1 u m [ 1 + u m + 3 2 u 2 m + 5 2 u 3 m + O ( u 4 m ) ] . (B3)  \nFollowing the standard steps, the trajectory equation u = u ( φ ) is  \nu = u m cos φ + u 2 m ( 1 -cos φ +sin 2 φ ) + u 3 m [( 2 -3 4 sin 2 φ ) cos φ -2(1 + sin 2 φ ) + 15 4 φ sin φ ] + O ( u 3 m ) , (B4)  \nUp to order O ( u 3 m ). The deflection angle, α , is obtained from the asymptotic behavior of solution (B4), given by  \nα = 4 u m + ( -4 + 15 4 π ) u 2 m + ( 122 3 -15 π 2 ) u 3 m . (B5)  \nAfter giving the deflection angle (B5), we re-drive the lens equation in this appendix. We first investigate image 1 shown in Fig. 1a, in which the source and image lie on the same side of the optical axis. In △ ASI, applying the sine theorem gives,  \nSI sin α = AI sin ∠ ASI ⇒ NI -SN sin α = OI -OA sin[ π/ 2 -( α -θ )] = OI -OA cos( α -θ ) , (B6)  \nwhere we use ∠ ASI = π -α -∠ AIS = π -α -( π/ 2 -θ ) = π/ 2 -( α -θ ) and then  \n( NI OI -SN OI ) cos( α -θ ) = ( 1 -OA OI ) sin α, (B7)  \nor, equivalently  \n( D S tan θ -D S tan β D S sec θ ) cos( α -θ ) = ( 1 -OA D S sec θ ) sin α. (B8)  \nThen the length of line OA can be expressed from the △ OAL by sine theorem,  \nOL sin ∠ OAL = OA sin ∠ OLA . (B9)  \nHere, the line AL divides ∠ OAS equally, such that the above equation becomes  \nD L sin( π -α ) / 2 = OA sin[ π -θ -( π -α ) / 2] . (B10)  \nUsing the expression of OA in terms of D L , θ , and α ,  \nOA = D L cos( α/ 2 -θ ) cos( α/ 2) . (B11)  \nWe finally arrived at  \n(tan θ -tan β ) cos θ cos( α -θ ) = sin α -2 D L D S cos ( α 2 -θ ) cos θ sin α 2 . (B12)  \nBy replacing θ by -θ , and β by -β (equivalently, replace NI -SN by NI+SN) in Eq.(B12), we get the lens equation for image 2,  \n-(tan θ -tan β ) cos θ cos( α + θ ) = sin α -2 D L D S cos ( α 2 + θ ) cos θ sin α 2 , (B13)  \nin which the source and image lie on the different sides of the optical axis.  \nTo solve the lens equation (B12, B13), we rescale the source position and image position as  \nβ = β 0 θ E , θ = θ E [ θ 0 + θ 1 ϵ + θ 2 ϵ 2 + O ( ϵ 3 ) ] , (B14)  \nwhere the bookkeeper ϵ is defined in Eq. (42). Then the lens equation at leading order is simplified as  \nβ 0 = θ 0 -1 θ 0 , (B15)  \nat the leading order, whose solution is  \nθ 0 = 1 2 [ β 0 ± √ β 2 0 +4 ] . (B16)  \nIn the second and third orders, the solutions to θ 1 and θ 2 are  \nθ 1 = ± 15 π 16(1 + θ 2 0 ) , (B17)  \nand  \nθ 2 = 1 3 θ 0 (1 + θ 2 0 ) 3 { 48 -675 256 π 2 -16 D 2 + [ 8(9 -3 D +2 D 2 ) -675 128 π 2 ] θ 2 0 +40 D 2 θ 4 0 -8(3 -9 D +4 D 2 ) θ 6 0 +8 D (6 -5 D ) θ 8 0 } . (B18)  \nThe plus-minus symbol corresponds to the image 1 and 2, respectively. The leading-order result (B16) is consistent with Ref. [31], but there is an extra plus-minus symbol in first-order correction (B17). And the second-order result (B18) is different, because the deflection point deviates from the lens plane.', 'Appendix C: Solution to the radial motion': 'In weak deflection limit, the radial potential reduces to  \nU r = ± √ 1 -x 2 { 1 -u m ( x 2 1 + x ) -u 2 m [ x 2 (2 + x ) 2 2(1 + x ) 2 ]} + O ( u 3 m ) . (C1)  \nwhere x ≡ u/u m , u = M/r ⩽ u m , u m = M/r m ≪ 1, with r m is the radius of pericenter. When the particle moves toward the pericenter, the above equation is integrated by  \nM -1 λ I = M -1 ∫ r r s U -1 r dr ≈ -1 u m ( ν x -ν s x s ) -( ν 1 + x -ν s 1 + x s ) + u m 2 [ xν (1 + x ) 2 -x s ν s (1 + x s ) 2 +3arcsin x -3 arcsin x s ] , (C2)  \nwith initial condition λ s = 0. For simplicity, we define ν = √ 1 -x 2 , and ν s = √ 1 -x 2 s . When arriving the pericenter, r = r m , where x = 1, the affine parameter is  \nλ m = 1 u m M x s [ 1 + x s ( u m -x s 2 )] , (C3)  \nAnd then, when moving backward, the integration should be  \nM -1 λ II = M -1 { λ m + ∫ r r m U -1 r dr } ≈ 1 u m ( ν x + ν s x s ) + ( ν 1 + x + ν s 1 + x s ) -u m 2 [ xν (1 + x ) 2 + x s ν s (1 + x s ) 2 +3arcsin x +3arcsin x s -3 π ] . (C4)  \nAnd the affine parameter at the observer position has been shown in Eq. (37).', 'Appendix D: The analytical expressions of Eqs. (31) and (32)': 'In this appendix, we present the analytical solution of Eqs. (31) and (32). The first three orders of the Dyson-like series are  \nM -1 I 1 = ± ( 1 + x s x -2 xx s ) ν ± u m 8 xx s (1 + x ) {[ 16(2 + x ) -(61 + 61 x -15 x 2 -7 x 3 ) x s +8 xx 2 s +7(1 + x ) x 2 s ] ν +15(3 νν s ± xx s ) arccos x } ± u 2 m 40 xx s (1 + x ) 2 [ -8(192 + 364 x +167 x 2 ) +(835 + 1709 x +1143 x 2 +326 x 3 +119 x 4 +42 x 5 ) x s +(543 + 1086 x +523 x 2 ) x 2 s +35 x (1 + x ) x 3 s +42(1 + x ) 2 x 4 s ] ν ∓ 3 u 2 m 8 xx s (1 + x )(1 + x s ) { ± 15 [3 + 2 x +(2 + x ) x s ] νν s +(1 + x ) [ -8 x -(8 + 13 x -4 x 2 ) x s -(8 + x -4 x 2 ) x 2 s +4 xx 3 s ] } arccos x ∓ ( x ↔ x s ) , (D1)  \nM -1 I 2 = ± u m 8 xx s { [ 16 x +(1 -3 x 2 ) x s -8 xx 2 s -3 x 2 s ] ν ± 3(5 νν s ∓ xx s ) arccos x } ± u 2 m 8 xx s (1 + x )(1 + x s ) [ -2(96 + 112 x + x 2 -7 x 3 ) -(179 + 151 x -60 x 2 -6 x 3 +5 x 4 ) x s +(139 + 199 x +47 x 2 -15 x 3 -5 x 4 ) x 2 s +(126 + 116 x -22 x 2 -7 x 3 ) x 3 s -(5 + 15 x +7 x 2 ) x 4 s -5(1 + x ) x 5 s ] ν + 3 u 2 m 8 xx s (1 + x )(1 + x s ) { -5 [ 3 + 5 x +3 x 2 +(5 + 7 x +3 x 2 ) x s +3(1 + x ) x 2 s ] νν s ± (1 + x ) [ 10 x +(10 + 11 x -5 x 2 ) x s +(10 -4 x -5 x 2 ) x 2 s -5 xx 3 s ] } arccos x ∓ ( x ↔ x s ) , (D2)  \nand  \nM -1 I 3 = ± u 2 m 40 xx s [ -2(32 + 15 x ) -x (19 -3 x ) x s +(107 + 15 x 2 ) x 2 s +15 xx 3 s +3 x 4 s ] ν + 3 u 2 m ( x + x s ) 8 xx s ( -5 νν s ∓ 2 ± xx s ) arccos x ∓ ( x ↔ x s ) . (D3)  \nThe upper sign corresponds to the case where v s < v < v m , and the lower one to v m < v < v o . ( x ↔ x s ) represents exchanging all x and x s in the formula.', 'Appendix E: Comparison between Eq. (43) and previous result.': "In this appendix, we compare our result on magnification with that provided by Ref. [31] [Eq. (76-79), with A 1 = 4, A 2 = 15 π/ 4, and A 3 = 128 / 3] and discuss the difference. In [31], the magnification is defined by the Jacobian matrix,  \nµ = ∣ ∣ ∣ ∣ sin βdβ sin θdθ ∣ ∣ ∣ ∣ -1 . (E1)  \nIn the leading order, these two results are consistent. But there are some missing corrections when [31] driving the higher-order corrections. Firstly, the images formed by the lensing are at the same and opposite sides to the source, with positive and negative angular momentum. Therefore the lens equation for these two cases is slightly different. This is shown in Eqs. (B12, B13) and Fig. 1. This leads to an extra ± symbol in the first-order correction, which is missed in Ref. [31]. Secondly, the deflection point (denoted by A) deviates from the lens plane (see Fig. 1), resulting in inconsistency in the second order. Thirdly, the lensed and unlensed signals are emitted in slightly different directions. As the point-source assumption in this work, another correction factor cos β/ cos( α -θ ) is needed to consider. After taking these three extra corrections, one obtains the same magnification through geodesic deviation and Jacobian matrix.  \n- [2] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (W. H. Freeman and Company, New York, 1971).\n- [3] S. Chandrasekhar, The mathematical theory of black holes (1983).\n- [4] P. Schneider, J. Ehlers, and E. E. Falco, Gravitational Lenses (1992).\n- [5] K. C. Wong et al. , The Astrophysical Journal Letters 789 , L31 (2014).\n- [6] X.-H. Liu, Z.-H. Li, J.-Z. Qi, and X. Zhang, The Astrophysical Journal 927 , 28 (2022).  \n[7] J. M. Ezquiaga  \net al.  \n, Phys. Rev. D  \n103  \n, 064047 (2021).  \n- [8] Q. Wang et al. , The Astrophysical Journal 969 , 119 (2024).\n- [9] A. Refregier, Annual Review of Astronomy and Astrophysics 41 , 645 (2003).\n- [10] L. V. E. Koopmans et al. , (2009), arXiv:0902.3186 [astro-ph.CO] .\n- [11] M. Grespan and M. Biesiada, Universe 9 , 10.3390/universe9050200 (2023).\n- [12] X. Guo and Y. Lu, Phys. Rev. D 102 , 124076 (2020).\n- [13] R. Takahashi and T. Nakamura, The Astrophysical Journal 595 , 1039 (2003).\n- [14] A. Mishra et al. , Monthly Notices of the Royal Astronomical Society 508 , 4869 (2021).\n- [15] A. K. Meena and J. S. Bagla, Monthly Notices of the Royal Astronomical Society 492 , 1127 (2019).\n- [16] L. Dai and T. Venumadhav, (2017), arXiv:1702.04724 [gr-qc] .\n- [17] Pagano, G., Hannuksela, O. A., and Li, T. G. F., Astronomy and Astrophysics 643 , A167 (2020).\n- [18] K. Haris et al. , (2018), arXiv:1807.07062 [gr-qc] .\n- [19] O. A. Hannuksela et al. , The Astrophysical Journal Letters 874 , L2 (2019).\n- [20] C. McIsaac et al. , Phys. Rev. D 102 , 084031 (2020).\n- [21] X. Liu, I. M. Hernandez, and J. Creighton, The Astrophysical Journal 908 , 97 (2021).\n- [22] R. Abbott et al. , The Astrophysical Journal 923 , 14 (2021).\n- [23] R. Abbott et al. , (2023), arXiv:2304.08393 [gr-qc] .\n- [24] R. K. L. Lo and I. Maga˜na Hernandez, Phys. Rev. D 107 , 123015 (2023).\n- [25] L. Yang et al. , Monthly Notices of the Royal Astronomical Society 509 , 3772 (2021).\n- [26] S. Hou, X.-L. Fan, and Z.-H. Zhu, Monthly Notices of the Royal Astronomical Society 507 , 761 (2021).\n- [27] X. Guo and Y. Lu, Phys. Rev. D 106 , 023018 (2022), arXiv:2207.00325 [astro-ph.CO] .\n- [28] R. A. Isaacson, Phys. Rev. 166 , 1263 (1968).\n- [29] R. A. Isaacson, Phys. Rev. 166 , 1272 (1968).\n- [30] S. Hou, X.-L. Fan, and Z.-H. Zhu, Phys. Rev. D 100 , 064028 (2019).\n- [31] C. R. Keeton and A. O. Petters, Phys. Rev. D 72 , 104006 (2005).\n- [32] M. Sereno and F. De Luca, Phys. Rev. D\n- 74 , 123009 (2006).\n- [33] E. Poisson, A Relativist's Toolkit: The Mathematics of Black-Hole Mechanics (Cambridge University Press, 2009).\n- [34] S. R. Dolan, International Journal of Modern Physics D 27 , 1843010 (2018).\n- [35] S. R. Dolan, (2018), arXiv:1801.02273 [gr-qc] .\n- [36] E. F. Boero and O. M. Moreschi, Monthly Notices of the Royal Astronomical Society 492 , 3763 (2019).\n- [37] E. Newman and R. Penrose, Journal of Mathematical Physics 3 , 566 (1962).\n- [38] S. Pineault and R. C. Roeder, Astrophys. J. 212 , 541 (1977).\n- [39] S. Seitz, P. Schneider, and J. Ehlers, Classical and Quantum Gravity 11 , 2345 (1994).\n- [40] E. Gallo and O. M. Moreschi, Phys. Rev. D 83 , 083007 (2011).\n- [41] C. C. Dyer, Monthly Notices of the Royal Astronomical Society 180 , 231 (1977).\n- [42] D. W. Berry and P. C. S. Costa, Quantum 8 , 1369 (2024).\n- [43] S. Frittelli, T. P. Kling, and E. T. Newman, Phys. Rev. D 63 , 023007 (2000).\n- [44] J. O. Shipley, Strong-field gravitational lensing by black holes (2019), arXiv:1909.04691 [gr-qc] .\n- [45] P. Dahal, Eur. Phys. J. Plus 138 , 205 (2023).\n- [46] A. I. Harte, General Relativity and Gravitation 51 , 10.1007/s10714-018-2494-x (2019).\n- [47] A. I. Harte, General Relativity and Gravitation 51 , 10.1007/s10714-019-2646-7 (2019).\n- [48] G. Cusin and M. Lagos, Phys. Rev. D 101 , 044041 (2020).\n- [49] C. Dalang, G. Cusin, and M. Lagos, Phys. Rev. D 105 , 024005 (2022).\n- [50] Z. Li, J. Qiao, W. Zhao, and X. Er, Journal of Cosmology and Astroparticle Physics 2022 (10), 095.\n- [51] K.-i. Kubota, S. Arai, and S. Mukohyama, Phys. Rev. D 109 , 044027 (2024).\n- [52] C. R. Keeton, (2002), arXiv:astro-ph/0102341 [astro-ph] .\n- [53] The superscript n represents the components along n -axis. It is similar for the m and ¯ m superscript in Eq. (11) and (12)."}