Evaluating the KirchhoffFresnel diffraction integral is essential in studying wave effects in astrophysical lensing but is often intractable because of the highly oscillatory integrand. A recent breakthrough was made by exploiting the PicardLefschetz theory the integral can be performed along the Lefschetz thimbles in the complex domain where the integrand is not oscillatory but rapidly converging. The application of this method however has been limited by both the unfamiliar concepts involved and the low numerical efficiency of the method used to find the Lefschetz thimbles. In this paper we give simple examples of the Lefschetz thimbles and define the flow lines that facilitate the understanding of the concepts. Based on this we propose new ways to obtain the Lefschetz thimbles with high numerical efficiency which provide an effective tool for studying wave effects in astrophysical lensing.
['Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
Acquiring the Lefschetz thimbles efficient evaluation of the diffraction integral for lensing in wave optics
{'Acquiring the Lefschetz thimbles: e ffi cient evaluation of the di ff raction integral for lensing in wave optics': 'Xun Shi 1 ⋆  \n1 South-Western Institute for Astronomy Research (SWIFAR), Yunnan University, 650500 Kunming, P. R. China  \n23 September 2024', 'ABSTRACT': "Evaluating the Kirchho ff -Fresnel di ff raction integral is essential in studying wave e ff ects in astrophysical lensing, but is often intractable because of the highly oscillatory integrand. A recent breakthrough was made by exploiting the Picard-Lefschetz theory: the integral can be performed along the 'Lefschetz thimbles' in the complex domain where the integrand is not oscillatory but rapidly converging. The application of this method, however, has been limited by both the unfamiliar concepts involved and the low numerical e ffi ciency of the method used to find the Lefschetz thimbles. In this paper, we give simple examples of the Lefschetz thimbles and define the 'flow lines' that facilitate the understanding of the concepts. Based on this, we propose new ways to obtain the Lefschetz thimbles with high numerical e ffi ciency, which provide an e ff ective tool for studying wave e ff ects in astrophysical lensing.  \nKey words: methods: numerical - software: development - gravitational lensing: strong - gravitational lensing: micro", '1 INTRODUCTION': "Wave optics description of electromagnetic wave propagation after a phase screen is given by the Kirchho ff -Fresnel di ff raction integral (e.g. Schneider et al. 1992; Born & Wolf 2019). How to e ffi -ciently evaluate the Kirchho ff -Fresnel di ff raction integral is thus the key problem in understanding di ff raction and interference phenomena that arise, e.g., in the case of gravitational or plasma lensing of coherent radio sources such as pulsars (Hewish et al. 1968), fast radio bursts (FRBs, Lorimer et al. 2007; Thornton et al. 2013), and gravitational waves (Abbott et al. 2016; Agazie et al. 2023; EPTA Collaboration et al. 2023; Reardon et al. 2023; Xu et al. 2023). Especially, scintillation and plasma lensing of radio pulsars and FRBs in the interstellar medium (e.g. Rickett 1990; Stinebring et al. 2001; Walker et al. 2004; Cordes et al. 2006; Brisken et al. 2010; Cordes et al. 2017) often involve structures comparable or smaller than the Fresnel scale. Thus, a particular interest in evaluating the Kirchho ff -Fresnel di ff raction integral exists in the e ff orts of understanding scintillation and plasma lensing caused by interstellar structures (Walker et al. 2004; Grillo & Cordes 2018; Jow et al. 2020, 2021; Shi & Xu 2021; Jow & Pen 2022; Jow et al. 2023; Shi 2024a,b).  \nThe Kirchho ff -Fresnel di ff raction integral is, however, non-trivial to calculate. Analytic solutions exist only for a few idealized situations (e.g. for point-mass lenses, see Deguchi & Watson 1986; Jow et al. 2020; Leung et al. 2023). The eikonal approximation, i.e. approximating the wave field as discrete stationary phase points and computing the interference among them, is precise at the highfrequency limit but does not capture the full wave e ff ects at finite frequencies. Due to the oscillatory nature of the integrand, numerical computation is in general challenging. Fast Fourier Transform (FFT) is in principle applicable for computing the di ff raction integral  \nnumerically for any lens shape. However, its application is usually limited to the perturbative wave optics regime where neither the frequency nor the fluctuations in the phase screen are high, otherwise, an exceedingly large number of sampling grids would be necessary (see e.g appendix A of Grillo & Cordes 2018). Between the regimes of validity of the eikonal approximation and the perturbative wave optics, there is a regime of non-trivial wave e ff ects (Shi 2024b). No viable method exists for this regime until the recent introduction of the Picard-Lefschetz theory (see Witten 2010) to radio astronomy by Feldbrugge et al. (2019, 2023).  \nThe Picard-Lefschetz theory is an exact and versatile approach for dealing with multidimensional oscillatory integrals which has recently been applied to fields ranging from condensed matter systems to quantum cosmology (e.g. Tanizaki et al. 2016; Feldbrugge et al. 2017; Mou et al. 2019). Its essence is to use Cauchy's theorem to deform the interval of the integration into the complex plane. Adding the imaginary dimension gives the freedom to choose the integration contour on which the integrand is non-oscillatory and the integral converges the most rapidly. The resulting contours are called the 'Lefschetz thimbles', which are a sum of the steepest descent contours connecting the saddle points of the integrand generalized to the complex domain. Integration along these thimbles is computationally e ffi cient and insensitive to numerical cuto ff s, making it suitable for practical applications.  \nDespite this success, the Picard-Lefschetz method has not yet been widely used. The main reason, apart from its conceptual difficulty, could be the limited numerical e ffi ciency in finding the Lefschetz thimbles for a general lensing configuration. Although performing the integration along the Lefschetz thimbles is numerically trivial, finding the Lefschetz thimbles is not. Feldbrugge et al. (2019) developed and advocated a numerical method called 'flowing the integration domain'. It progressively distorts the integration domain into the complex plane using a di ff erential equation. This method  \nhas been adopted by almost all subsequent astrophysics works using the Picard-Lefschetz method (Feldbrugge & Turok 2020; Feldbrugge 2023; Shi & Xu 2021; Jow et al. 2021; Jow & Pen 2022; Suvorov 2022; Jow et al. 2023, 2024). In particular, a numerical algorithm to realize this flow method for N -dimensional integrations was detailed in Jow et al. (2021). However, as we shall demonstrate in this paper, an infinite number of flow steps is required for convergence to the Lefschetz thimbles for almost all points on the original integration domain. In this 'flow method', the entire Lefschetz thimbles are contributed by a tiny fraction of the original integration domain, namely the vicinity of a number of discrete points. Thus, a lot of refinement of these small intervals is necessary during the flow process, leading to a low numerical e ffi ciency. With this understood, the numerical method of finding the Lefschetz thimbles can be immediately improved upon.  \nThe purpose of this paper is threefold: First, we would like to promote the familiarity with the Lefschetz thimbles by giving simple examples where their shapes are easy to understand. These examples can also serve as test cases for verifying numerical convergence. Second, we would like to aid the understanding of Lefschetz thimbles by defining the 'flow lines', which will help us show that each piece of Lefschetz thimble is flown from with an infinitesimal interval on the real axis. Third, we describe ways to e ffi ciently obtain the Lefschetz thimbles numerically.", '2 THE DIFFRACTION INTEGRAL AND THE LEFSCHETZ THIMBLES': "The wave amplitude received by an observer for a point source with a unit magnitude is given by the Kirchho ff -Fresnel di ff raction integral (e.g. Schneider et al. 1992; Born & Wolf 2019, see also Kramer et al. 2024 for a di ff erent approach)  \nE ( β ) = GLYPH<18> ν 2 π i GLYPH<19> Z exp GLYPH<2> i νϕ GLYPH<3> d 2 x . (1)  \nThis wave amplitude E ( β ) is also referred to as the 'transmission factor' of an intermediate phase screen acting as the lens (Schneider et al. 1992). It is contributed by all possible wave paths connecting the source, the observer, and a point x on the phase screen according to the Principle of Huygens. Here, β is the location of the source line-of-sight on the phase screen. The coordinates x and β are made dimensionless by scaling with the size of the lens a lens following the convention of Shi & Xu (2021) and Jow et al. (2023). The dimensionless parameter ν is addressed as the 'reduced frequency' since it is proportional to the observing frequency. In a lensing system, it is given by  \nν = 2 π (1 + z lens) a 2 lens λ ¯ D , (2)  \nwhere λ is the wavelength, and ¯ D = D l D ls / D s is a combination of the distances to the lens D l, the source D s, and that between the source and lens D ls. On cosmological scales, these distances should be interpreted as angular diameter distances, and one needs to include the factor with the lens redshift z lens to account for the redshift of the light. From the above definition, one can derive that the reduced frequency is the square of the lens size - Fresnel scale ratio, ν = a 2 lens / r 2 F . It is the key parameter that determines the importance of the wave e ff ects (see e.g. Shi 2024b).  \nThe phase function  \nϕ ( x ) = ( x -β ) 2 2 -κψ ( x ) (3)  \nFigure 1. Lens shapes used in this paper: a rational lens (magenta line) and a k max = 4 simG lens (blue solid line). The Gaussian lens (black dashed line) and the k max = 2 simG lens (green dotted line) are shown for comparison. As k max increases, the simG lens rapidly approaches the Gaussian lens.  \n<!-- image -->  \nincludes the geometrical phase delay and the dispersive phase delay from a lens potential with an amplitude κ and a shape ψ ( x ) (see Fig. 1 for shapes used in this paper). We have adopted the convention that a converging lens (e.g. a gravitational lens) has a negative κ value. Positive κ values correspond to diverging lenses, which can represent overdense plasma lenses (e.g. Clegg et al. 1998; Pen & King 2012; Er & Mao 2014; Cordes et al. 2017; Er & Rogers 2018, 2019; Wagner & Er 2020; Shi & Xu 2021).  \nFrom Eqs. 1-3, one can see that how oscillatory the integrand is depends on the reduced frequency ν and the dimensionless lens amplitude κ . For lensing systems with relatively large ν or κ values e.g. plasma lensing of pulsars by the structures in the interstellar medium (e.g. Hill et al. 2005; Brisken et al. 2010; Sprenger et al. 2022), the integrand is highly oscillatory and the di ff raction integral is intractable with traditional numerical methods (Grillo & Cordes 2018).  \nThe Picard-Lefschetz theory has been recently introduced to evaluate di ff raction integrals in astrophysical lensing by Feldbrugge et al. (2019, 2023). It provides an exact approach to computing highly oscillatory integrals for a wide range of lenses, and is highly complementary to traditional numerical methods such as FFT.  \nBy generalizing the integrand to the complex domain x → z = x + i y and separating the real and imaginary parts of the phase function,  \ni ϕ ≡ h + i H , (4)  \nthe Picard-Lefschetz theory deforms the integration contour into a sum of the steepest descent contours of the h -function where the integrand is non-oscillatory and rapidly converging. The resulting contours are called the 'Lefschetz thimbles' - a sum of smooth subcontours connecting the saddle points of the integrand.  \nEach piece of Lefschetz thimble J i is associated with one image zi . By its definition, J i is the steepest descent contours connected to the image zi . It is also all of the following:  \n- (i) fix point of the downward flow of the h -function;\n- (ii) pathline of the downward flow starting from the image zi ;\n- (iii) a subset (with descending h ) of the constant phase contours connected to the image zi .  \nSee Fig. 2 for an illustration of the Lefschetz thimbles.  \nIn addition to the Lefschetz thimbles associated with real saddle points i.e. the geometric images, there will generally be contributions from Lefschetz thimbles associated with saddle points with  \nnon-negligible imaginary parts i.e. the imaginary images (Jow et al. 2021; Shi 2024a), e.g., one marked by the black cross at Im(x) ≈ -1 . 5 in Fig. 2.", '3 LEFSCHETZ THIMBLE SHAPE EXAMPLES': 'Here we show examples of the Lefschetz thimbles for a few simple cases. For simplicity, we shall limit our discussion to 1D lenses in this and the next sections, and discuss a generalization to 2D situations in Section 6.', '3.1 Free propagation': 'We first use the simplest case - image in the absence of any lens to illustrate how an image is linked to the Lefschetz thimble. In this case, the phase  \ni ϕ = i ( x + i y -β ) 2 2 = -( x -β ) y + i ( x -β ) 2 -y 2 2 . (5)  \nThere exists a unique real image at x 0 = β , y 0 = 0, with h 0 = H 0 = 0.  \nThe Lefschetz thimble is the contour with constant imaginary phase H = H 0 and h < h 0 across the image. The single line satisfying these conditions is y = x -β with slope k = 1. Along this Lefschetz thimble, the integral becomes  \nE ( β ) = (1 + i) r ν 2 π i Z exp GLYPH<16> -ν x 2 GLYPH<17> d x . (6)  \nThe integrand is a Gaussian distribution with width σ lens / √ 2 ν = r F / √ 2. This matches the common wisdom that the Fresnel scale characterizes the size of the image in free propagation.', '3.2 In the vicinity of a lensed image': 'Next, let us consider the case of a lens that creates additional lensed images. In the existence of a lens, the shape of the Lefschetz thimbles still asymptotes to y = x -β far away from the lens at x →±∞ .  \nWe now seek the expression for the slope k of the thimble near an image in the complex plane zi . Since the image is a stationary phase point with ϕ \' ( zi ) = 0, the phase can be expanded around the image as  \nϕ ≈ ϕ ( zi ) + ϕ \'\' ( zi ) δ z 2 2 . (7)  \nLetting ϕ \'\' ( zi ) = ϕ \'\' real , i + i ϕ \'\' imag , i , and writing δ z = (1 + i k ) δ x , we get  \nδ h = Real i ϕ \'\' (zi) δ z 2 2 ! = " -ϕ \'\' real , i k -ϕ \'\' imag , i 1 -k 2 2 # δ x 2 , (8)  \nδ H = Imag i ϕ \'\' (zi) δ z 2 2 ! = " ϕ \'\' real , i 1 -k 2 2 -ϕ \'\' imag , i k # δ x 2 . (9)  \nWe can then obtain the value of k from the condition that δ H = 0 along a Lefschetz thimble.  \nFor a real image, ϕ \'\' imag , i = 0 and thus k 2 = 1. From the condition δ h = -ϕ \'\' real , i k δ x 2 < 0, one can derive k = 1 for an image with the same parity with the source which has a magnification parameter µ i ≡ 1 /ϕ \'\' real , i > 0, and k = -1 for an image with the opposite parity and µ i < 0.  \nFigure 2. Flow lines (white contours) for points on the real axis for a lensing system with a lens ψ ( x ) = 1 / (1 + x 2 ) with an amplitude κ = 10 placed at the origin, and the source at β = 5. The images are marked with black crosses, and the poles are marked with orange dots. The Lefschetz thimble is shown as the red contour. The contours of the H -function are shown in the background.  \n<!-- image -->  \nFor an imaginary image with ϕ \'\' imag , i , 0, δ H = 0 gives  \nk = -ϕ \'\' imag , i ϕ \'\' real , i ± v t       ϕ \'\' imag , i ϕ \'\' real , i       2 + 1 , (10)  \ni.e., the slope depends on the shape of the phase function around the image in the complex domain and has no bound.', '4 THE FLOW LINES': "The works of Feldbrugge et al. (2019, 2023) have established a link between the original integration contour - the real axis and the contour of Lefschetz thimbles: the former can 'flow' into the latter with a downward flow γℓ of the h -function defined by  \n∂γℓ ( z ) = -∇ zh [ γℓ ( z )] , (11)  \n∂ℓ  \ni.e., the flowed integration contour X ℓ will converge to a set of steepest descent contours i.e. the Lefschetz thimbles as the number of flow steps ℓ →∞ ,  \nlim →∞ X ℓ = J . (12)  \nℓ  \nHowever, to consider the map between these two contours, one still needs to know how intervals on these contours map to each other. For this purpose, we define 'flow lines' for points on the real axis, with a flow line tracking the trajectory of a point as it flows. In other words, the flow line of a point is the pathline of the downward flow of the h -function starting from this point.  \nFig. 2 presents the flow lines for points on the real axis for a lensing system with a lens of amplitude κ = 10 at β = 5. The corresponding integral is R ∞ -∞ exp h i ν h ( x -5) 2 / 2 -10 / (1 + x 2 ) ii d x . One  \ncan see that all flow lines end at either infinity or a pole, and a piece of Lefschetz thimble J i connects either an infinity and a pole, two poles, or two infinities.  \nThe flow lines (white contours) help define the optimal generalized integration contour for an interval ( x 0 , x 1) on the original integration axis. When the flow lines associated with both x 0 and x 1 end at the same infinity, e.g. in the case of x 0 = -10 and x 1 = -9 in Fig. 2, the integration over ( x 0 , x 1) can be performed e ffi ciently along the flow line of x 0 from x 0 to infinity and then again from infinity to x 1 along the flow line of x 1.  \nRemarkably, a whole piece of Lefschetz thimble J i corresponds only to an infinitesimal interval on the real axis. The location of this interval ai is where the real axis crosses the steepest ascending contour of the image. For a finite interval that contains ai , the optimal generalized integration contour then consists of J i in addition to the flow lines of the two endpoints of the interval. The flow lines for the special points ai are the most special. They flow first onto an image location zi , and from there bifurcate into flows along the thimble piece J i in two directions reaching infinities and / or poles.", '5 FINDING THE LEFSCHETZ THIMBLES NUMERICALLY': 'That the Lefschetz thimble J i is on the flow line of one point on the real axis has implications for the numerical method of finding J i .  \nThe numerical method favored by Feldbrugge et al. (2019, 2023) and inherited by latter studies of the wave e ff ects (Feldbrugge & Turok 2020; Feldbrugge 2023; Shi & Xu 2021; Suvorov 2022; Jow et al. 2021; Jow & Pen 2022; Jow et al. 2023, 2024) is to flow the integrand using Eq. 11. However, the flow lines demonstrate that the convergence of this flow to the Lefschetz thimbles (Eq. 12) occurs only at ℓ → ∞ . Most points sampled on the real axis will not contribute to the final Lefschetz thimble at all. They only asymptote to J at ±∞ where J approaches the free propagation solution. Numerically, to find a Lefschetz thimble piece J i one needs to repeatedly refine intervals in the vicinity of the flow line of ai . A lot of computation is spent on flowing the other intervals which are finally removed because their contributions turn out to fall below a threshold. This imposes a limitation on the numerical e ffi ciency of this flow method in finding the Lefschetz thimbles.  \nWith the current flow method, it is particularly di ffi cult for the integration contour to converge to the Lefschetz thimbles in regions far away from the real axis in the complex domain. However, the precise shape of the Lefschetz thimbles in these regions is increasingly important for the evaluation of the di ff raction integral at smaller reduced frequency ν , and such a shape is non-trivial for a lensing system with a lens well-separated from the source (i.e. with | β | ≫ 1). Most computations currently performed with the Picard-Lefschetz method avoid this problem by limiting to relatively high ν and small | β | values (i.e. with small angular o ff sets between the source and the lens). In such a regime, most line intervals quickly dip below a threshold in their contribution to the integral and are removed from the flow to improve numerical e ffi ciency. Although the small | β | , relatively high ν regime is relevant in many applications, it is not the only regime of interest. In the study of pulsar scintillation, wave effects in the large | β | regime are of particular interest. This regime is challenging for the current flow method, and the typical large lens amplitude κ for pulsar scintillation (e.g. Shi 2021) makes the perturbative methods invalid already at small reduced frequencies (Jow  \net al. 2023). A more e ffi cient approach to finding the Lefschetz thimbles is particularly needed for this regime.', '5.1 The improved flow method': 'To see how one can improve on this original method, let us revisit the di ff erent conceptual perspectives of a Lefschetz thimble piece J i (see Sect. 2). The original flow method utilizes the first perspective (i), that J i is the fixed point of the downward flow of the h -function. The second perspective (ii), that J i is the pathline of the downward flow starting from zi , gives a straightforward way to improve on it: Instead of starting the flow from sampled points on the whole real axis, one just needs to start from points in the vicinity of each effective image. The pathlines i.e. flow trajectories of these points will yield the Lefschetz thimbles J . Clearly, this improved flow method can save a lot of computing power and memory compared to the original flow method.  \nOne advantage of the original flow method is that it does not require an identification of the e ff ective images. The improved flow method, however, does require this identification. The identification can be achieved by first flowing points in the vicinity of all stationary phase points with the upward flow of the h -function to trace the ascending contours, and selecting the e ff ective images as those with an ascending contour that crosses the real axis. Considering this additional step, the improved flow method could still be more e ffi cient than the original flow method for some lensing configurations.', '5.2 The constant phase contour method': 'The third perspective (iii), that J i is a subset (with descending h ) of the constant phase contours connected to the image zi , provides a completely di ff erent way to find the Lefschetz thimbles. Taking advantage of the fact that the phase function is non-oscillatory on a piece of Lefschetz thimble J i , one can find J i starting from globally mapping the constant phase contour H J i = Hi . This results in a very e ffi cient method to find the Lefschetz thimbles which we would like to promote in this paper and have applied to recent works (Shi 2024a,b). The algorithm is sketched below (see also a demonstration in Fig. 3):  \n- (1) get the saddle points of the phase function i.e. the images;\n- (2) get the relevant Hi values at the images;\n- (3) get the contour H J i = Hi associated with each Hi value;\n- (4) cut the contour at images, poles, and integration domain boundaries to obtain 2 descending (with h < hi ) + 2 ascending (with h > hi ) contour pieces linked to each image;\n- (5) select the e ff ective images that are either real or with an ascending contour passing the real axis;\n- (6) connect descending contours of all e ff ective images to form the overall Lefschetz thimbles J .  \nUnlike the flow method which is a dynamic method that utilizes predominantly the h -function, this constant phase contour method is a static method that relies on the global information of the H -function. The idea of this method has been mentioned in Feldbrugge et al. (2019) and Feldbrugge et al. (2023), but they thought it would have low numerical e ffi ciency and disfavored it compared to the flow method. However, given that fast contour-finding algorithms are available in most computational languages, the numerical e ffi -ciency of this method can be extremely high. With our sample script written in Python run on a laptop, it takes less than 2 seconds for the algorithm to locate the Lefschetz thimbles for the lensing configuration presented in Figure A2 of Jow et al. (2023). The numerical  \nFigure 3. An example of optimal integration contour found using the constant phase contour method for the radial integration of a 2D integration (see Eq. 13). The source at β x = 0, β y = 5 is lensed by a converging 2D approximate Gaussian lens (see Fig. 1) with an amplitude κ = -10. The polar coordinate of the 2D integration is taken to be θ = π / 6 for this example. With the constant phase contour method, one first locates all stationary phase points (white crosses), finds contours with H J i = Hi (dashed lines), and breaks them into steepest ascending and descending contours. Then, one identifies the e ff ective images from the fact that one of their associated steepest ascending contours crosses the real axis. Finally, one connects the steepest descending contours of the e ff ective images to form the Lefschetz thimbles (the blue solid line). The optimal integration contour for the radial integration (black dotted line) connects the r = 0 point to the Lefschetz thimble with the flow line for r = 0. Contours of the h -function are shown in the background.  \n<!-- image -->  \ne ffi ciency of this method is not sensitive to the range of integration intervals, and thus is particularly useful for lensing configurations with large source-lens angular o ff set | β | and / or small reduced frequency ν . Another advantage of this new approach is that the precision of the resulting Lefschetz thimbles is set by the resolution of the grid for contour-finding, i.e. the precision is known and easily controllable.  \nNote that the constant phase domain has a dimension of 2 N -1 for an N-dimensional integral while the Lefschetz thimbles have a dimension of N . Their dimensions match only for one-dimensional integrals. We show how this method can be developed into an e ffi -cient method for a general 2D integral in the following.', '6 GENERALIZATION TO 2D INTEGRALS': "A great advantage of the original flow method is that it can be easily applied to oscillatory integrals of any dimension. The methods presented in this paper both require all stationary phase points in the generalized complex domain to be found as a first step. This is numerically trivial for a 1D integral but is much harder for a higher dimensional integral with a general, asymmetric lens.  \nHowever, as noticed by Feldbrugge & Turok (2020) and Tambalo et al. (2023), there is an easy way out for 2D integrals. When written  \nin the polar coordinate,  \nE ( β ) = ν 2 π i Z 2 π 0 d θ Z ∞ 0 r e i νϕ ( r ,θ, β ) d r , (13)  \nthe integral in θ is over a finite range and can be performed with standard numerical techniques 1 . Thus, the problem reduces to an evaluation of a 1D oscillatory integral in r . The only caveat is that the lower limit of integration is 0 rather than -∞ . One just needs to connect the Lefschetz thimbles with the flow line for r = 0 to form a closed integration contour (see Fig. 3).  \nThis solves the general problem of studying single-screen astrophysical lensing in the wave optics framework since 2D integrals provide an adequate description there. We show an application of the constant phase contour method to compute the di ff raction pattern produced by a 2D approximate Gaussian lens in the appendix.  \nFinally, we note the restriction of lens shapes that can be treated with the Picard-Lefschetz method. Since an analytic continuation into the complex plane is a key procedure in the Picard-Lefschetz method, the lens shape must be such that the phase function can be analytically continued into the complex plane. For lensing by point masses, the phase involves a logarithm function which makes the generalization into the complex plane more intricate. However, the  \nlogarithm potential represents a special case where the result of the Kirchho ff -Fresnel di ff raction integral can be expressed analytically in terms of a special function (Nakamura & Deguchi 1999). The Gaussian lens shape ψ ∝ e -x 2 also poses a di ffi culty by having an infinite number of saddle points in the complex plane and an essential singularity at infinity. One can, nevertheless, approximate the Gaussian lens using a rational function ψ ( x ) ∝ 1 / P k max k = 0 2 -k x 2 k k ! to arbitrarily high precision (this lens shape is referred to as the 'simG' lens, see Fig. 1). Similar approximations can be made for any lens shape.", '7 CONCLUSION': "The recently introduced Picard-Lefschetz method has enabled the evaluation of the oscillatory Kirchho ff -Fresnel di ff raction integral for a general lensing system. This paper is meant to clear the main obstacles to its application by clarifying the concept of the Lefschetz thimbles (the new integration contours in the complex domain) and proposing new e ffi cient ways to find them.  \nWe show that in the case of free propagation, the Lefschetz thimbles in the complex plane z = x + i y generated from each integration axis is a linear line with slope one: y = x -β . In the presence of a lens, the Lefschetz thimbles are still asymptotic to this line far away from the lens. In the vicinity of an image, the slope of the Lefschetz thimble is k = ± 1 around real images and can have an arbitrary value determined by the shape of the phase function around imaginary images.  \nThe 'flow lines' are defined to show the pathlines of the downward flow of the h -function (real part of the phase i ϕ ). With these pathlines, it is clear that the Lefschetz thimbles are mapped from tiny intervals on the real axis, each corresponding to an image. Thus, the originally proposed 'flowing the integration domain' method for finding the Lefschetz thimbles is far from optimal. Although the flow method claims to let the integration domain converge to the Lefschetz thimbles, it requires an infinite number of flow steps and the convergence occurs only at ±∞ for almost all intervals on the original integration domain. This renders this method numerically ine ffi -cient and makes it di ffi cult to apply to lensing at very small reduced frequencies and / or lensing systems with large angular separations between the lens and the source.  \nWe propose new methods of finding the Lefschetz thimbles that either start the flow from points in the vicinity of the images or use the constant phase contour method. The latter method which is based on the global information of the H -function (complex part of i ϕ ), is the prime method for 1D integrals and is particularly e ffi cient and accurate. We have sketched the algorithm for this method and shown that it can also be developed into an e ffi cient method for 2D integrals.  \nWith these insights into the Lefschetz thimbles and the e ffi cient ways to locate them, the oscillatory nature no longer poses a challenge to the integration of the Kirchho ff -Fresnel di ff raction integral. This provides an e ff ective tool for the study of wave e ff ects for a wide range of lensing configurations.", 'DATA AVAILABILITY STATEMENTS': 'No new data were generated or analyzed in support of this research.', 'ACKNOWLEDGEMENTS': 'X.S. thanks the referees for the helpful comments and suggestions. This work is supported by NSFC No. 12373025.', 'REFERENCES': '- Tanizaki Y., Hidaka Y., Hayata T., 2016, New Journal of Physics, 18, 033002  \nThornton D., et al., 2013, Science, 341, 53 Wagner J., Er X., 2020, arXiv:2006.16263, p. arXiv:2006.16263 Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43 Witten E., 2010, arXiv:1009.6032, p. arXiv:1009.6032 Xu H., et al., 2023, Research in Astronomy and Astrophysics, 23, 075024', 'APPENDIX A: DIFFRACTION BY A 2D LENS': "As an application of the constant phase contour method to the evaluation of the Kirchho ff -Fresnel di ff raction integral, we compute the di ff raction pattern produced by a 2D lens. The lens is an axisymmetric, converging lens with an approximate Gaussian ('simG') shape with k max = 4 (see Fig. 1) and an amplitude κ = -10.  \nIn Fig. 1 we show the di ff raction pattern as a function of the 2D source location β x and β y for three di ff erent reduced frequencies. The upper panels are results computed with the Picard-Lefschetz method using the numerical approach outlined in Section 6. The lower panels are results computed with the eikonal approximation. The chosen frequencies are within the range where wave e ff ects are nontrivial (Shi 2024b). In this range, the results considering full wave e ff ects di ff er significantly both from their eikonal approximations and also from the unlensed cases. Notably, there exists a caustic at | β | = 5 . 1 outside which there exists only one real image, and thus the eikonal results show no flux variation outside the caustic. Within the caustic, the flux under the eikonal approximation is computed from the interference of the three real images created by the lens. The Picard-Lefschetz method, on the other hand, captures the full wave e ff ects and shows a clear di ff raction pattern both inside and outside the caustic that di ff ers from the eikonal results.  \nThe numerical speed for finding the Lefschetz thimbles with the constant phase contour method depends on the complexity of the lens. The required time increases with the number of stationary phase points. For the k max = 4 simG lens with 17 stationary phase points (Fig. 3), our sample code takes ∼ 10 seconds for each radial integral of Eq. 13 using one core. The required time for each pixel in the di ff raction pattern also depends on the sampling of the angular integral. The necessary θ sampling heavily depends on the reduced frequency ν . For a higher ν , the integrand over θ is more oscillatory. We have used 120 θ -points for ν = 3. This is close to the maximum number of θ -points needed since for ν ≳ 10 the eikonal approximation would work well and there is little need for full integration of the Kirchho ff -Fresnel di ff raction integral.  \nFigure 1. Di ff raction patterns computed using the numerical method outlined in Section 6 (upper panels). Shown are the magnifications of a point source ( | E | 2 ) at various source locations β x and β y after being lensed by a 2D approximate Gaussian lens (see Fig. 1) with an amplitude κ = -10. As a comparison, the di ff raction patterns computed using the eikonal approximation are shown in the lower panels. The reduced frequency ν is set to 0 . 3, 1, and 3 for the columns rows, respectively.  \n<!-- image -->"}