Numerical evaluation of oscillatory integrals via automated steepest descent contour deformation

Steepest descent methods combining complex contour deformation with numerical quadrature provide an efficient and accurate approach for the evaluation of highly oscillatory integrals. However, unless the phase function governing the oscillation is particularly simple, their application requires a significant amount of a priori analysis and expert user input, to determine the appropriate contour deformation, and to deal with the non-uniformity in the accuracy of standard quadrature techniques associated with the coalescence of stationary points (saddle points) with each other, or with the endpoints of the original integration contour. In this paper we present a novel algorithm for the numerical evaluation of oscillatory integrals with general polynomial phase functions, which automates the contour deformation process and avoids the difficulties typically encountered with coalescing stationary points and endpoints. The inputs to the algorithm are simply the phase and amplitude functions, the endpoints and orientation of the original integration contour, and a small number of numerical parameters. By a series of numerical experiments we demonstrate that the algorithm is accurate and efficient over a large range of frequencies, even for examples with a large number of coalescing stationary points and with endpoints at infinity. As a particular application, we use our algorithm to evaluate cuspoid canonical integrals from scattering theory. A Matlab implementation of the algorithm is made available and is called PathFinder.


Introduction
In this paper we consider numerical evaluation of the integral where Γ is a contour in C, possibly starting and/or ending at infinity, f and g are functions of a complex variable, and ω > 0 is a frequency parameter.Such integrals arise in numerous application areas, particularly in wave phenomena and quantum mechanics, and are generally challenging to evaluate numerically, especially when ω is large, because the presence of the exponential factor e iωg (z)  means that the integrand may undergo rapid oscillations and/or significant variations in amplitude along the integration contour.
1 arXiv:2307.07261v1[math.NA] 14 Jul 2023 When f and g are analytic, Cauchy's theorem provides the possibility of deforming the integration contour so as to make numerical evaluation easier.This is the basis of steepest descent (SD) methods, in which one aims to deform Γ onto a contour, or, more typically, a union of contours, which we term the steepest descent (SD) deformation, on which [g(z)] is constant, so that the exponential factor e iωg(z) is no longer oscillatory.By the Cauchy-Riemann equations, these contours coincide with the steepest descent curves of − [g(z)], and they connect endpoints of the original integration contour, valleys at infinity (sectors in which the integrand decays rapidly as |z| → ∞), and stationary points of g, which are points ξ ∈ C at which g (ξ) = 0.1 Along each SD contour, away from stationary points the integrand typically decays exponentially, with the rate of decay increasing with increasing ω, and as ω → ∞ the value of the integral is dominated by local contributions close to the endpoints of Γ and any stationary points traversed by the SD deformation.In the asymptotic steepest descent method (described e.g. in [3,19]), one exploits this to obtain an asymptotic expansion for the integral, valid as ω → ∞, by performing a local Taylor expansion of the integrand around the endpoints and relevant stationary points, and reducing the local integrals along the SD contours to simpler integrals that can be expressed in terms of known special functions.
In the numerical steepest descent (NSD) method (described e.g. in [7, §5]) one evaluates the integrals along the SD contours numerically.This involves numerically tracing an appropriate segment of each SD contour in the SD deformation and applying suitable numerical quadrature rules to evaluate the associated contributions to the integral.
In principle, NSD is a highly accurate and efficient method for evaluating integrals of the form (1) for moderate or large ω.Indeed, under appropriate assumptions, the NSD method outputs approximations which, for a fixed number of quadrature points N , are asymptotic to (1) as ω → ∞, with the asymptotic accuracy improving with increasing N (see, e.g., [7,Thm 5.7]).Furthermore, if f and g are sufficiently well behaved it can also be the case that the NSD approximation converges to (1) as N → ∞, for fixed ω > 0, with a cost that remains bounded as ω → ∞.
In practice, however, applying the NSD method to an integral of the form (1) often requires significant expert user input.This is because: (P1) Determining the SD contour deformation corresponding to a given g and Γ requires careful a priori analysis.
(P2) Parametrizing SD contours from or near stationary points, and evaluating integrals along them, is fraught with numerical difficulties, especially when stationary points are close to other stationary points or endpoints of Γ.
The issues described in (P1) and (P2) are particularly troublesome when one wishes to evaluate (1) for multiple instances of a phase function g(z) = g(z, c) depending on a set of parameters c ∈ C q .This is because, firstly, the number and location of the stationary points, and the nature of the SD deformation, have to be determined for each different value of c, and, secondly, stationary points may coalesce as c approaches certain regions in parameter space, leading to a non-uniformity in the accuracy of the resulting NSD approximations.The problem of stationary point coalescence in the context of NSD was studied in detail in [10] in the special case of the cubic phase function g(z, c) = z 3   3 − cz, for c ∈ C, which for c = 0 has a pair of order one stationary points which coalesce as c → 0 (at z = 0) into a single stationary point of order two for c = 0.2 In this case, the SD deformation and contour parametrization were carried out manually by analytically inverting the phase (illustrating (P1)), but the resulting integrals were found to be nearly singular for small c, leading to poor accuracy of standard NSD approximations (illustrating (P2)).It was shown in [10] how to construct a family of non-standard quadrature rules for this integral which perform uniformly well for c ≈ 0 using complex-valued Gaussian quadrature, producing quadrature nodes that in general lie off the SD deformation.In principle, similar rules could be developed for more complicated coalescences involving higher order stationary points and/or endpoints of Γ.However, for each type of coalescence a bespoke quadrature rule would have to be developed, and a general catalogue of such rules is not yet available in the literature.
In contrast to [10], our aim is not to develop an optimized method for a specific instance of (1), but rather to present a relatively simple algorithm that can evaluate (1) accurately, for a general class of f and g, without the need for expert user input or a priori analysis, even in the case of coalescing stationary points, thus addressing problems (P1) and (P2).Our specific focus in this paper is on the case where f is entire and g is a polynomial.The extension of our approach to more general cases where f and/or g have pole or branch point singularities is the subject of ongoing research.Necessarily, in aiming for generality and robustness we will sacrifice some efficiency.However, our method is designed to be rapidly convergent as N → ∞ with approximately ωindependent cost, and the fact that this is realised in practice is demonstrated by extensive numerical experiments in §5.
Our algorithm follows the basic principles of NSD, combining complex contour deformation with numerical quadrature.However, in contrast to standard NSD our algorithm does not trace SD contours directly from stationary points.Instead, stationary points are enclosed in a bounded "non-oscillatory region" within which the integrand is guaranteed to undergo at most a fixed number of oscillations.The original contour Γ is replaced by a "quasi-SD deformation" comprising a union of straight-line contours in the non-oscillatory region, for which numerical quadrature is straightforward, and SD contours outside the non-oscillatory region, on which standard NSD quadrature techniques can be applied.By excluding a neighbourhood of the stationary points from the region in which SD contours are traced, our algorithm avoids the problems mentioned in (P2) associated with stationary-point/stationary-point and/or stationarypoint/endpoint coalescence.This not only "uniformizes" the accuracy of our algorithm compared to standard NSD, but it also enables us to tackle the problem (P1) by automating the contour deformation step.For the latter, we first perform low accuracy SD contour tracing outside the non-oscillatory region to build a graph describing the global connections (via SD contours) between the endpoints of Γ, the different components of the non-oscillatory region, and the valleys at infinity, and then determine the quasi-SD deformation using a short-est path algorithm, before refining the accuracy of the SD contour tracing at the quadrature stage.
One other problem with standard NSD is that it typically degenerates as ω → 0, because the quadrature points diverge to infinity [7, §5.2.4].This issue has been addressed in the special case g(z) = z for bounded Γ in [2,6]; however, it remains an open problem for general g(z).Our algorithm is well-behaved in the limit as ω → 0 for general polynomial g(z), since it reduces to standard non-oscillatory quadrature for sufficiently small ω for any bounded Γ.
Our algorithm is implemented in the open-source Matlab code "PathFinder", available at github.com/AndrewGibbs/PathFinder[8].The basic user input to the code is a function handle for the amplitude f, the coefficients of the polynomial phase g, endpoints a and b (complex numbers, or angles in the case of infinite endpoints), the frequency parameter omega, and a parameter N controlling the number of quadrature points to be used.Approximating the integral (1) using PathFinder can be done with the following Matlab command: Here 'infcontour' is an optional input for which the user should supply a Boolean array Examples of PathFinder code will be given in §5.Advanced users can also adjust a small number of other tuning parameters, whose role will be discussed during the presentation of our algorithm.An outline of the paper is as follows.In §2 we provide a detailed description of our algorithm, first presenting an overview of the main steps, and then providing details of how each step is realised in PathFinder.In §3 we present some theoretical results underpinning our approach.In §4 we discuss some further implementation details, and in §5 we exhibit numerical results demonstrating the performance of our algorithm on a range of challenging integrals with large ω and complicated stationary point configurations.
We end this introduction by remarking that integrals with coalescing stationary points are of fundamental importance in numerous applications, including the study of optics and high frequency (short wavelength) acoustics, where they describe the wave field in the vicinity of geometrical singularities (or "catastrophes") in the classical ray-theoretic framework, Kelvin's celebrated ship-wave problem, and the theory of molecular collisions in quantum mechanics and theoretical chemistry.A catalogue of such integrals, along with links to relevant literature, can be found in [1, §36].In §5.5 we show how PathFinder can be applied to accurately calculate these types of integrals.
We start with some definitions and basic facts.Let for some J ∈ N, J ≥ 1, and α j ∈ C, j = 0, . . ., J, with α J = 0. Then g has at most J − 1 stationary points, which are the solutions of We denote the set of all stationary points by P stat .We define the valleys at infinity to be the sectors of angular width π/J centred on the angles These have the property that if z = re iθ with θ ∈ (v − π/(2J), v + π/(2J)) for some v ∈ V then e iωg(z) → 0 as r → ∞.For each η ∈ C \ P stat there exists a unique SD contour γ η beginning at η and ending either at a stationary point ξ ∈ P stat or at a valley v ∈ V, on which g(z) = g(η) for z ∈ γ η (see, e.g., [4]).
We let P endp denote the set of finite endpoints of the integration contour Γ, which could have zero, one or two elements.We assume for now that any infinite endpoint of Γ is at one of the valleys v ∈ V; see §4.4 for extensions.
We now provide a high-level overview of our algorithm.The following steps will be explained in more detail in sections 2.1-2.6. 1. Compute the set of stationary points P stat (the solutions of (4)).
2. For each ξ ∈ P stat , select a radius r ξ > 0 for which the function e iωg (z)  is considered "non-oscillatory" on the closed ball Ω ξ of radius r ξ centred at ξ.These balls may overlap.However, if two balls overlap significantly, indicating near coalescence, one of the stationary points (along with its associated ball) is removed from the set P stat .This removal process continues recursively until no pair of balls is judged to overlap too much.
We call {Ω ξ } ξ∈Pstat the non-oscillatory balls, and their union the non-oscillatory region.
3. Find the local minima of |e iωg(z) | on the boundary of the non-oscillatory region Ω.We call these points exits, and denote by P exit the set of all exits.
4. For each η ∈ P exit ∪ (P endp \ Ω), trace the SD contour γ η from η, and determine whether We call points z ∈ ∂Ω determined in case (i) entrances, and denote by P entr the set of all entrances.
5. Construct a graph G with a vertex for each of the elements of P stat , P endp , P exit , P entr and V. Add edges between the vertices of G as follows: • For each ξ ∈ P stat , add an edge between each pair of elements of (P stat ∪ P endp ∪ P exit ∪ P entr ) ∩ Ω ξ .• For each pair ξ, ξ ∈ P stat , ξ = ξ , for which Ω ξ ∩ Ω ξ = ∅, add an edge between ξ and ξ , if not already added in the previous step.• For each η ∈ P exit ∪ (P endp \ Ω), add an edge between η and the entrance z ∈ P entr or the valley v ∈ V to which the SD contour γ η leads.
Find the shortest path (in the graph-theoretic sense) between the vertices corresponding to the endpoints of Γ.
6. Generate quadrature nodes and weights for the evaluation of each of the contour integrals corresponding to the edges in the shortest path.For an edge between two points in the non-oscillatory region, use a straight-line contour.For an edge between an exit or an endpoint of Γ to an entrance or a valley, use a refined version of the SD contour traced in step 4. The union of all the contours corresponding to the edges of the shortest path defines the "quasi-SD deformation" of the original integration contour.
Finally, use the quadrature nodes and weights to approximate the integrals over the contours in the quasi-SD deformation and sum them to obtain an approximation of the original integral (1).
In Figures and the parameters ω = 40, a = −1.5, b = 2, N = 10, using the default parameter set for PathFinder (see Table 4.1).For this choice of g there is one order 2 stationary point and 4 order one stationary points.In Figure 2.1 we plot these stationary points, along with their non-oscillatory balls, and the SD contours traced from the exits.Such plots can be generated in PathFinder by adding the optional 'plot' flag.The ball centred at the stationary point ξ = −i contains two entrances, reached by SD contours from the balls above.In Figure 2.2 we plot the graph G, using the optional PathFinder input flag 'plot graph'.This graph, in addition to edges corresponding to the SD contours shown in Figure 2.1, contains edges corresponding to contours between points in the nonoscillatory region, including connections within the two overlapping balls.The shortest path between a and b, which is highlighted with thick lines in Figure 2.2, corresponds to the quasi-SD deformation, the integral over which is equal to (1) by Cauchy's Theorem.The integral is discretised using N quadrature points on each contour in the quasi-SD deformation that makes a non-negligible contribution to the integral (see §2.6) -these points are plotted in Figure 2.1 in red.
The process of computing all the SD contours and the selection of a subset thereof via the shortest path algorithm addresses problem (P1).Surrounding stationary points by balls, and only tracing SD contours outside the balls, means that we avoid having to determine the local structure of the SD contours and compute integrals along them near stationary points, addressing problem (P2).
In the following subsections we provide further details of how we carry out the steps outlined above in PathFinder.

Step 1 -Computing stationary points
Computing the stationary points of g (the roots of g (z)) requires us to find the complex roots of the polynomial (4).In our implementation we compute stationary points using the Matlab roots command, which applies a companion matrix approach.We note that obtaining highly accurate values for the positions of stationary points is not critical to our algorithm, since the stationary points are enclosed in the non-oscillatory region and we never trace SD contours from them.Indeed, the difficulty in distinguishing numerically between multiple roots and roots of higher order contributes to the motivation for considering such nonoscillatory regions.

Step 2 -Defining the non-oscillatory region
The non-oscillatory region Ω was defined in (6) to be a union of balls centred at the elements of P stat .We choose the radii of the balls as follows.First fix some user-defined constant C ball > 0.Then, given ξ ∈ P stat , define This definition enforces an upper bound on the number of oscillations within each ball.Accordingly, the region Ω shrinks to the stationary points as ω → ∞ and expands to fill the whole complex plane as ω → 0.
In our implementation we approximate r ξ numerically as follows.Let N ball ∈ N be a user-defined parameter.For each n ∈ {1, . . ., N ball } we consider the ray {z = ξ + re i2πn/N ball , r > 0}, and compute the smallest positive root r n > 0 of the function u n (r ball , which is a polynomial in r of degree 2J.For this root-finding problem we use the Matlab roots command; in case this command produces no positive real roots (because of stability issues) we resort to a bisection approach instead.We then take as our approximation to r ξ the positive number max n∈{1,...,N ball } r n .
When elements of P stat are close it is natural to amalgamate their respective non-oscillatory balls.To do this systematically we adopt an iterative approach.Let δ ball > 0 be a user-defined parameter.
• For each pair ξ 1 , ξ 2 ∈ P stat compute   4.1).Here we observe stationary points (black stars) surrounded by balls (grey), SD contours traced from exits and finite endpoints (black lines), and quadrature points allocated along the appropriate contours in the quasi-SD deformation (red points).The "region of no return" (see §3.2) around the valleys at infinity is also shaded grey.• If min ξ1,ξ2 d ξ1,ξ2 < δ ball let ξ 1 , ξ 2 be a pair realising the minimum.Remove from P stat whichever of ξ 1 , ξ 2 has the smaller associated ball radius (or choose arbitrarily between them if r ξ1 = r ξ2 ).
• Repeat the previous step until either min ξ1,ξ2 d ξ1,ξ2 ≥ δ ball , or there is only one element of P stat remaining.

Step 3 -Determining the exits
The exits associated with each ξ ∈ P stat are defined to be the local minima on For each ξ ∈ P stat the function − g(z) restricted to the ∂Ω ξ is a trigonometric polynomial.Using this fact, in our implementation we determine the local minima of − g(z) on ∂Ω ξ by finding the roots of the derivative of − g(z) in the angular direction (which is also a trigonometric polynomial) by the companion matrix approach of [5, §2.2], and keep only the real roots corresponding to local minima.We discard all those minima corresponding to points inside ξ ∈Pstat,ξ =ξ Ω • ξ , and add the remaining minima to the set P exit .

Step 4 -Tracing the SD contours
Given η ∈ P exit ∪ (P endp \ Ω), the SD contour γ η beginning at η is the unique curve on which g(z) is constant, with − g(z) decreasing along γ η .It can be parametrized in terms of a parameter p ≥ 0 as z = h η (p), where h η (p) is defined implicitly by Differentiating ( 9) with respect to p gives which is a first order ODE initial value problem for h η (p).By solving (10) numerically one can trace the contour γ η until it either (i) enters the nonoscillatory region Ω, or (ii) one can be sure that it will tend to a valley v ∈ V, without entering Ω.For (ii) we appeal to the theoretical result in Theorem 3.3, which provides a "region of no return" R v associated with each valley v ∈ V for which it is guaranteed that if an SD contour enters R v it will never leave R v , and will converge to v.
Staying away from stationary points ensures that the factor 1/g in the righthand side of (10) does not get too large.
In our implementation we trace the SD contour using a predictor-corrector approach, combining a forward Euler step for (10) and a Newton iteration for (9), to generate approximations h where the total number of steps n max is determined as part of the algorithm, as discussed below.
As the initial value we take h we first apply a forward Euler step for the ODE (10), with adaptive step length where δ ODE ∈ (0, 1) is a user-specified parameter.The first argument of the minimum is included to ensure stability of the solver -note that F (h) = − ig (h) (g (h)) 2 and we might expect instability if the local step length were as large as 2/|F (h The second argument is included to ensure that the solver "slows down" as it approaches the non-oscillatory region Ω, so that we can detect whether γ η enters Ω or not.To ensure that |h This also ensures that h (n+1) η remains far enough from P stat , so that (10) doesn't get too large.
After each forward Euler step, we correct h by running a Newton iteration to enforce (9) (with p = p n+1 fixed), until the Newton step size , P stat ), for some userspecified tolerance δ coarse > 0.
We repeat this process for n = 0, 1, 2, . . .until either ∈ Ω ξ for some ξ ∈ P stat , in which case we add z = h to the set P entr of entrances.Note that in general the point z = h (n) η will lie inside Ω • ξ rather than on ∂Ω ξ , but will be closer to ∂Ω ξ the smaller δ ODE is; or in which case, by Theorem 3.3, γ η converges to the valley v.Here the "region of no return" R v is defined by where and, for r > 0 and θ ∈ (0, π/(2J)), For further explanation of the meaning of R v see §3.2 below.
A necessary condition for a point z to lie in R v is that |z| ≥ r * , where r * > 0 is the unique positive solution of the polynomial equation G(r * , π/(4J)) = 0, i.e. the solution of Having found r * once and for all (using the Matlab roots commnad), to check that a point z lies in R v we first check that |z| ≥ r * .If so, we then check that | arg z − v| 2π < π/(2J).If so, we then check that G(|z|, | arg z − v| 2π ) > 0. The point of introducing r * is so that we don't compute G(|z|, | arg z − v| 2π ) unless absolutely necessary.
In either case, tracing of the SD contour stops and we set n max = n for this contour.

Step 5 -Finding the shortest path
The construction of the graph G requires no further explanation.To find the shortest path in G between the endpoints of the original contour Γ we apply the standard Dijkstra shortest path algorithm [13, §10.6.2].

Step 6 -Evaluating the contour integrals
The quasi-SD contour deformation corresponding to the graph-theoretic shortest path between the endpoints of G calculated during step 5 involves integrals over three types of contour: Type 1: Straight line contours between points in the non-oscillatory region; Type 2: Infinite SD contours from exits/endpoints to valleys; Type 3: Finite SD contours from exits/endpoints to entrances.Some of these contours will make a larger contribution to the value of the original integral (1) than others.It is natural to neglect contours that make a negligibly small contribution.In our implementation, we only compute the contribution from a contour γ in the quasi-SD deformation if at least one of the finite endpoints η of γ satisfies |e iωg(η) |/M > δ quad , where δ quad ≥ 0 is a small, user-specified parameter and where the maximum is taken over all ξ ∈ P stat ∪ P endp ∪ P exit appearing in the shortest path corresponding to the quasi-SD deformation.
In our implementation, for Type 1 contours we use Gauss-Legendre quadrature, for Type 2 contours we use either Gauss-Laguerre quadrature (which is the default choice in PathFinder) or truncated Gauss-Legendre quadrature, and for Type 3 contours we use (possibly truncated) Gauss-Legendre quadrature, as detailed below.By default our implementation uses the same number N of quadrature points on each contour in the quasi-SD deformation whose contribution we compute, regardless of the type of integral (we comment on this in §3.3.4).Accordingly, if N cont is the number of these contours then the total number of quadrature points used in the algorithm, N tot , is given by

Evaluation of integrals over Type 1 contours
Let z 0 , z 1 ∈ Ω, and let γ be the straight-line contour in C starting at z 0 and ending at z 1 , parametrized by Given N ∈ N, let t Leg m and w Leg m , for m = 1, . . ., N , denote the nodes and weights for standard N -point Gauss-Legendre quadrature on [−1, 1].Our quadrature approximation to the integral over γ is then:
By tracing contours outside of Ω, the contours remain a positive distance from P stat .This ensures that f is analytic in a complex neighbourhood of [0, ∞).By default, PathFinder evaluates the integral on the right-hand side of ( 17) by Gauss-Laguerre quadrature.Let t Lag m and w Lag m , for n = 1, . . ., N , denote the standard Gauss-Laguerre nodes and weights on [0, ∞).Our quadrature approximation to the integral over γ is then: To evaluate f (t Lag m ) we need accurate computations of h η (t Lag m /ω) for m = 1, . . ., N .For this, for each m we run a Newton iteration on (9) with p = t Lag m /ω fixed, until the magnitude of the increment is smaller than a user-specified tolerance δ fine > 0. Typically we take δ fine to be considerably smaller than the tolerance δ coarse used in the Newton iteration in step 4, since when carrying out quadrature we require higher accuracy in our approximation of the SD contour than is required for determining the global structure of the quasi-SD deformation in step 4. As the initial guess for the Newton method we use a piecewise linear interpolant of the points {(p 0 , h )} computed in step 4, where n max denotes the total number of steps taken in the ODE solve in step 4 before the contour tracing algorithm terminated.If p nmax < t Lag N /ω then before running the Newton iteration we first need to restart the contour tracing algorithm of step 4 to extend the SD contour until p nmax ≥ t Lag N /ω, so that there are points to interpolate between.
As an alternative, one can evaluate the integral over a Type 2 contour using truncated Gauss-Legendre quadrature, as suggested in [16].To activate this alternative in PathFinder one should add the optional input 'inf quad rule', 'legendre'.In this case we truncate the integral to for some P > 0, then apply Gauss-Legendre quadrature on [0, P ], to obtain the approximation where we compute h η (z [0,P ] (t Leg m /ω)) (which is required for the evaluation of f (z [0,P ] (t Leg m ))) by the same Newton iteration discussed above for h η (t Lag m /ω).For the truncation point P we take where which describes the point at which the magnitude of the exponential part of the integrand drops below δ quad times its maximum value M on the quasi-SD deformation.

Evaluation of integrals over Type 3 contours
Let η ∈ P exit ∪ (P endp \ Ω) be such that the SD contour γ from η leads to an entrance z ∈ P entr .In this case we apply (possibly truncated) Gauss-Legendre quadrature as in formulas ( 19) and ( 20), but now with where p nmax is defined as in §2.4 and L is defined as in §2.6.2.
In the case where the minimum is attained by p nmax /ω, so that the whole contour is considered, a potential inconsistency arises, because the application of the higher accuracy Newton iteration described in §2.6.2 for the calculation of h η (z [0,P ] (t Leg m /ω)) corresponds implicitly to a slight shifting of the endpoint of the contour γ away from the entrance z = h (nmax) η added to the graph G in step 5. To avoid this inconsistency, in our implementation, in step 4, whenever the contour tracing terminates in case (i), we run a Newton iteration on the final point h (nmax) η with the high accuracy tolerance δ fine , before adding it to the list of entrances P entr .Note that this may mean that h (nmax) η lies very slightly outside Ω.

Theoretical results
In this section we collect some theoretical results that motivate the design of our algorithm.

Removal of stationary points
In §2.2 we described our algorithm for removing stationary points from the set P stat when they are close.When removing stationary points and their associated non-oscillatory balls, we need to ensure that the removed stationary points still lie inside one of the remaining non-oscillatory balls, so that we don't encounter any stationary points along the trajectory in our ODE solve for the SD contour tracing (see the discussion in §3.3.2 below).In this section we provide a sufficient condition on the parameter δ ball for this to be guaranteed.Proposition 3.1.Suppose that in the removal algorithm of §2.2, n stationary points have been removed from P stat .Then for any stationary point ξ that was removed, there exists ξ ∈ P stat such that |ξ − ξ | ≤ nδ ball r ξ .
Proof.We proceed by induction on n.The result is trivially true for n = 0. Assume that it is true after the removal of n points, and suppose that the (n+1)st point is now to be removed.Let ξ 1 , ξ 2 denote the pair of points selected as realising min ξ1,ξ2 d ξ1,ξ2 , and without loss of generality suppose that ξ 2 is the point to be removed (so that r ξ1 ≥ r ξ2 ).Then |ξ 2 −ξ 1 | ≤ δ ball r ξ1 ≤ (n+1)δ ball r ξ1 , so the claimed property holds for ξ 2 .Furthermore, by the inductive hypothesis, for each point ξ previously removed, there exists ξ ∈ P stat such that |ξ − ξ | ≤ nδ ball r ξ .If ξ = ξ 2 then ξ will still be present in P stat after the removal of ξ 2 , and |ξ − ξ | ≤ nδ ball r ξ ≤ (n + 1)δ ball r ξ .On the other hand, if ξ = ξ 2 then ξ will not be present in P stat after the removal of ξ 2 , but ξ 1 will be, and by the triangle inequality completing the inductive step.

Region of no return for SD contours
The following result establishes a region of no return: once an SD contour enters this region, we can say with certainty which valley it will converge to.The idea behind this result is that in the region of no return the highest degree term α J z J of the polynomial g is sufficiently dominant over the lower degree terms that the SD contours inside the region converge to the same valley as those corresponding to the monomial phase α J z J .Theorem 3.3 (Region of no return).Let g, V and R v , for v ∈ V, be as in (3), ( 5) and (11).The regions R v , v ∈ V, contain no stationary points of g.Furthermore, if an SD contour enters R v for some v ∈ V, it never leaves R v .
Proof.That R v contains no stationary points follows because if G(r, θ) > 0 then and R > 0 we define the sector with | • | 2π defined as in (12).We also define the function which for each fixed θ is a polynomial in R of degree J − 1.
We claim that if θ ∈ (0, π/(2J)) and G(R, θ ) > 0, then if an SD contour enters S v (R, θ ) it never leaves S v (R, θ ).To prove this, we show that if an SD contour intersects ∂S v (R, θ ) then the direction of descent always points into S v (R, θ ).Since ∂S v (R, θ ) is the union of the sets ijα j e ijθ r j−1 , and [iα J e iJθ ] = [iα J e iJ(v±θ ) ] = ∓|α J | sin(Jθ ) we have that The function φ(r) has the property that if R > 0 and φ(R) > 0 then φ(r) > 0 for all r ≥ R. To see this, note that and that the term in brackets is a strictly decreasing function of r, which tends to −∞ as r → 0 and to J|α J | sin(Jθ ) > 0 as r → ∞.Hence a sufficient condition for (26) is that φ(R) > 0, i.e.

Quadrature error
In §2.2 we defined the non-oscillatory region as a union of balls on which the exponential e iωg(z) undergoes a bounded number of oscillations.Here we show that the definition (8) strikes a balance between the accuracy of our quadrature approximations to the integrals outside and inside this region.

Quadrature in the non-oscillatory region
The Type 1 straight line contour integrals between points in the non-oscillatory region are evaluated using Gauss-Legendre quadrature, as detailed in §2.6.1.To assess the accuracy of this we note the following theorem, which is a simple consequence of the standard error analysis presented in [15,Chap. 19].
Theorem 3.4.Let z 0 , z 1 ∈ C. Suppose that γ is a straight-line contour in C starting at z 0 and ending at z 1 and that there exists ρ > 0, C > 0 and ξ * ∈ C such that f is analytic and bounded in z [z0,z1] (B ρ ), where B ρ is a standard Bernstein ellipse (relative to [−1, 1]) and z [z0,z1] is defined as in (15), and Let I and Q denote the left-and right-hand sides of (16), respectively.Then, for some C > 0, depending only on ρ, Proof.Noting that the result follows from [15,Thm 19.3].
Theorem 3.4 motivates the definition of the non-oscillatory region in (8).Indeed, if the assumptions of Theorem 3.4 hold with ρ and C independent of ω then the bound (30) guarantees ω-independent exponential convergence for ω bounded away from zero.However, even when ( 8) is satisfied, the relationship between ξ * , ρ, C and ω is beyond our control in general because the ellipse may extend beyond the non-oscillatory region, so that C > C ball .Thus we cannot control the factor e C entirely based on condition (8).
Still, the bound (30) shows that the quadrature error decreases with increasing N .The precise rate of decrease depends on a balance between the decay of ρ −2N and the growth of e C and f L ∞ (z [z 0 ,z 1 ] (Bρ)) for increasing ρ.We quantify this in the special case of monomial phase in §3.3.3.

Quadrature for the SD contours
For Type 2 or Type 3 integrals along SD contours we use either Gauss-Laguerre or (possibly truncated) Gauss-Legendre quadrature, as detailed in §2.6.2 and §2.6.3.We expect these rules to converge rapidly to the true value of the integral as the number of quadrature points N tends to infinity, provided that the integrand is analytic and bounded in a suitable region of the complex p plane.
Theorem 3.5.Suppose that f is analytic inside and on the parabola P ρ := {z ∈ C : √ −z = ρ} for some ρ > 0, where the branch cut is along the positive real axis and √ −z is real and positive on the negative real axis, that f grows at most algebraically as z → ∞ inside the parabola, and that the integral is finite.Let I and Q denote the left-and right-hand sides of (18), respectively.Then This result implies that our Gauss-Laguerre quadrature approximation should converge root-exponentially as N → ∞, provided that f is sufficiently wellbehaved at infinity.The presence of singularities in the complex p-plane limits the size of ρ, and hence the convergence rate.We know from ( 17) that our integrand is singular at points p ∈ C where g (h η (p/ω)) = 0, i.e.where h η (p/ω) = ξ for some stationary point ξ.Since we only trace SD contours outside the nonoscillatory region (which contains the stationary points), we know that there cannot be singularities on the SD contour itself.If the start point η lies on an SD contour emanating from a stationary point ξ then we expect there to be a singularity in the p-plane at p = ω [g(ξ) − g(η)] < 0. We show in §3.3.3 that in the special case of monomial phase this singularity lies at p = −C ball , which implies root-exponential convergence independent of ω for ω bounded away from zero.Determining the locations of the other possible singularities in the complex p-plane is more challenging, since it involves study of the (multivalued) inverse of g.We leave further theoretical investigation of this to future work.

Results for monomial phase
It is instructive to consider the special case of a monomial phase g(z) = z J for some J ∈ N. In this case there is a single stationary point of order J − 1 at ξ = 0, and g(0) = 0. Following the prescription (8), we obtain a ball radius We first consider a Type 1 integral in the non-oscillatory region.For simplicity we choose f (z) ≡ 1.Specifically, we consider the evaluation of the integral r0e iθ 0 e iωg(z) dz, for some θ ∈ [0, 2π].Taking ξ * = 0, we can apply Theorem 3.4 with any ρ > 1, and the resulting scaled and translated Bernstein ellipse surrounding [0, r 0 e iθ ] is tightly contained in the disc |z| ≤ sr 0 , where ρ and s are related by Hence condition (29) is satisfied, independently of θ, with which is independent of ω but dependent on J.When s is large, we have ρ ≈ 4s, and in this regime the error bound provided by (30) for Gauss-Legendre quadrature is approximately proportional to As a function of s, with J and N fixed, this quantity is minimised where its s-derivative vanishes, which occurs where Accordingly, the error bound is approximately proportional to Thus we expect super-exponential convergence as N → ∞ for fixed J.However, we expect the convergence to be slower the larger J is.
The integrand has a branch point at p = −C ball , but we note that the distance between the branch point and the positive real p-axis equals C ball , which is independent of both ω and J.
For truncated Gauss-Legendre the relevant theory can be found in [15,Chap. 19] (and see also [16]).Due to the branch point at p = −C ball , as N → ∞ we obtain exponential convergence to the integral over the interval [0, P ], where P is given by either (21) or (23).In the case where P = L, by the definition of L in (22), we expect the truncation error to have relative order δ quad .

Number and distribution of quadrature points
PathFinder uses a fixed number N of quadrature points on each contributing contour, and that number is the same both for integrals within and outside the non-oscillatory region, i.e., for Gauss-Legendre and Gauss-Laguerre quadrature.Thus, increasing the single parameter N provides a way of uniformly improving accuracy.
The theoretical results in this section (specifically, Theorems 3.4 and 3.5) imply that the precise rate of improvement with respect to N depends on the type of integral being approximated.They suggest even that a different strategy for the distribution of quadrature points may be superior.Indeed, exponential convergence of Gauss-Legendre for Type 1 integrals in the non-oscillatory region is not balanced with root-exponential convergence of Gauss-Laguerre for Type 2 integrals outside.Similarly, convergence rates of Gauss-Laguerre and truncated Gauss-Legendre outside the non-oscillatory region are different.Our choice of a fixed parameter N is inspired on the one hand by simplicity, and on the other hand by the lack of robust methods to optimize parameters in alternative schemes.For example, we have shown in §3.3.3 that the convergence rate of Gauss-Legendre for Type 1 integrals may depend on the order of nearby stationary points.While this can be quantified precisely for the case of monomial phase, it is not at all clear how to generalise this analysis when a cluster of multiple stationary points is present.Hence, stationary point order is a quantity that we deliberately do not explicitly compute, estimate or rely on in any way.Implicitly, of course, it plays a big role, and it does so mainly via the definition of the ball of the radius in (8).
The main practical benefit of the theoretical analysis of quadrature error in this section is the guarantee that N is a robust parameter for improving accuracy.Concerning possible future improvements, rather than attempting to optimize the quadrature point distribution a priori, we believe a more promising development would be the ability to invoke standard adaptive quadrature schemes along the contours for a given function f .However, it should be borne in mind that quadrature forms just one step in our algorithm, and that the other steps (particularly the SD path tracing) incur a non-negligible cost overhead, that should also be considered when trying to further optimize performance.

The case J = 1
In the case J = 1 (linear phase) there are no stationary points, and our algorithm simplifies dramatically.Furthermore, the SD contours are simply parallel straight lines in the direction of the single valley at angle π/2 − arg(α 1 ), and there is no need to trace them numerically.Hence when J = 1 PathFinder skips the ODE contour tracing step and exploits the exact characterization of the SD contours mentioned above.

Specifying infinite endpoints
In the description of our algorithm in §2 we made the assumption that any infinite endpoint of the contour Γ should be at a valley v ∈ V. PathFinder is actually more flexible than this.The user is permitted to specify an infinite endpoint at any θ ∈ [v − π/(2J), v + π/(2J)] and the code will automatically adjust this to equal v.The case θ = v ± π/(2J) is delicate because the highest order term in the phase does not provide exponential decay along the contour.Nonetheless, we include it, because in applications one often encounters this case, with the integral converging conditionally (under appropriate assumptions on f ) and the contour deformation to v being justified by Jordan's Lemma.

Numerical results
In this section we present numerical results illustrating the performance of our algorithm and its implementation in PathFinder.All results in this section were produced using PathFinder Version 1.0 [8].

A "generic" example
We begin by illustrating the performance of PathFinder on the integral where, to convey the message that our approach is applicable to truly "generic" amplitudes and polynomial phase functions, the coefficients of f and g are chosen to be the first 5 digits of e and the first 10 digits of π, respectively.This can be approximated by PathFinder via the Matlab code (cf.( 2)) PathFinder(-1,1,@(z) 2*z.^4+7*z.^3+z.^2+8*z+2,...
[3 1 4 1 5 9 2 6 5 3],omega,N) In Figure 5.1 we plot the quasi-SD deformations and quadrature point distributions (using the PathFinder 'plot' option) for (32) for ω ∈ {0.01, 1, 5, 50} and N = 10.As explained in §2.2, for smaller ω the non-oscillatory balls are larger, and can overlap, while for larger ω they shrink around the stationary points.In more detail, in Figure 5.1(a) (ω = 0.01), ω is small enough that both endpoints are inside the same non-oscillatory ball.Hence the integral is treated as non-oscillatory and is approximated by Gauss-Legendre quadrature along a single straight-line contour.In Figure 5.1(b) (ω = 1), ω is still small enough that many of the balls overlap, and the quasi-SD deformation comprises two SD contours (one from an exit and one from an endpoint) plus four straightline contours in the non-oscillatory region.In Figure 5.1(c) (ω = 5), ω is large enough that only two balls overlap, and the quasi-SD deformation comprises five SD contours (two from endpoints, two from exits to valleys, and one from an exit to an entrance), plus four straight-line contours in the non-oscillatory region.Finally, in Figure 5.1(d) (ω = 50), ω is so large that none of the balls overlap, and the quasi-SD deformation comprises eight contributing SD contours (two from endpoints and six from exits to valleys), plus three straight-line contours in the non-oscillatory region.However, in this case the two SD contours and one straight-line contour associated with the stationary point near 0.2 + 0.5i are judged to make a negligible contribution to the integral, so are not  assigned any quadrature points.We emphasize that this intricate behaviour is fully automated, with no expert input required from the user.In Figure 5.2(a) we plot the error in the PathFinder approximation of (32), compared to reference values computed using the Julia QuadGK package when ω < 500, and using PathFinder with N = 500 when ω ≥ 500.For fixed ω we observe rapid convergence as N → ∞, at a rate that appears independent of ω.In Figure 5.2(b) we show the associated computation times, which remain bounded as ω increases.

Coalescence and the Airy function
which is of the form (1) with f ≡ 1, ω = 1 and g(z; x) = −i(z 3 /3 − xz).Up to a change of variable this is the same example for which, as mentioned in §1, a bespoke, complex Gaussian quadrature rule was developed in [10].Ai can be approximated by PathFinder via the Matlab code Ai = @(x) PathFinder(-pi/3,pi/3,[],... where the input [] for f indicates that f ≡ 1.In Figure 5.3 we plot the quasi-SD deformations, along with the distribution of quadrature points for N = 20, for the evaluation of Ai(x) at x ∈ {−5, −1, −0.5, 0, 5}.Here one observes in detail how our algorithm deals with stationary point coalescence, as the non-oscillatory balls overlap, merge, then split.In Figure 5.3(a) the quasi-SD deformation comprises four SD contours from exits, plus two straight-line contours inside balls (which do not go via stationary points).In Figure 5.3(b) the balls overlap and this changes to two SD  contours from exits plus three straight-line contours inside balls (which go via both stationary points).In Figure 5.3(c) the balls overlap enough that both stationary points are contained in both balls, so we get two SD contours from exits plus just two straight-line contours inside balls (which go via only one of the stationary points).In Figure 5.3(d) the balls have merged completely and in addition to the two SD contours from exits there is just one straight-line contour inside a ball (which does not go via the stationary point).In Figure 5.3(e) the balls have split again, but we see the same deformation structure as in Figure 5. 3(d).Again, we emphasize that these calculations are fully automated.In Figure 5.4 we show the accuracy of the PathFinder approximation for this example as a function of x ∈ [−10, 4], for different N .Our reference is the built-in Matlab command airy.We note that between x = −3 and x = 0 the error for the smaller values of N undergoes some jumps.These are due to the fact that near stationary point coalescence the topology of the quasi-SD deformation, the number of contours constituting it, and hence the total number of quadrature points along it (recall ( 14)), all change discontinuously as a function of x (as illustrated in Figure 5.3).However, as N increases we see a clear, approximately exponential decrease in the error, and, although the rate of decrease depends slightly on x (because of the factors mentioned above), for N = 30 we achieve approximately 10 −13 error uniformly across the interval.

A high order stationary point -comparison with Mathematica's implementation of Levin quadrature
We now consider the integral which has a stationary point of order 8 at the origin.The integral (34) can be approximated by PathFinder via the command PathFinder(-1,1,@(z) sin(z),[1 0 0 0 0 0 0 0 0 0],omega,N) Figure 5.5 shows the quasi-SD deformation and quadrature point distribution obtained by PathFinder for ω = 100, 000 and N = 50.There are small contributions from the endpoints, but the main contribution comes from the ball containing the stationary point.
In the Mathematica documentation [18, pp75-86], it is stated that oscillatory integrals with monomial phase functions such as (34) can be evaluated efficiently using the built-in Mathematica function NIntegrate, via its implementation of Levin quadrature (which is described, e.g., in [7, §3.3]).To do this one can use the Mathematica command:

Coalescence to a high order stationary point
We now investigate the robustness of our algorithm in the presence of a large number of coalescing stationary points.Specifically, we consider the integral 1 −1 e iω(z 7 /7−r 6 z) dz, where r ≥ 0 is a parameter controlling the coalescence.For r > 0 there are 6 stationary points with |ξ| = r, namely the solutions of ξ 6 = r 6 , and for r = 0 there is a single stationary point of order 6.To evaluate this integral in PathFinder for a given r, one can use the command PathFinder(-1,1,[],[1/7 0 0 0 0 0 -r^6 0],omega,N) In Figure 5.7 we plot the quasi-SD deformations and quadrature point distributions for a some different r values, showing how the balls first intersect and then merge as r → 0. In Figure 5.8 we show convergence (with respect to a PathFinder reference with N = 500) and CPU times (averaged over 100 runs) for fixed r = 0.01.We see that both the error and the CPU time are essentially independent of ω in this case.In Figure 5.9 we plot errors for two fixed N values N = 10, 50, as a function of r.We observe that as r → 0, the error stays bounded.For N = 10 the error jumps up between r = 10 −3 and r = 10 −2 , at a point depending on ω.This represents the point at which the balls around the stationary points merge, resulting in a reduction of N tot , and hence a reduction in accuracy.But after this point we observe no further reduction in accuracy as r → 0. We remark that for sufficiently small r > 0 the six stationary points are numerically indistinguishable, but this isn't a problem for our algorithm because in that case the problem will be treated identically to that of a monomial phase.

Canonical cuspoid integrals and their generalisations
In this section we show how our algorithm can be applied to the computation of some of the canonical integrals catalogued by Berry and Howls in [1, §36], which, as mentioned already in §1, are of fundamental importance in numerous application areas including optics, acoustics and quantum mechanics.In this context, our algorithm is related to that of [11], where an adaptive contour deformation approach was applied to evaluate the cuspoid integrals considered in §5.5.1.The algorithm in [11] is similar in spirit to our approach, in that it deforms the integration contour so that it terminates in valleys at infinity, and splits the contour into portions close to stationary points plus portions away from stationary points.However, in contrast to our approach, the algorithm in [11] does not attempt to trace SD contours, and hence is susceptible to rounding errors associated with the "violent" behaviour of the exponential factor e iωg(z) when one is not on a true SD contour -see [11, §2].Furthermore, while the algorithm in [11] was specialised to the case of integration over the real line, our algorithm can handle much more general contours, as we illustrate in §5.5.2.

Cuspoid integrals
The so-called "cuspoid integrals" listed in [1, §36.2.4] are all of the form (1) with polynomial phase g and unit amplitude f ≡ 1, unit frequency ω = 1, and integration along the real line.Our algorithm is ideally suited to the evaluation of these integrals, and to demonstrate this we compute two of them.In the notation of [1, §36], we consider the cusp catastrophe integral where P is the Pearcey function, and the swallowtail catastrophe integral Both exhibit coalescence of stationary points on certain algebraic varieties (see [1, §36.5(ii)]) on which both the first and second derivatives of the phase function vanish.In the case of (36) this occurs when and for (37) this occurs when  1]. Computation times on a small desktop computer (Intel i7-4790, 32GB RAM) were less than a minute for the cusp (which required the computation of 10000 instances of (36), averaging 0.005s per instance) and less than an hour for the swallowtail (which required 250000 instances of (37), averaging 0.01s per instance).

Generalisations
In [9] the authors considered a family of generalisations of certain canonical cuspoid integrals, with integration no longer over the real line, but rather over a complex contour starting and ending at valleys at infinity, and possibly with a non-unit amplitude function.
A specific aim of [9] was to investigate the relevance of such integrals to the study of the so-called "inflection point problem", a canonical problem in wave scattering originally introduced over 50 years ago by Popov in [12].This problem, which remains unsolved in closed form, concerns two-dimensional timeharmonic wave propagation near a boundary with an inflection point, and seeks a solution for the wave field near the inflection point that describes the transition from an incoming "whispering gallery wave" supported on the concave portion of the boundary, to outgoing "creeping waves" along the convex portion of the boundary, along with a scattered "searchlight" beam (for details and further references see [14]).
In this context, in [9, §3.3] the authors studied the family of integrals where f (t) is some amplitude to be specified, and Γ ij denotes any contour from valley v i to valley v j , where v j := (2(j − 1) + 1/2)π/5, j = 1, . . ., 5.These integrals have stationary point coalescence on the cubic curve y + 4x 3 /27 = 0, which suggests that, by appropriately choosing f and Γ ij , they might exhibit certain features of the solution of the inflection point problem.Indeed, in [9, §4] it was shown that as x → −∞ near the cubic curve, the integral A 32 has the character of an incoming whispering gallery type wave, and that, as x → +∞ near the cubic curve, the integral A 52 has the character of an outgoing creeping wave.However, plots of the resulting fields could not be presented in [9] due to the lack of a suitable numerical evaluation method and implementation.Using PathFinder we are able to remedy this.In Figures 5.11 We only plot A 52 above the cubic curve y + 4x 3 /27 = 0, because below this curve A 52 becomes exponentially large (cf.[9, Fig. 12(i)]).In Figures 5.11Here one observes the predicted incoming whispering gallery type behaviour of A 32 near the top of Figure 5.11(c) between x 0 = −2 and x 0 = −1, with oscillations giving way to an exponentially small field in the caustic shadow, and the predicted creeping wave type behaviour of A 52 near the bottom of Figure 5.11(d) between x 0 = 1 and x 0 = 2, with waves propagating along the cubic curve, shedding rays tangentially.In ongoing and future studies we plan to use PathFinder to further investigate the properties of integrals of the form (40), and generalisations involving different choices of f and higher degree phase functions (see [9]), which we hope may shed new light on the inflection point problem and related problems in high frequency wave propagation.
(a) SD contours and quasi-SD deformation (b) Zoom of boxed region near the origin

Figure 2 . 1 :
Figure 2.1: Output of algorithm when applied with phase (7), ω = 40, a = −1.5, b = 2, N = 10, and the default parameter set (see Table4.1).Here we observe stationary points (black stars) surrounded by balls (grey), SD contours traced from exits and finite endpoints (black lines), and quadrature points allocated along the appropriate contours in the quasi-SD deformation (red points).The "region of no return" (see §3.2) around the valleys at infinity is also shaded grey.8 (a) Graph G (b) Zoom of the graph G

Figure 2 . 2 :
Figure 2.2: The graph G corresponding to the problem considered in Figure 2.1.The thick line represents the shortest path between the endpoints, which in this case are both finite.The balls (shaded grey) are included for ease of comparison with Figure 2.1.The lower figure zooms in on the centre of the upper figure, showing the multiple edges that are constructed inside the balls.

Figure 5 . 4 :
Figure 5.4: Accuracy of PathFinder approximation of Ai(x) for different N .
NIntegrate[Sin[x]Exp[omega*I*x^9],{x,-1,1}, Method->{"LevinRule","Kernel"->Exp[omega*I*x^9]}] In Figure 5.6(a) we show a plot of the relative accuracy of our PathFinder approximation, compared to the Mathematica approximation (using the default settings), as a function of ω, for different N values.For all three N values the accuracy of our approximation is approximately uniform in ω, and for N = 50 our approximation agrees with Mathematica's to approx 13 digits.In Figure 5.6(b) we report the corresponding computation times (averaged over 100 (a) Accuracy (compared to Mathematica) (b) Timings

400x 3 −Figure 5 .
Figure 5.10 shows plots of the magnitude of (36) and (37) (the latter over the plane z = −7.5),computed using PathFinder with the default settings and N = 50.The plots agree qualitatively with those presented in [1, Figs 36.3.1 & 36.6.5],and, for (36), agree quantitatively (to all five decimal places presented) with the values presented in [11, Table1].Computation times on a small desktop computer (Intel i7-4790, 32GB RAM) were less than a minute for the cusp (which required the computation of 10000 instances of (36), averaging 0.005s per instance) and less than an hour for the swallowtail (which required 250000 instances of (37), averaging 0.01s per instance).