Acquiring the Lefschetz thimbles: efficient evaluation of the diffraction integral for lensing in wave optics

Evaluating the Kirchhoff-Fresnel diffraction integral is essential in studying wave effects in astrophysical lensing, but is often intractable because of the highly oscillatory integrated. 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 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.


INTRODUCTION
Wave optics description of electromagnetic wave propagation after a phase screen is given by the Kirchhoff-Fresnel diffraction integral (e.g.Schneider et al. 1992;Born & Wolf 2019).How to efficiently evaluate the Kirchhoff-Fresnel diffraction integral is thus the key problem in understanding diffraction 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 Kirchhoff-Fresnel diffraction integral exists in the efforts of understanding scintillation and plasma lensing caused by interstellar structures (Walker et al. 2004;Grillo & Cordes 2018;Jow et al. 2020Jow et al. , 2021;;Shi & Xu 2021;Jow & Pen 2022;Jow et al. 2023;Shi 2024a,b).
The Kirchhoff-Fresnel diffraction 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 effects 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 diffraction integral ⋆ E-mail: xun@ynu.edu.cnnumerically 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 effects (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. (2019Feldbrugge et al. ( , 2023)).
The 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 efficient and insensitive to numerical cutoffs, making it suitable for practical applications.
Despite 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 efficiency 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 differential equation.This method has 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. 2023Jow et al. , 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 efficiency.With this understood, the numerical method of finding the Lefschetz thimbles can be immediately improved upon.
The 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 efficiently obtain the Lefschetz thimbles numerically.

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 Kirchhoff-Fresnel diffraction integral (e.g.Schneider et al. 1992;Born & Wolf 2019, see also Kramer et al. 2024 for a different approach) This 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 where λ 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 effects (see e.g.Shi 2024b).
The phase function includes 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.From 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 diffraction integral is intractable with traditional numerical methods (Grillo & Cordes 2018).
The Picard-Lefschetz theory has been recently introduced to evaluate diffraction integrals in astrophysical lensing by Feldbrugge et al. (2019Feldbrugge et al. ( , 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.
By generalizing the integrand to the complex domain x → z = x+ iy and separating the real and imaginary parts of the phase function, the 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.Each piece of Lefschetz thimble J i is associated with one image z i .By its definition, J i is the steepest descent contours connected to the image z i .It is also all of the following: (i) fix point of the downward flow of the h-function; (ii) pathline of the downward flow starting from the image z i ; (iii) a subset (with descending h) of the constant phase contours connected to the image z i .
See Fig. 2 for an illustration of the Lefschetz thimbles.
In 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 non-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.

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.

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 There exists a unique real image at x 0 = β, y 0 = 0, with h 0 = H 0 = 0.The 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 The integrand is a Gaussian distribution with width σ lens / 2. This matches the common wisdom that the Fresnel scale characterizes the size of the image in free propagation.

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 → ±∞.
We now seek the expression for the slope k of the thimble near an image in the complex plane z i .Since the image is a stationary phase point with ϕ ′ (z i ) = 0, the phase can be expanded around the image as , and writing δz = (1 + ik)δx, we get We can then obtain the value of k from the condition that δH = 0 along a Lefschetz thimble.For 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. For an imaginary image with ϕ ′′ imag,i 0, δH = 0 gives i.e., the slope depends on the shape of the phase function around the image in the complex domain and has no bound.

THE FLOW LINES
The works of Feldbrugge et al. (2019Feldbrugge et al. ( , 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 i.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 ℓ → ∞, However, 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.
Fig. 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 can 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.
The 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 efficiently 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 .
Remarkably, a whole piece of Lefschetz thimble J i corresponds only to an infinitesimal interval on the real axis.The location of this interval a i is where the real axis crosses the steepest ascending contour of the image.For a finite interval that contains a i , 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 a i are the most special.They flow first onto an image location z i , and from there bifurcate into flows along the thimble piece J i in two directions reaching infinities and/or poles.

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 .
The numerical method favored by Feldbrugge et al. (2019Feldbrugge et al. ( , 2023) and inherited by latter studies of the wave effects (Feldbrugge & Turok 2020;Feldbrugge 2023;Shi & Xu 2021;Suvorov 2022;Jow et al. 2021;Jow & Pen 2022;Jow et al. 2023Jow et al. , 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 a i .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 efficiency of this flow method in finding the Lefschetz thimbles.
With the current flow method, it is particularly difficult 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 diffraction 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 offsets 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 efficiency.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 et al. 2023).A more efficient approach to finding the Lefschetz thimbles is particularly needed for this regime.

The improved flow method
To see how one can improve on this original method, let us revisit the different 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 z i , 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.
One advantage of the original flow method is that it does not require an identification of the effective 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 effective images as those with an ascending contour that crosses the real axis.Considering this additional step, the improved flow method could still be more efficient than the original flow method for some lensing configurations.

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 z i , provides a completely different 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 = H i .This results in a very efficient 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): (1) get the saddle points of the phase function i.e. the images; (2) get the relevant H i values at the images; (3) get the contour H J i = H i associated with each H i value; (4) cut the contour at images, poles, and integration domain boundaries to obtain 2 descending (with h < h i ) + 2 ascending (with h > h i ) contour pieces linked to each image; (5) select the effective images that are either real or with an ascending contour passing the real axis; (6) connect descending contours of all effective images to form the overall Lefschetz thimbles J.
Unlike 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 Hfunction.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 efficiency and disfavored it compared to the flow method.However, given that fast contour-finding algorithms are available in most computational languages, the numerical efficiency 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 Figure 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 = H i (dashed lines), and breaks them into steepest ascending and descending contours.Then, one identifies the effective 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 effective 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.
efficiency 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 offset |β| 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.
Note that the constant phase domain has a dimension of 2N − 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 efficient method for a general 2D integral in the following.

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.
However, as noticed by Feldbrugge & Turok (2020) and Tambalo et al. (2023), there is an easy way out for 2D integrals.When written in the polar coordinate, r,θ,β) dr , the integral in θ is over a finite range and can be performed with standard numerical techniques1 .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).This 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 diffraction pattern produced by a 2D approximate Gaussian lens in the appendix.
Finally, 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 logarithm potential represents a special case where the result of the Kirchhoff-Fresnel diffraction integral can be expressed analytically in terms of a special function (Nakamura & Deguchi 1999).The Gaussian lens shape ψ ∝ e −x 2 also poses a difficulty 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 ψ 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.

CONCLUSION
The recently introduced Picard-Lefschetz method has enabled the evaluation of the oscillatory Kirchhoff-Fresnel diffraction 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 efficient ways to find them.
We show that in the case of free propagation, the Lefschetz thimbles in the complex plane z = x + iy 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.
The '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 inefficient and makes it difficult to apply to lensing at very small reduced frequencies and/or lensing systems with large angular separations between the lens and the source.
We 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 efficient and accurate.We have sketched the algorithm for this method and shown that it can also be developed into an efficient method for 2D integrals.
With these insights into the Lefschetz thimbles and the efficient ways to locate them, the oscillatory nature no longer poses a challenge to the integration of the Kirchhoff-Fresnel diffraction integral.This provides an effective tool for the study of wave effects for a wide range of lensing configurations. .Diffraction 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 diffraction 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.

Figure 1 .
Figure 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.

Figure 2 .
Figure 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.

Figure 1
Figure1.Diffraction 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 diffraction 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.