Lensing of gravitational waves: efficient wave-optics methods and validation with symmetric lenses

Gravitational wave (GW) astronomy offers the potential to probe the wave-optics regime of gravitational lensing. Wave optics (WO) effects are relevant at low frequencies, when the wavelength is comparable to the characteristic lensing time delay multiplied by the speed of light, and are thus often negligible for electromagnetic signals. Accurate predictions require computing the conditionally convergent diffraction integral, which involves highly oscillatory integrands and is numerically difficult. We develop and implement several methods to compute lensing predictions in the WO regime valid for general gravitational lenses. First, we derive approximations for high and low frequencies, obtaining explicit expressions for several analytic lens models. Next, we discuss two numerical methods suitable in the intermediate frequency range: 1) Regularized contour flow yields accurate answers in a fraction of a second for a broad range of frequencies. 2) Complex deformation is slower, but requires no knowledge of solutions to the geometric lens equation. Both methods are independent and complement each other. We verify sub-percent accuracy for several lens models, which should be sufficient for applications to GW astronomy in the near future. Apart from modelling lensed GWs, our method will also be applicable to the study of plasma lensing of radio waves and tests of gravity.

Gravitational lensing, the deflection of waves by gravitational fields, has become an essential tool for exploring the Universe's structure and contents. The rich phenomenology of gravitational lensing [1] has enabled many applications, from inferring cosmological parameters to testing dark matter models. Progress on these fronts has relied exclusively on observations across the electromagnetic spectrum. However, recent advances in GW astronomy [2,3] open up a new arena for gravitational lensing. In time, searches of lensed GWs [4,5] are likely to turn into conclusive detections and novel applications.
Differences between gravitational and electromagnetic radiation from astrophysical sources make GW lensing qualitatively distinct and complementary to electromagnetic observations. GWs emit coherently and at much lower frequencies. GWs detectable by LIGO-Virgo-Kagra (LVK) have wavelengths more than three orders of magnitude longer than the lowest frequency radio waves permitted by the Earth's ionosphere. This difference may allow the observation of WO effects [6,7], which emerge when the wavelength is comparable to the time delay produced by the lens multiplied by the speed of light. WO effects are frequency dependent and their detection would allow an accurate determination of the lens properties [8][9][10][11][12]. Moreover, WO lensing of GWs could serve to identify stellar-scale microlenses [13][14][15][16][17] and test dark matter scenarios [6,[18][19][20][21][22][23][24].
Unfortunately, accurately computing lensed waveforms in the WO regime is challenging. A closed-form analytical expression exists only for the simplest point-mass lens [25], and series expansions have been developed for more general lenses [26]. These solutions have been used widely to study WO lensing. However, even these simple expressions are costly to evaluate in practice, particularly at high frequencies. General predictions require conditionally convergent integrals of rapidly oscillating functions over the lens plane. Previous works used direct integration [6,27], Levin's algorithm method [28,29], sampling the Fermat potential over contours [16,30] (related to our first method) or discretely [14,15,17], discrete FFTconvolution [31], and Picard-Lefschetz theory [32,33] (related to our second method). While these methods have been used to study complex lenses (e.g. Refs [14][15][16][17]), they have been validated (e.g. cross-validated with independent calculations) only for simple examples.
Here we describe methods to obtain WO predictions and cross-validate them on several lens models. In Section II, we present the diffraction integral. Section III presents expansions valid in the low-and high-frequency limits before turning to general, numerical algorithms. In Section IV we solve the Fourier transform of the integral by adaptively sampling contours of equal time delay. Then, in Section V we analytically continue the integration variable to make the integral manifestly convergent. Finally, in Section VI we validate the accuracy of both methods and discuss their performance. We will explore the phenomenology of GW lensing separately [34]. Throughout this paper, we will work in a unit system with c = 1.

II. WAVE OPTICS REGIME OF GRAVITATIONAL LENSING
In this Section, we will review the equations governing gravitational lensing in the WO regime. In order to focus on the mathematical problem we will not provide a detailed derivation of the quantities involved. Readers are referred to Refs. [8,34] for details. Our goal is to evaluate the diffraction integral, which we will give in dimensionless form: See Ref. [35] for a derivation. The integration is over the lens plane, with the coordinates rescaled by a dimensionful scale ξ 0 (e.g. a characteristic scale of the lens), so x is dimensionless. The impact parameter y is rescaled by η 0 ≡ D S ξ 0 /D L , where D S , D L are the angular diameter distances to the lens and the source, respectively.
Here we introduced the dimensionless frequency which is given in terms of a redshifted effective lens mass: The factor d eff ≡ D L D LS (1+z L )D S also depends on the angular diameter distance between the lens and the source D LS . For a point lens, M Lz coincides with the total mass of the lens (i.e. setting ξ 0 to be the Einstein radius), but this is not true for extended lenses.
The integral depends on the Fermat potential : Here ψ is the lensing potential, which depends on the matter distribution projected on the lens plane and whose derivative gives the deflection angle (Eq. (5) below). We conventionally shift by ϕ m (y), the global minimum value of the Fermat potential. From here on, we will suppress ϕ m (y) in our formulas and assume that it is added to make the minimum arrival time equal to zero. When necessary, we will introduce it back. We will consider several lens models in this work, summarized in Table I. All of them are spherically symmetric, leading to an axi-symmetric projected mass and lensing potential ψ(x) = ψ(x) (here and in the following x ≡ |x|). First, we consider the point lens, whose analytic solution will help us test the accuracy of different methods in Sec. VI. We will additionally consider three extended lenses: the Singular Isothermal Sphere (SIS), and two one-parameter extensions. SIS lenses follow from a matter profile ρ ∝ 1/r 2 and are often employed to model lensing by galaxies. Our first extension, the generalized-SIS (gSIS), has an arbitrary slope ρ ∝ 1/r 1+k (0 < k < 2) and can be used to model steeper or shallower lenses [35][36][37]. Its central density diverges, but the enclosed mass up to some radius remains finite if k < 2. Our second extension, the Cored Isothermal Sphere (CIS), has a central core of physical radius r c = x c ξ 0 [38,39]. Therefore, the matter density ρ ∝ 1/(r 2 + r 2 c ) is finite at the centre. Details about these lenses and their phenomenology will be provided separately [34]. Figure 1 shows the amplification factor for an SIS lens with impact parameter y = 0.3. The full WO solution was obtained by regularized contour flow, discussed in Sec. IV, matched to an analytic expansion at high frequencies. The remaining curves correspond to the expansions presented in Sec. III, each with limited range of validity: Geometric Optics (GO) and its next-order refinement (bGO) (Sec. III A) are good descriptions at high frequency, while the series expansion (Sec. III B) is a good approximation only for w ≪ 1. GO remains bounded at all frequencies, while bGO and the series expansion diverge at low/high frequencies, respectively.

III. ANALYTIC EXPANSIONS
We now present analytic expansions valid in the highand low-frequency limits, Subsections III A, III B, respectively.

A. Geometric optics & beyond
In the high-frequency limit, following the same arguments leading to the stationary-phase approximation for path integrals, only the neighbourhoods of extrema of the Fermat potential (4) contribute to the amplification factor (1). Each extremum is associated with an image J, located at a position x J in the image plane where the lens equation is satisfied (here α(x J ) ≡ ∇ x ψ(x J ) is the deflection angle and ∇ x is the gradient computed with respect to x). The geometric optics regime emerges from a quadratic expansion of the Fermat potential around each image so that the diffraction integral can be performed analytically.
The GO amplification factor (1) receives contributions from each image J where the magnification is evaluated on the image position x J (the second equality above applies to the specific case of axially-symmetric lenses). In the above expressions, ϕ J is the Fermat potential of the J-th image, ϕ ,ij ≡ ∂ i ∂ j ϕ is its Hessian matrix and α(x) ≡ |α(x)|. As we are working in the twodimensional lens plane, i, j, · · · ∈ {1, 2}, corresponding to the x 1 and x 2 coordinates (the Cartesian components of x). The Morse Phase [8,35,40] depends on the type of image as .
(8) Minima, saddle points and maxima of the time delay function are also known as type I, II and III images, respectively [41].

Beyond geometric optics
Beyond GO (bGO) corrections can be obtained as a series expansion in 1/w. We now review the leading bGO correction, following [27] and focusing on axiallysymmetric lenses (see also [42]). First of all, we expand the lensing potential ϕ around each image x J up to quar- For a symmetric lens, the quadratic term inx i is diagonal and can be written as (here y is taken to be along the where primes denote radial derivatives and we defined a J ≡ (1 − ψ ′′ J )/2 and b J ≡ (1 − ψ ′ J /x J )/2. At this point, in the diffraction integral Eq. (1) we shift and rescale x i to z i ≡ √ wx i . In terms of z i , the quadratic term at the exponent is w independent. The cubic and quartic terms of Eq. (9), multiplied by w, scale instead as O(w −1/2 ) and Core size xc O(w −1 ) respectively. For large w, we can then Taylor expand the exponential and keep terms up to O(w −1 ). At order O(w 0 ) we recover the GO result (from the quadratic part). The term of order O(w −1/2 ) vanishes since it leads to an odd integrand. The first correction comes instead at order O(w −1 ), where we have two distinct contributions: one from squaring the term (ϕ J ) ,ijkx ixjxk and another from (ϕ J ) ,ijklx ixjxkxl . Terms with higher powers ofx contribute at order O(w −3/2 ) or higher, and can thus be neglected at sufficiently large frequencies.
After performing these two integrals, one is left with the following simple result Here ψ (n) ≡ d n dx n ψ. Equation (11) shows that the leading-order GO result is a good approximation provided that ∆ J /w ≪ 1 for all images. 1 Note also that non-analytic features in the Fermat potential (e.g. cusps) produce other O(w −1 ) contributions without a corresponding GO image [27]. We will now address the contribution of non-analytic features in specific cases.

Contribution from the cusp
The leading terms in the GO expansion, Eq. (6), arise from the stationary points of the Fermat potential and capture the high-w contributions to the amplification factor. Nonetheless, other locations in the lens plane can induce corrections at subleading order in the ∼ 1/w n expansion and might be comparable to the bGO terms. In particular, they can arise from singular points of the lens equation (cusps in the lensing potential). See [27] for a similar discussion on cusp contributions to F (w).
In this Subsection we are going to discuss these contributions for the lens models featuring a central cusp (gSIS and SIS lenses, Table I). 2 In particular, we focus on the strong lensing regime, where y can be taken as a small number.
Let us consider the gSIS lens in the limit of large w. For this lens, we distinguish two behaviours depending on the value of the slope k. For 0 < k < 1 (broad profiles) the lens equation is smooth at the lens' centre and a central image forms in the strong-lensing regime. In other terms, the deflection angle α is bounded as the ray approaches the centre of the lens. In the complementary interval 1 ≤ k < 2 (narrow profiles) the lens equation is singular at the centre and no image forms. In both these cases we isolate the contribution to F (w) from the centre by truncating the integration range from x ∈ (0, R c ), for some radius R c small enough for the GO images not to be enclosed. The range x > R c , at high w, is then dominated by the GO expansion around the minimum and/or saddle (depending on y). In the lower integration interval, we have instead where ϕ c ≡ y 2 /2 − ϕ m is the time delay associated to the lens centre (here we re-introduced the minimum time delay ϕ m ) and J ν (z) is the Bessel function of the first kind (obtained after performing the angular integral). The integrand, in the limit w ≫ 1, peaks around x = 0 once we rotate the integration line into the complex plane. To see this, first notice that J 0 (wyx)e iwx 2 /2 ≃ 1 for small x (we will motivate better why x can be taken small aposteriori ). By writing x e −iwψ(x) = e log x−iwψ(x) ≡ e Ω , we can locate the peak as the stationary point x s of Ω: Here we defined A ≡ 1/(2 − k), which is a positive quantity. Equation (14) is solved for x s = (iw) −A . Notice that x s becomes smaller for larger w, making our approximation adequate in this limit (in particular the Gaussian part at the exponent can be neglected since x 2 s ≪ ψ(x s )). Therefore, F c (w) for w ≫ 1 can be obtained using a saddle-point approximation around x s . However, we prefer to take a slightly different approach that yields very similar results: we evaluate J 0 (wyx)e iwx 2 /2 in Eq. (13) at the peak x s , while performing the exact integration over for e Ω . Since the integral is highly localized for w ≫ 1, the calculation can be simplified by taking R c → ∞, making exponentially small errors. We obtain This formula is valid when the limits w ≫ 1 and wyx s ≪ 1 are satisfied (therefore, for small enough y). 3 The full F (w) is then given by the sum of Eq. (15) and the usual GO expansion for the other images. In the following, we will refer to this expansion as resummed GO (rGO). We can better understand the behaviour of F c (w) by first looking at the SIS case k = 1. Here, neglecting again the Gaussian and the Bessel function, we have F c (w) ≃ i/we iwϕc . This can be interpreted as an additional bGO contribution from the cusp x = x c = 0, with time delay ϕ c and with a vanishing GO term (i.e. not accompanied by an image). More in general, for narrow (broad) profiles, F c (w) decays faster (slower) than 1/w. For broad profiles, there is a caveat in the previous derivation at very large w: the Gaussian part can start contributing significantly to the integral, thus leading to the usual GO expansion for the central image. Therefore, for 0 < k < 1, F c in Eq. (15) is a better approximation than bGO for the central image only in the range 1 ≲ w ≲ ∆ c , while for ∆ c /w ≪ 1 bGO performs better (here ∆ c is the bGO coefficient of the central image, Eq. (12)). This issue does not arise for k > 1, since here there is no central image.
From the discussion above we conclude that WO effects from the cusp are relevant even when no central image forms. As we will study in [34], this has interesting implications for parameter estimations with GWs.

B. Low-frequency expansion
We are now interested in understanding the behaviour of the amplification factor in the limit of small w. In this limit, GO fails and one has to resort to other methods to obtain good approximations.
For small w, F (w) approaches 1 since the wavelength becomes much larger than the lens' characteristic scale, and the wave is unperturbed by the lens. Here we would like to motivate that corrections to F (w) ∼ 1 in this limit correspond to an expansion in powers of the lensing potential ψ(x). A physical motivation can be given as follows. If the wavelength is much larger than the typical scale of the lens (i.e. Einstein radius), then the impact parameter's value cannot be precisely resolved. This implies that the impact parameter should be irrelevant (at least at leading order) in this low-frequency limit. Thus, we could imagine performing the calculation for F (w) with y ≫ 1 (i.e. impact parameter much larger than the scale of the lens, set by ξ 0 ) but still much smaller than the wavelength ∝ 1/w. In this case one can expand the diffraction integral in powers of the lensing potential (see [43] on the conditions for the applicability of this approximation).
We can also see explicitly that this procedure gives a sensible series expansion in w: higher powers on ψ(x) lead to subleading terms in w. For simplicity, we show this for axis-symmetric lenses.
First, we perform a rotation of the integration contour to make the integrals manifestly convergent. Similarly to Eq. (13), the amplification factor can be written as follows where in the second line we rescaled the radial variable and rotated the integration contour by 45 degrees in the complex plane: x = e iπ/4 z/ √ w. 4 Note that the Gaussian part dictates the leading behaviour at infinity of the integrand, since the Bessel function only grows exponentially and we assume for convergence that |ψ| grows more slowly than the argument of the Gaussian (lim x→∞ |ψ|/x 2 = 0, see Eq. (19) below). Thus, convergence is now manifest. However, this choice for the integration contours is not optimal: the Bessel function and the lensing potential can make the integrand exponentially large at intermediate values of z, for w large. On the other hand, for small w, Eq. (16) is easy to evaluate, giving us a tool to explore F (w) in the WO regime. We will discuss how to obtain the optimal integration contour for more general lenses in Sec. V.
After this first step, we can expand the integrand in powers of ψ(x(z)). Equation (16), expanded in ψ(x(z)) up to quadratic order, has the form Obtaining higher terms in this expansion is also trivial. If ψ is well behaved for large z, then the integrals in Eq. (17) are localized around z = 1 and can be estimated through a steepest-descent calculation. 5 See [37] for a similar expansion, but in the weak-lensing regime. Let us first consider y = 0 for simplicity. With the z = 1 approximation we obtain This is a series in power of wψ e iπ/4 / √ w : higher powers in the expansion are suppressed for small w provided that This condition is physically sensible since it is equivalent to the requirement for ψ(x) to grow less than the quadratic part of the lensing potential at infinity (in this particular direction of the complex plane). Equivalently, in terms of the density profile ρ, the requirement translates to ρ decaying faster that ∝ 1/r at infinity. All the analytic lens models we consider in this work satisfy this requirement, hence this expansion is applicable in our cases.
The conclusion that Eq. (17) is a good expansion for small w is not spoiled in the case of y ̸ = 0. To see this, we first notice that y enters only the combination wy 2 in Eq. (17). Hence, if wy 2 is smaller than 1, we can expand the Bessel function in a series around zero (the integral is localized around z = 1 and the argument of J 0 remains small). In the limit w ≪ 1 this is possible, with the mild requirement of y ≪ 1/ √ w. Additionally, we can also notice that y will not enter at first order in the w expansion of F (w). Indeed, the series expansion shows that effects due to y are suppressed by additional powers of w.

Leading corrections for symmetric lenses
After these general results, we can focus on the low-w behaviour for the specific lenses shown in Table I. For the simplest cases, we can directly integrate Eq. (17) without assuming small y.
Point lens: In this case ψ(x) = log x and the integral in Eq. (1) has a closed-form solution [8]: where is the confluent hypergeometric function. We can nonetheless employ the low-frequency expansion and then compare with the formula above, expanded in the same limit. Equation (17) (neglecting w 2 ψ 2 terms) gives where Ei(z) is the exponential integral. We checked that indeed Eq. (20) expanded for small w, but fixed wy 2 , yields Eq. (21). It is also useful to further expand Eq. (21) for small wy 2 : where γ E is the Euler's constant. Let us briefly digress and comment about the analytic properties of Eq. (22) as a function of w. We can first note that the diffraction integral, Eq. (1), is analytic for Im w > 0, while possible non-analyticities can appear on the real axis and in the lower half of the complex plane. The analyticity property for Im w > 0 is a consequence of causality: while the lensed signal can have a distorted waveform due to diffraction, no signal can travel beyond the light cone (see [44,45] for a related discussion). The above properties can be checked for the point lens, using its closed-form solution Eq. (20). The latter formula is indeed analytic for Im w > 0, has poles due to the Gamma function for Im w < 0 and a branch cut due to the log w for negative w. Notice that by construction F (w) satisfies a reality condition F * (w) = F (−w). Non-analyticity of F (w) and the reality condition imply that F (w) must be a complex number for w real. Otherwise it would be either real or imaginary. By direct inspection, Eq. (22) is consistent with these properties. In particular, a branchcut appears due to log w. Similar considerations apply to the lens models of Tab. I.
SIS lens: Here ψ(x) = x and keeping up to the ψ 2 term in Eq. (17) we obtain (again keeping wy 2 fixed) After expanding the expression above for small y we obtain Interestingly, in this case the leading behaviour in w differs from the one of the point-lens. The dependence on √ w can be understood from Eq. (17). If we call s ≡ √ wy and use ψ(x) = x, we see that the only dependence on w is through √ w. Together with the fact that s will not appear at leading order, we obtain the correct dependence of F . The slope as a function of √ w is related to the steepness of ψ(x) far from the centre of the lens. This can be seen by considering a generalized SIS lens, with a generic slope.
Generalized SIS lens: The gSIS lens, as already discussed, is described by the lensing potential ψ(x) = x 2−k /(2 − k), where the slope is parametrized by 0 < k < 2. The case k = 1 reduces to the SIS.
The first-order expansion in w leads to where Γ(z) is the Gamma function, L α (z) are the Laguerre polynomials and we introduced again A = 1/(2 − k). Expanding for small y leads simply to This expression agrees with what is found in [37]. Notice that the leading term in w depends on the slope k: steep lenses (k close to zero) have a weaker dependence on w. CIS lens: This lens represents a SIS lens with a central core of size x c . Its lensing potential is given in Table  I. As we have argued, in the limit of low w we are mainly sensitive to the profile of the lens far away from the centre. Thus, we expect to have a weak dependence on x c in this limit. We do not know a closed expression for the integrals of ψ(x) in Eq. (17). However, we can approximate them by noticing that since the integral is peaked around z ∼ 1 and we are interested in w ≪ 1, we can simply expand ψ e iπ/4 z/ √ w for large argument. Interestingly, In other words, the lens looks like an SIS plus a central point lens with negative mass. Notice that to correctly capture the dependence on x c we need to include the second-order correction from Eq. (17). By doing so, we obtain More explicitly, further expanding for small y we have This result shows that the leading behaviour in w is due to the SIS part, while the presence of the core is only visible at subleading order. From Eq. (28), we notice that increasing x c leads to a smaller |F (w)|. This is consistent with our numerical results. Later, we will test the accuracy of this approximation against the full WO result (Sec. VI B, Fig. 9). Throughout this discussion we suppressed the minimum time delay ϕ m (see discussion below Eq. (4)). To reintroduce this parameter, all the expanded expressions for F (w) in this Section should be multiplied by the phase e −iwϕm .

gSIS series expansion
The amplification factor for simple lens models can be written in a useful series representation, by expanding the integrand in powers of the lensing potential ψ(x). In particular, we will be able to obtain such representation for the SIS and gSIS lens models (see Tab. I). This provides an additional independent test of our numerical methods, which will be very valuable to validate our results. Let us first consider a generic axially-symmetric ψ(x) and later specialize to particular functional forms. To proceed, we start from Eq. (1) and perform the angular integral, which yields the usual Bessel function J 0 (wyx) (already encountered in Eq. (17)). After this, we expand e −iwψ(x) in powers of ψ(x) and notice that when ψ(x) is a power-law function of x, each integral in the series expansion can be performed analytically. In the case of a gSIS lens we have where once again A = 1/(2 − k). Notice that for low w, this expansion reduces to our approximation obtained in Eq. (25). Moreover, in the particular case of k = 1 (SIS) we re-obtain the series representation first derived in [26]: In order to recover the oscillatory features of F (w) accurately for large-enough w, the series in (29) and (30) need to be truncated at relatively high values of n. In our comparisons this truncation is made after reaching a 10 −15 precision. As an example, for w = 200 and y = 0.3 this is obtained at around n = 740. For lower w at fixed y, convergence is reached at lower n. As a result, the features due to the GO, bGO and rGO terms are not manifest from this series expansion.
The series (29) is impractical for most applications, since the evaluation of the Laguerre polynomials is slow. Moreover, for imaginary argument they grow exponentially in n. Therefore, many terms in the series are required to reach convergence, even at moderate w. Thus, as we will see, these results are outcompeted by numerical implementations of the diffraction integral in practical applications. Additionally, this series expansion is difficult to generalize to other lens models such as the CIS lens, since in that case we do not have a closed-form solution for the integrals.

IV. REGULARIZED CONTOUR FLOW
We will now turn to general methods to numerically solve Eq. (1). This Section presents the regularized contour flow, a calculation performed by Fourier transforming the integral in Eq. (1) and evaluating the resulting time-domain integral on contours of equal time delay. Each set of contours is then flowed adaptively to a different value of the time delay until "hitting" a critical point, when the contour ends (several contours end at saddle points). The total integral is transformed back to frequency space by means of a fast-Fourier transform (FFT), after splitting the result into regular, smooth and singular contributions, associated to GO results. We follow Ulmer & Goodman [30] (see also [46]) but use a different regularization used for the saddle points. The steps in the method are described in Fig. 2, the precision parameters and their default values in Table II. A. Adaptive Sampling in the Time Domain We will first compute the amplification factor in time-delay space. Fourier-transforming the integrand in Eq. (1) yields where δ D (x) is the Dirac-delta function. This expression is the primitive of the Green's function: the timedomain lensed waveform is given by h(t) = dτ dĨ(τ − t)/dτ h 0 (τ ). One can compute the functionĨ(τ ) by binning in time delays, cf. Ref. [14,15]. Instead, in this approach we are going to reduce the 2D integral into a 1D integral over the contours where the argument of the Dirac delta is zero: where s is the arc-length distance that parametrizes the k-th contour γ k of constant τ (it is possible to reparameterize the contour so the integrand in Eq. (32) reduces to du, cf. Ref. [30,Eq. 9]). The sum is over all contours contributing to a given τ (see below).
Contours of constant Fermat potential are orthogonal to the gradient of ϕ, hence given by To sampleĨ(τ ), we flow the contours as i.e. each point is displaced (linearly) in the direction of the Fermat potential's gradient, with ∆τ positive/negative if τ is increasing/decreasing. The integral I(τ ) is sampled adaptively on each family of contours. We choose the step ∆τ depending on the rate of varia-tionĨ, according to the following prescription Here η, η ′ are precision parameters and dĨ/dτ, d 2Ĩ /dτ 2 are computed numerically from the previous iterations. The variation ∆τ is kept within a minimum and maximum values. The algorithm stops if a prescribed number of iterations is reached. At each step τ i we refine the contour so the distance between its nodes is small compared to the local curvature radius of the contour (see Table II). Contours are flown until reaching a critical point, where they either shrink to a point (maxima, minima and cusps) or become nondifferentiable (saddle points). Therefore, we first find the critical points x J such that ∇ x ϕ(x J , y) = 0 (i.e. GO images), Eq. (5). Since we work with symmetric lenses we search only along the lens-source direction x 1 ≡ x · y/y ( Fig. 2 panels A, B). For a given value of τ there can be zero, one or multiple contours depending on the lens configuration. The total number of contours changes when τ crosses the values of the time delay ϕ J associated to critical points and the integrand in Eq. (34) becomes singular. Hence, critical points are associated with singularities and discontinuities ofĨ(τ ). Contours do also end in non-regular points of the lensing potential, which are not associated to GO images. 6 In the weak-lensing regime there is a single type I image at the minimum of ϕ: the contours then flow be-tween the minimum and infinity, where they approach circles/ellipses centered around the source. Strong lenses produce multiple contours in certain ranges of τ . For the SIS, gSIS and CIS in the strong-lensing regime there are three regions to consider: 1) from infinity down to the saddle point, 2) from the minimum up to the saddle point and 3) from the maximum down to the saddle point. Contours that begin near minimum/maximum are initiated with a small radius around the critical points. The contour that asymptotes to infinity is initiated at a large radius around the source (see Table II). These regions and their contribution toĨ(τ ) are shown in Fig. 2 (panels B and C). Here d l,m is the distance between the l,m-th node points and R l is the local curvature radius of the contour evaluated at the l-th node.

B. Time-domain Regularization and GO Counterterms
The frequency-domain amplification factor can be computed from Eq. (31) via inverse Fourier transform. We will perform this operation via fast-Fourier transform, for which we interpolateĨ(τ ) on an equally spaced grid with N FFT points and spacing δτ , whose range and density are determined by the minimum and maximum values of w that we are interested in. In order to avoid boundary effects, we apply a Tukey window function to the time-domain signal. We extend τ towards negative values, such that the window function is ≃ 1 over all ϕ I corresponding to GO images.
One difficulty is dealing with the discontinuities and singularities inĨ(τ ) associated to GO images, as they produce components at arbitrarily high frequency. The discretization needed for the FFT causes aliasing of the frequencies w > w max , contaminating the computed signal at high w (Fig. 2, panel D).
To avoid numerical artefacts it is convenient to treat the contribution of stationary points separately [30]. Hence we split the integral into a regular and a singu- The l.h.s. is obtained by evaluating Eq. (32) numerically, as described above. The singular contributions have closed-form expressions, which we give below. The WO amplification factor follows by Fourier transforming The regular contribution is the FFT fromĨ reg (τ ). The singular contributions are the GO amplification factors for each image, Eq. (6), or related to them. Note that the above splitting is arbitrary and valid as long as the time and frequency domain terms are consistent. Therefore, we can add any such terms in order to make the computation more robust. We will now discuss these terms for different types of critical points. Type I/III images (minima/maxima of the Fermat potential) correspond to discontinuities inĨ(τ ). In this casẽ where θ is the Heaviside step function and +/− corresponds to a minimum/maximum with time delay ϕ J . The corresponding GO contribution reads The discontinuity is interpreted as a family of contours ceasing to exist at the extremum. In the case of a cusp (x → 0 in the SIS and gSIS) the contour ceases to exist but no discontinuity forms because µ J = 0. A type II image (saddle point) with time delay ϕ J produces a logarithmic divergence inĨ(τ ) [30] with C an integration constant. This calculation assumes that the quadratic approximation around the saddlepoint, ϕ ≃ (∆x 1 /a) 2 − (∆x 2 /b) 2 , is valid for arbitrarily large separations ∆x. 7 The corresponding GO contribution reads 7 Including only a finite region around the saddle point yields where δτ defines the limit of the contour around the saddle point. Equation (40)  This term contributes at arbitrarily large values of τ , causing spurious low-frequency behaviour upon FFT. One can avoid this issue by windowing the singular con-tributionĨ S sing (τ ) → W (τ, T )Ĩ S sing (τ ) in Eq. (40). Choosing W (τ, T ) = e −|τ −ϕ J |/T (43) preserves the singular behaviour and avoids the low frequency problems if δτ ≪ T ≪ δτ N FFT , where N FFT is number of sampled points. This choice produces a closed-form expression for the frequency domain where Note that lim T →∞ I ± = ∓ i w γ E + log(w) ∓ i π 2 , recovering Eq. (42).
Panel C of Fig. 2 shows the integral in the time domain, including both the regular and the different singular contributions. Panel D shows the different contributions to F (w): at low frequencies the singular part (corresponding to GO predictions) is compensated by the regular contribution and recovers the low-frequency limit F (w → 0) → 1. Without splitting the singular part, F (w) loses precision at high frequencies and eventually becomes unreliable (solid gray line).
After the regularization, discontinuities remain in derivatives ofĨ(τ ) at the critical points τ = ϕ J . When Fourier transformed, discontinuous d dτĨ (τ ) corresponds to corrections ∝ 1/w in the amplification factor. These terms are precisely the bGO and cusp contributions discussed in Sec. III A. Discontinuities extend to any derivative d n dτ nĨ (τ ), with corresponding corrections ∆F (n) ∝ w −n as higher order bGO terms. Numerically, the discontinuous derivatives cause aliasing in the FFT and are a source of error, although the scaling with w makes these terms subdominant in the computation of F (w) (see discussion in Sec. VI A). Eventually, our regularization method could be extended to split discontinuities in derivatives ofĨ(τ ) as higher precision is required.
Finally, let us mention that computing derivatives of the amplification factor accurately requires handling discontinuities in d dτĨ (τ ). This is because an additional derivative (e.g. with respect to the lens parameter Θ) promotes the w n terms in F to ∼ w n+1 in dF/dΘ: our regularization removes ∼ w terms, but leaves terms ∼ w 0 , which contribute significantly to aliasing. As computing derivatives of F (w) is important in some applications of WO lensing (e.g. Fisher-matrix forecasts [34]), we discuss a workaround in Appendix A.

V. COMPLEX DEFORMATION
Another method to improve the convergence of the diffraction integral (1) relies on analytically continuing the integration variable into the complex plane. We follow a variation of the procedure described in [32] (see also [47][48][49]). Initially, we briefly review how to obtain a good integration contour in the complex plane for highly-oscillatory 1D integrals. Then, we generalize the discussion to the more relevant case of 2D diffraction integrals with non-analytic features, as appear in gravitational lensing.

A. Flow of the integration domain
To outline the method, we start with a prototypical oscillatory integral, which resembles (1) in one dimension where f (x) and g(x) are assumed to be analytic in the full complex plane z = x + iy ∈ C (later we will mention how to deal with possible branch-cuts or poles away from the real line). 8 Clearly, the limit of large w corresponds to the saddle-point, or GO, approximation for Eq. (46) (when dealing with complex functions, critical points can only be saddles). First, let us recall how to make sense of this type of oscillatory integrals, which are not manifestly convergent (see e.g. [50] for more details on this procedure). The integral (46) can be extended to the complex plane z = x + iy by deforming the contour of integration away from the real line (because of Cauchy's theorem the final answer does not change). In particular, the contour can be decomposed into so-called steepest-descent paths, each associated with a saddle point. (Some saddles might be irrelevant in this decomposition and neglected, as they are not encountered when moving the real-line integration contour into the complex plane). The steepestdescent path z(λ) of the function f (x) associated with an isolated and non-degenerate saddle z c (i.e. f ′′ (z c ) ̸ = 0) satisfies the flow equation where λ ∈ R parametrizes the path and z * is the complex 8 The distinction between f (x) and g(x) is ambiguous in Eq. (46). To partially fix the ambiguity, we assume that the convergence properties of the integral at infinity are solely determined by f (x). As we are going to elaborate, convergence is established by deforming the integration contour into the complex plane (in particular, by tilting the contour above or below the real line, depending on the properties of f (x)). conjugate of z. From the equation above, it follows that Therefore, the real part of the exponent of the integrand (46) monotonically decreases along the steepestdescent path: the integral is now manifestly convergent, as the integrand decays exponentially. 9 On the other cedure, introduced in [32] and discussed below, is what we will use in our applications. We start by considering again Eq. (47). In this flow equation now we impose as initial condition at λ = 0 that z(0) = x, where x is a generic point on the initial integration contour C λ=0 (in the case of Eq. (46), C λ=0 = R). Let us call C λ the set of points z(λ) at the "time" λ. As a result of Morse theory, the contour C λ for λ → ∞ converges to a steepest-descent contour (Eq. (47) represents a smooth deformation of the initial path C λ=0 ). Therefore, just by repeat use of Eq. (47) for all the points of the initial domain of integration, we can write the original integral for large enough λ as where the exponent is now real (up to a z-independent imaginary part) and not oscillatory. Standard numerical techniques can now be applied to this integral. Notice that C λ does not depend on w , therefore the path needs to be computed only once. This property applies to GW lensing, where w appears only linearly at the exponent. In situation where the Fermat potential depends on the frequency, as for wave diffraction in dispersive media, this is no longer the case.

B. Extension to realistic lenses
The procedure just outlined can be generalized to higher dimensions and applies to the diffraction integral (1). For our practical purposes, however, it is simpler to reduce F (w) to a set of one-dimensional integrals and apply the procedure above. To achieve this, we first write the integral over x in polar coordinates x and θ as The integral in θ is over a finite range and can be performed with standard numerical techniques. On the other hand, the integral over x is highly oscillatory and is suitable for the analytic continuation procedure. The contour-deformation method has to be applied for various values of θ, until the sampling of points is dense enough to guarantee numerical convergence. The only difference with Eq. (46) is the lower limit of integration, which stops as x = 0. In order to apply Cauchy's theorem, the initial and final contours C λ=0 and C λ must close (up to a semicircle at large |x|, which is negligible). However, if the point z(0) = 0 is evolved according to the flow, it will in general move away from the origin, thus leaving the sum of the initial and final contours open. To avoid this issue, we decide to alter the flow equation (47) in such a way as to force the points close to the origin not to evolve. To do so we rescale the flow variable λ → h z0 λ, where h z0 is a function that depends on the initial position z 0 ≡ z(0). Note that changing the flow equation does not modify the final result due to Cauchy's theorem. The flow equation then becomes For convenience we choose h z0 to be h z0 = θ ε (z 0 − δ), where θ ε (x) ≡ 1 2 [tanh(x/ε) + 1] is smooth, and converges to the Heaviside step function for ε → 0. The function h z0 interpolates between h z0 ≃ 0 for points z 0 < δ and h z0 ≃ 1 for z 0 > δ. The parameter ε sets how sharp the transition is. In the application of the following sections we set ε = 10 −3 and δ = 10 −2 .
An additional complication arises when considering realistic lenses, which typically feature non-analytic lensing potentials (even the simplest example, the point lens ψ = log x, has a branch-cut ). Fortunately, our procedure is not significantly altered as long as we are dealing with branch cuts (and possibly poles, but we are not going to encounter them in our lensing models). Away from the branch cut the exponent is analytic, and the flow equation can be applied without modifications. It is possible however that some points are driven towards the branch cut during the flow. In such cases, to avoid them crossing the cut, we decide to stop the flow. This can be implemented case by case (depending on the location of the cuts) by modifying the flow. Thus, we multiply the righthand side of Eq. (51) by a function b(z(λ)) with the requirement that for z(λ) approaching the cut, b(z(λ)) → 0 smoothly, while b(z(λ)) ≃ 1 everywhere else. For the case of a point lens we have a branch cut for Re z < 0. Following the same logic as for h z0 , with this lens we choose b(z) = 1 − θ ε (φ − π δ ) − θ ε (−φ − π δ ), where z is written in polar coordinates z = re iφ and we defined π δ ≡ π − δ. With this choice, the evolution of the contours is halted as φ approaches ±π. For the function b(z), in the next section we will use the values ε = 1/200 and δ = 10 −1 .
Since we are modifying the flow equation, some of the nice properties of Eq. (47) are partially lost. In particular, it will no longer be true that the imaginary part remains constant along the final contour C λ (hence some mild oscillations can reappear). In practice, for the cases we will consider, this is not an issue since the modifications only affect the contour close to the origin, while leaving the behaviour at large |z| unaffected (the convergence properties of the integral are thus preserved).
We also notice that the integration contour at infinity converges to the 45-degrees line arg x = π/4, Re x > 0 (for large x the Fermat potential is dominated by its quadratic part). The contour is however deformed as x approaches the origin, due the Lensing potential. Therefore, the method outlined in this Section generalizes Eq. (16) and optimizes the choice for the contour.
Let us see how this procedure is applied in a particular lens model, the SIS lens. First, we fix a value for the angular variable θ in Eq. (50). We take it to be θ = π/4 in this example (Fig. 3, panel A). Then, we evolve the integration contour from the positive real line x > 0 to the complex plane. Since this lens model does not introduce branch cuts or poles, the evolution of the path is obtained using Eq. (51). Then, it is stopped at some given flow time λ = T and the path is truncated at some large value of |x|. The evolution of the paths is shown in Fig. 3, Panel B. We can notice that all the paths cross the real line at a fixed location: this corresponds to a saddle point of the exponent in Eq. (50) (this is however different from the saddles corresponding to the GO approximation, since we are fixing θ here). Once the paths are obtained, we can evaluate the integrand in Eq. (50) (that we call F) on each path. This quantity is now also a function of w. As an example, we show F for fixed w and θ in Fig. 3, Panel C. Clearly, the problematic oscillations in F get damped very quickly as the flow progresses.
For large values of λ the integral becomes very localized around the saddle point. Finally, one has to repeat the steps A-C for different values of θ. The integral of F over x as a function of θ is shown in Fig. 3, Panel D: for moderate w, this function is not oscillating too rapidly and can be integrated easily.
In the results presented in the following sections, we use the following settings for the complex-deformation method. The contours are evaluated for n 1 = 25 values of θ, uniformly distributed between 0 and 2π (for symmetric lenses one can limit to the range 0 to π). The integrand is then sampled over a larger number of values of θ, n 1 · n 2 , with n 2 = 25. We can sidestep evaluating n 1 · n 2 contour flows because the final integral is independent on the choice of the contour. Therefore, for a given angle, we use the contour evaluated with the nearest value of θ.
For the flow of each contour, we sample the initial contour along x > 0 with n x = 340 points, not uniformly distributed but concentrated towards x = 0. The flow equation is implemented in python, using the odeint function of the scipy package on default settings. After applying the flow equation, the final contour is interpolated over n interp x = 1000 points.

VI. ACCURACY AND PERFORMANCE
We will now discuss the accuracy of the algorithms described in Sec. IV and V and their convergence to systematic expansions, Sec. III. We first compare the results for a point lens (Sec. VI A) and then against each other (Sec. VI B) for the other axially-symmetric lenses in Table I. We end by discussing the performance of the different methods (Sec. VI C). All the comparisons are made with impact parameter y = 0.3; similar conclusions are reached for different values of y, sufficiently far from caustics. The precision parameter used in the contour method are summarized in Table II.

A. Comparison to point lens
In order to assess the goodness of our numerical methods, we can compare to the point-lens model (ψ(x) = log x), where the diffraction integral is known analytically (see Eq. (20)).
We compare the result from the contour method of Sec. IV and the complex-deformation method of Sec. V with Eq. (20) for y = 0.3 in Fig. 4, where also the results of GO and bGO are reported. The contour method is most accurate in the intermediate-w regime, up to w ∼ 10, while the complex-deformation method remains good even at w ∼ 10 2 .
Let us now discuss the comparison with the contour method more in details, since it will be the primary method used in future applications, [34]. We will discuss the complex-deformation method's performance in the following subsection. At low frequencies (w ≲ 0.1) the numerical result from the contour method stops being accurate. This is related to the way the numerical calculation is performed. Indeed, the signal is obtained through a Fourier transform from the time domain signal I(τ ). Therefore, the low-frequency errors are related to the numerical truncation of the integral in Eq. (32) for large time delays, which depends on the windowing of I(τ ). Higher precision can be achieved at low frequencies by extending the integral (32) to larger τ s, at the expense of making the numerical evaluation slower.
In the opposite regime, for large w, we also lose accuracy. The appearance of error in this regime can be understood in the following way. At high w, the signal can be written as follows where the sums are over the different images J and the higher-order GO corrections that scale as ∼ w −n . Here for n = 0 we recover GO (c (J) 0 = 1) and for n = 1 we have instead bGO (c (J) 1 = i∆ J ). The term F (nmax) reg (w) represents the regular WO contribution, not captured by the GO up to order n max . We can notice that all the GO terms, when Fourier-transformed to the time domain τ give some "singular" features. In particular for n = 0, as we already discussed in Sec. IV B, we can have θ-function discontinuities or log divergences inĨ(τ ). As discussed there, there are also discontinuities/singularities on the n-th derivative ofĨ(τ ) with respect to τ , for any n. Due to finite numerical accuracy, such sharp features pollute the signal at arbitrary high w when transformed back to frequency space. Having understood this, we have a strategy for potential future improvements in the accuracy of our code. Indeed, these additional GO contributions could be subtracted before performing the inverse Fourier transform, in the same spirit of what is already done in the case of n = 0. By performing this procedure up to n = n max , we expect the residuals against the full result to scale as ∼ 1/w 1+nmax . Of course, other sources Methods comparison for a point lens (y = 0.3). Top: Absolute value of the amplification factor |F (w)|. The contour and complex method are not shown explicitly, as they overlap with the exact solution. Bottom: Differences relative to the exact solution. The contour method (blue) always performs at the sub-percent level, while the complexdeformation method (orange) at low frequencies is around one order of magnitude better. The GO approximation (light green) converges very quickly, remaining sub-percent for w > 10, whereas the bGO approximation (red) converges even faster.
of error might then become dominant. We expect that these remarks also apply to other lens models.

B. Comparison between methods
Here we compare the numerical results from the contour method and the complex-deformation method introduced respectively in Sec. IV and V. Again, we will consider the extended lenses described in Table I. For the SIS lens, as we have discussed in III B 2, a series representation for the integral is available and can be used for comparison with our numerical methods. A comparison with the series Eq. (30) is shown in Fig. 5, again for y = 0.3. We can see that the agreement for both our numerical methods is below the permille level in the range 10 −1 ≲ w ≲ 10.
At low frequencies the contour method starts to fail, as seen from the non-smooth curve for F . On the other hand, the complex-deformation method works best in this regime. It is more accurate since the (typically oscillatory) angular integral over θ in Eq. (50) is not par- ticularly computationally demanding. 10 Going instead to higher w reverses the situation: the complex-deformation method becomes slower and is outperformed by the contour one. The main reasons are the increasing oscillations in the angular integral and the fact that the computation has to be performed frequency-by-frequency (the contour method instead evaluates F directly at all frequencies).
In our code, the accuracy of the complex-deformation method can be improved by increasing the sampling of the path C λ and of the integral over θ. This is at the cost of a longer evaluation time.
Let us compare different methods for computing the gSIS amplification factor. We test the analytic result of Eq. (29) against the contour and complex-deformation methods. Additionally, we also check the improvement of the rGO approximation of Eq. (15) against the GOonly result. The results are shown in Fig. 6. As we can  see, the general results found for the SIS lens in Fig. 5 are preserved here, with similar orders of magnitude for the accuracy against the analytical result. Also, we can notice the slightly better agreement between the rGO compared with only GO (the bGO curve, not shown, gives a similar residual as the GO one). With these results, we can establish that the contour and complex-deformation methods are accurate enough to study the gSIS lens. The final lens we consider is the CIS. In this case, as already discussed in Sec. III B 2, a series representation is not available. Therefore, in Fig. 7 we compare the complex-deformation method against the con-  Comparison between the low-w expansion F CIS of Eq. (28) and the complex-deformation method of Sec. V for a CIS lens with y = 0.3. Different curves correspond to different values of the core size xc. The analytic expansion reaches per-mille accuracy for w ≲ 10 −2 , while the expansion breaks down around w ∼ 1, as expected.
tour method directly. Moreover, comparisons with the GO and bGO results are also shown. As we can see, the two numerical methods have the lowest residuals around w ∼ 1, whilst at lower and higher frequencies the results degrade. As for the other lenses, at low w the contour method loses accuracy. At higher frequencies instead, the complex-deformation method become less reliable, making the residual bigger. We can notice that the trends are very similar to those of Fig. 5, 6, for the SIS and gSIS lenses. In comparing against GO and bGO we find similar features as for other lenses: as expected the residuals scale as ∼ 1/w and 1/w 2 respectively. From this observation, we can argue that the contour method remains accurate below the percent level for 1 ≲ w ≲ 300, at least. If that was not the case, the residual w.r.t. GO and bGO would saturate at high frequencies. For CIS there is no rGO curve, as the centre of the lens is regular and there is no cusp contribution.
We tested our results also at different impact parameters, obtaining similar agreements between methods. However, accuracy degrades as we approach a caustic where two images merge. Since images are closer in this situation, a lower precision is expected: resolving them requires higher resolution in the time delays (for the contour method) and higher resolution in the image plane (complex-deformation method). In the case of the CIS, a caustic occurs at an impact parameter y rc (see eq. (64) in Ref. [34] for the explicit expression). A comparison at a larger impact parameter is shown in Fig. 8.
For SIS and gSIS lenses, the series representations (Eq. (30) and (29) respectively) reduce to the low-w expansions (Eq. (24) and (26)) for small enough w. Therefore, the comparisons against the numerical methods in Fig. 5, 6 implicitly show the goodness of this approximation. To perform a similar test for the CIS lens, we compare the low-w expansion Eq. (28) against the complexdeformation method in Fig. 9. The low-frequency approximation performs at sub-percent level for w ≲ 10 −1 , while per-mille accuracy is reached for w ≲ 10 −2 . Notice that the contour method would lose accuracy in this region. These results have a mild dependence on the value of x c , with more accurate results for larger x c .

C. Performance
The contour method is by far the fastest computationally. Our implementation in python (optimized with Numba, after pre-compiling) computes F (w) in the range w min ∼ 0.01, w max ∼ 1000 on a 12-cores laptop (i7-10750H CPU) in 323 ms, a range similar to the one used for Figs. 4-7 (the numbers refer to the point lens, the extended lenses require similar execution times). 11 The contour method is faster than computing the exact solution for the point lens, (20): it takes 5.1 s on the same machine sampling the same range of w over 10 4 points (a factor 1/10 fewer values of w than quoted above). A single evaluation of the exact solution (20) takes 0.47 ms at w ∼ 0.01 and 5.7 ms at w ∼ 1000.
GO and bGO calculations are very inexpensive, requiring 12 and 15 ms for the same lens, respectively. Because the contour method employs regularization, the GO amplification factor is already computed.
Complex deformation is efficient at low w but each dimensionless frequency needs a separate computation. Moreover, the angular integral becomes highly oscillatory for large w and needs to be sampled very finely to obtain the desired accuracy, at the cost of increasing the evaluation time considerably. With this method, one needs first to evolve the contours of integration (Panel B of Fig. 3) through the flow equation for different values of the angle θ. This evaluation takes around ∼ 300 ms for a single angle and 100 points in the x integration path. Parallelization can be used to speed up the evaluation over multiple angles (for ∼ 50 values of θ the overall evaluation takes ∼ few s). After this step, the evaluation of F (w) at a single value of w takes around ∼ 300 ms. Most of the computational time is spent evaluating the integrand function over the 2D domain of integration. Evaluation over multiple ws is then parallelized: this generally provides a factor ∼ 5 gain in speed. The situation worsens if w is increased too much since the number of sampled angles θ needed for a precise evaluation has to increase (the timing will scale almost linearly with the number of angles).
Given the different execution times between both methods, employing the contour method in applications such as parameter-estimation for GWs is more convenient. The complex-contour method remains a valuable tool for the validation of our results, in particular at low frequencies. For instance, it allows us to test the low-w behaviour derived in Sec. III B for different lens models.

VII. CONCLUSIONS
In this work we have developed, implemented and validated methods to compute gravitational lensing predictions in the WO regime. We first outlined the general framework for WO computations and introduced several lenses that serve as examples for comparison (Sec. II). We described systematic expansions valid at high and low frequencies (Sec. III). The geometric optics (GO) expansion is valid for w ∝ GM Lz f → ∞ (Sec. III A). It gives the amplification as a sum over images, each of which carries a magnification, time delay and Morse phase. Finite frequency corrections ∝ 1/w can be systematically included, with terms that depend on derivatives of the lensing potential at the GO images. Similar 1/w corrections stem from non-smooth features of the lensing potential, which are not associated to GO images and vanish at sufficiently high frequency.
While GO is a local expansion around critical points, a low-frequency expansion highlights the dependence on the lens's global properties and asymptotic behaviour. The w ≪ 1 limit (Sec. III B) allows us to explore the deep WO regime and leads to simple expressions for our example lenses. The leading-order corrections depend on the asymptotic of the lensing potential, a trend clearly seen in the convergence towards free propagation F ≃ 1 as w → 0: • The point lens shows the fastest convergence, with F − 1 ∝ w + O(w 2 ). This follows from ψ ∝ log(x) being the slowest possible asymptotic growth of |ψ|.
• Extended lenses approach the unlensed case more slowly, as F − 1 ∝ w k/2 + O(w k/2+1 ) for ψ ∝ x k , and k < 2 needed to keep the enclosed mass finite.
This difference can be explained by the larger projected mass within a region of radius ∝ 1/w, which dominates the diffraction integral at low frequencies. Hence, we expect that an extended but isolated lens (e.g. truncated at finite radius) will recover the point lens convergence at sufficiently low w. The point-lens limit can not be recovered by the gSIS, as it requires k → 2, leading to divergent mass. While the low-frequency series can be computed, convergence requires including many terms, even at moderately high frequencies. Intermediate frequencies depend on the properties of the lensing potential across the lens plane, which are more efficiently computed numerically.
We developed and presented two numerical methods to compute WO lensing. The regularized contour-flow method (Sec. IV) solves the Fourier-transform of the diffraction integral (1) by adaptively sampling equal-time contours of the Fermat potential. Subtracting the singular contributions and then adding their appropriate terms after transforming back to the frequency domain significantly reduces numerical noise at high frequencies. The complex-deformation method (Sec. V) analytically continues the integration variable in the complex plane. A well-defined process to flow the integration contour allows us to solve convergent, non-oscillatory integrals.
Both methods are complementary to each other. The method of contour flow allows a very fast computation of the amplification factor over a whole range of frequencies using FFT. The computation is efficient enough for GW parameter estimation with LVK data, although further optimization might be necessary to study complex lenses described by many parameters. A main shortcoming of contour flow is that it requires knowing the initial and final conditions for each contour. While the endpoints are solutions to the lens equation, setting up a calculation for a complex setup with many images can become very involved. In contrast, complex deformation does not require knowledge of solutions to the lens equation. However, the method is costly, as each value of the w needs to be computed independently. Moreover, high frequencies require a very fine sampling of the angular integral (this is not the case for contour flow, where one only needs to increase N FFT to reach higher w).
We performed the first cross-validation of these numerical methods for different lensing potentials. We can achieve sub-percent accuracy over a broad range of frequencies with both methods. These results are optimal at intermediate frequencies, and can be matched to analytic expansions whenever the method fails or becomes too costly. Besides comparing with the exact point-lens solution, we found excellent agreement with the series solution for the SIS and gSIS. Having two independent methods allows us to validate the predictions for lenses for which no systematic solution is known, as we demonstrated for the CIS. This will be important when consid-ering more realistic and involved lens models.
Our main conclusions can be summarized as follows 1. High-and low-frequency limits capture different properties of the lens. The large w limit depends on the critical points and non-analytic features, while the low w limit depends on the asymptotics of the lensing potential. Intermediate w requires knowledge of the entire lensing potential.
2. Predictions in any regime can be obtained efficiently by combining numerical methods at intermediate frequencies with analytical approximations at low/high frequency.
3. We find sub-percent agreement to exact solutions, as well as between different methods in their regime of validity. Sub-percent accuracy holds for all the lenses we considered without fine-tuning the precision parameters.
4. Numerical methods offer complementary advantages: contour flow gives frequency-dependent predictions very fast (≲ 1 s on a laptop). Complex deformation does not require prior knowledge of solutions to the lens equation.

These methods offer insights into WO phenomena.
Complex deformation is particularly useful for understanding the low-frequency limit. The contour method gives a transparent interpretation of GO and bGO (as non-analytic features in the timedomain integral or its derivatives) and its splitting from other WO effects.
Our methods provide a baseline for computing gravitational lensing predictions in the WO regime. However, many improvements are possible and will be desirable in the near future. A clear direction will be the extension to more complex gravitational lenses. This requires allowing the methods to work with tabulated values for the lensing potential (rather than closed-form expressions), streamlining steps such as finding the limits of contours or the sampling of integrals in terms of precision parameters. A goal of this program is to make the methods available and integrate them into a public software tool [51,52]. We also envision further developments, such as including a robust computation of amplification factor derivatives to improve lens parameter inference. Our tools can also be generalized to other applications, such as plasma lensing [31,53], studying GW polarization effects [54][55][56] and testing gravitational theories [57][58][59][60].
The sub-percent accuracy we have demonstrated is sufficient for the near-term future of GW observations. In particular, it allows us to model lensed waveforms by LVK and even 3rd generation GW observatories, where we expect signal-to-noise (SNR) ratio ≲ 100. Higher accuracy might be possible by better adjusting the settings of our calculations, although refinements, such as including 1/w terms in the regularization scheme, will eventu-ally become necessary. These improvements will eventually be required for analysing higher SNR sources, such as massive black hole binaries (SNR ∼ 10 3 − 10 4 ) that will be observed by LISA and other space-borne detectors. Now and in the future, these methods will facilitate many applications, from searching WO effects in GW data to novel probes of the matter distribution in the universe and fundamental physics.