An asymptotic expansion of the hyberbolic umbilic catastrophe integral

We obtain an asymptotic expansion of the hyperbolic umbilic catastrophe integral Ψ(H)(x,y,z):=∫-∞∞∫-∞∞exp(i(s3+t3+zst\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Psi ^{(H)}(x,y,z):= \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\exp (i(s^3+t^3+zst$$\end{document}+yt+xs))dsdt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+yt+xs))\mathrm{{d}}s\,\mathrm{{d}}t$$\end{document} for large values of |x| and bounded values of |y| and |z|. The expansion is given in terms of Airy functions and inverse powers of x. There is only one Stokes ray at argx=π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\arg x=\pi $$\end{document}. We use the modified saddle point method introduced in (López et al. J Math Anal Appl 354(1):347–359, 2009). The accuracy and the asymptotic character of the approximations are illustrated with numerical experiments.


Introduction
Catastrophe theory studies and classifies phenomena characterized by sudden shifts in behavior arising from small changes in circumstances [24]. Catastrophe integrals This [5,Chap. 36] appear in the description of many physical phenomena, specially those related to wave propagation (see [10,11,15,25,27,28] or references therein for a detailed information). They also have an important mathematical application in the uniform asymptotic approximation of oscillatory integrals [23]. They are oscillatory integrals with several nearly coincident stationary phase or saddle points and they are classified according to the number of free independent parameters that describe the type of singularities arising in catastrophe theory and dimensionality of the integral. The canonical integrals associated to cuspoid catastrophes (simple integrals) are defined in the form [5,Eq. 36.2.4] K (x, y, z) := ∞ −∞ exp (i K (t; x)) dt, with K (t; x) := t K +2 + K m=1 x m t m ([5, Eq. 36.2.1]). The most interesting cases are K = 1, contains only one free parameter, it is related to the fold catastrophe and has two coalescing stationary points (Airy integral); K = 2, contains two free parameters, it is related to the cusp catastrophe and involves three coalescing stationary points (Pearcey integral), and K = 3, depends on three free parameters, corresponds to the swallowtail catastrophe and involves four coalescing stationary points (swallowtail integral). On the other hand, the most important canonical integrals related to umbilic catastrophes (double integrals) are the elliptic umbilic catastrophe [5,Eqs. 36 Both, the elliptic umbilic catastrophe integral and the hyperbolic umbilic catastrophe integral, depend on three free parameters and they are related to the elliptic and the hyperbolic umbilic catastrophe, respectively. These are double-sheet caustic surfaces that coalesce at a unique point corresponding to an umbilical point, an isolated point in the elliptic umbilic, and a finite-angled corner in the hyperbolic umbilic [2]. The catastrophes of internal dimension two, the so-called umbilics, have a directing role in breaking phenomena in hydrodynamics (breaking of waves and breaking of jets). In biology, they govern the morphology of engulfing phenomena like phagocytosis, neurulation, etc., and, in reproduction, the emission and spreading of gametes [26]. Mathematically speaking, the canonical integrals associated to cuspoid and umbilic catastrophes are oscillatory integrals with several nearly coincident stationary phase or saddle points. They have physical applications in the description of surface gravity waves [15,28], optics, quantum mechanics, and acoustics (see [5,Sec. 36.14] and references therein), playing a fundamental role in the uniform asymptotic approximation of oscillatory integrals. Current knowledge of these integrals, such as symmetries, bifurcation sets, zeros, convergent series expansions, differential equations, or leadingorder asymptotic approximations, among others, can be found in Chapter 36 of the Digital Library of Mathematical Functions [5]. For the particular case of the hyperbolic umbilic catastrophe, the bifurcation set of codimension three is given by (see [5,Eq. 36.4.12]) where the + sign labels the cusped sheet and the − sign labels the sheet that is smooth for z = 0 (see Fig. [5, Fig. 36. 4.4]).
An important practical problem in dealing with these integrals is represented by their numerical evaluation. Convergent expansions, methods based on the numerical integration of certain differential equations, or complex contour quadrature techniques are some of the techniques used so far (see [5,Sect. 36.15]). However, these techniques are not useful, in general, for large values of the variables. For example, convergent expansions may not include large values of the variables in the region of convergence or the convergence of the series can be very slow, and complex contour quadrature techniques are not usually well adapted to highly oscillatory integrals. Then, a more efficient approach, such as the use of asymptotic expansions, is recommended, not only for numerical computations but also from an analytic point of view. In recent papers, complete asymptotic expansions of the Pearcey [17,18] and swallowtail integrals [12][13][14] have been derived. In this paper, we are concerned with the asymptotic analysis of the hyperbolic umbilic catastrophe integral (H) (x, y, z) given in (1).
The hyperbolic umbilic catastrophe integral has important applications in semiclassical collision theory and in optics. More specifically, the hyperbolic umbilic describes diffraction by a thin water-droplet lens with a circular boundary, hanging vertically [3,Eq. 19], [7,20]. It is used in the scattering theory: the two-dimensional problem of atom-rigid rotor rotationally inelastic scattering [16] and atom-crystal scattering [1]. The hyperbolic umbilic also represents the optical pattern in the neighborhood of the focus in the limit of vanishing wavelength [21, Eq. 2.1]. For more applications and detailed information we refer to [3,4,8,9,21,29] and references therein.
In [5, Eq. 36.11.8], we can find the leading-order asymptotic of (1) along a symmetry line Berry and Howls calculated and displayed in [3] the Stokes surface for the diffraction patterns of the hyperbolic umbilic singularities for x, y, and z real. The surface in the full x, y, z space has two branches, smooth apart from finite-angled creases on the complex whiskers and a cusped edge at the cusp line of the bifurcation set. In [29], the authors propose two methods for the computation of the hyperbolic umbilic catastrophe integral: the integration of various differential equations satisfied by this function and a quadrature formula which evaluates the integral directly. This second method is based on a similar idea applied in [6] for the computation of the elliptic umbilic catastrophe integral. Borghi proposes in [8] a computational approach based on expanding a simple integral representation of (H) (x, y, z) as a convergent power series and the application of the Weniger transformation to avoid the slow convergence of the series for non-small values of the control parameters x, y, and z. In spite of the acceleration, the approximation fails for large values of the control parameters. Then, in this paper, we analyze the asymptotic behavior of (1) for |x| large and bounded values of |y|, |z|. It is worth noting that, because of the symmetry (H) (x, y, z) = (H) (y, x, z), the approximation for large |x| and bounded |y|, |z| solves the problem for large |y| and bounded |x|, |z|.
The standard saddle point method is difficult to apply to the families of catastrophe integrals because of the complexity of the phase function (see Sections 2.2 and 2.4 for further details). The modified saddle point method [19] simplifies the computations considerably and permits a straightforward derivation of an asymptotic expansion. It has been successfully applied to derive asymptotic expansions of some integrals associated to cuspoid catastrophes [12-14, 17, 18]. The hyperbolic umbilic catastrophe integral is defined by a double integral (1), but it can be transformed into a one-dimensional integral. Then, in this paper, we consider appropriate simple integral representations, and we apply the modified saddle point method [19] to obtain a complete asymptotic expansion of (H) (x, y, z).
In the following section, we analyze the asymptotic behavior of (H) (x, y, z) when |x| is large and |y|, |z| are bounded. We provide a complete and detailed application of the method [19]. By symmetry arguments, the approximation is also valid for large values of |y| and bounded values of |x| and |z|. In Section 3, some numerical experiments show the accuracy and the asymptotic character of the approximation. Throughout all the paper we use the principal argument arg w ∈ (−π, π] for any complex number w and square roots are assumed to take their principal value.

An appropriate integral representation and saddle points
We consider the integral representation of the hyperbolic umbilic catastrophe (1) given in [5,Eq. 36 where the path of integration L := (∞ exp (5πi/6) , 0]∪[0, ∞ exp (πi/6)) is depicted in Fig. 1, and Ai(w) is the Airy function [22]. In order to apply the modified saddle point method [19], we consider the change of variableu = s |x|  (2) is the union, at the origin t = 0, of two half straight lines that form an angle ∓π/6 with the negative and the positive real axes, respectively where we have defined the phase function f (s) := i(s 3 + 3e iθ s) and θ := arg x. This phase function has two saddle points: s ± := ±ie i θ 2 . We know that the asymptotically relevant saddle points are those ones (only s + , only s − , or both) for which the integration path L can be deformed into a steepest descent path (or union of steepest descent paths) that contains these points.

Steepest descent paths
On the one hand, the standard saddle point method requires the analytic expression of the steepest descent paths of f (s) at the saddle points, expression that is not straightforward to obtain. On the other hand, we know from [19] that the asymptotic analysis of the integral (2) does not require the computation of the steepest descent paths of the complete phase function f (s) at the saddle points s ± , but the steepest descent paths of the main part of f (s) at s ± , that may always be computed in a systematic and straightforward manner, as they are nothing but straight lines [19].
At the two saddle points s ± , the first non-vanishing derivative of f (s) is f (s). Then, denoting by φ ± the phase of f (s ± ), and by f ± (s) the Taylor polynomial of degree 2 of f (s) at the saddle point Then, the phase function may be rewritten in the form At every saddle point s ± , the main part of f (s) is just the Taylor polynomial of degree 2 of f (s) at s ± , f ± (s). Then, at each saddle point s ± , we have that the steepest descent paths of f ± (s) are the following straight lines through the saddle points [19]:

Deformation of the integration path
Following [19], the next step consists in the deformation of the integration path L to a new path appropriate for the asymptotic analysis: the path must contain certain portions ± of some of the two steepest descent straight lines ± . By Cauchy's residue theorem we can deform L → L 2 ∪ ∪ L 1 , where L 1 and L 2 are residual portions of the original path L, and is a portion of one or two appropriate steepest descent paths ± . The precise form of the deformation of the original path L depends on θ (−π < θ ≤ π ) and is given by

Approximate computation of (6)
The next point in our analysis is the approximate computation of the right-hand side of (6) when , L 1 , and L 2 are the paths described previously. The application of the standard saddle point method at this step would require a change of variable difficult to implement that it is not necessary with the modified technique. Then, in the application of our method, we must have into account the following observations: • ( f (s)) is a decreasing function in both paths, L 1 and L 2 , as |s| → +∞.
Then, from the above observations, we conclude that

Computation of the integrals (9)
The last step in our analysis is the computation of the integrals (9) over the paths ± . For that purpose, at every steepest descent path ± , we separate the phase function f (s) in the form specified in (4), as the sum of the main part f ± (s) plus a residual term of order three ± (x, y, z) ∼ u for − in (10), we find with w + := 1, w − := −i, and Consider now the Taylor expansion of g ± (x, y, z, u) at u = 0, where, for convenience, we have written the Taylor coefficients of g ± (x, y, z, u) at , y, z). It can be checked that g ± (x, y, z, u) is a solution to the second-order differential equation in the variable u: Introducing (12) in (13) and identifying the coefficients of each power u n , n = 0, 1, . . ., we obtain that the sequence {b ± n (x, y, z)} n=0,1,2,... is the solution to the sixth-order recurrence relation with initial values b ± 0 (x, y, z) = Ai Introducing the Taylor expansion (12) in (11), we find

The resulting asymptotic expansion
Therefore, from (8) we finally derive the asymptotic expansion for large |x| of the hyperbolic umbilic catastrophe integral in the form with b ± 2n (x, y, z) given in (14). In particular, the leading-order approximation is given by Observation. The expansion (15) is not a genuine Poincaré expansion, since the coefficients b ± 2n (x, y, z) depend on the asymptotic variable x. More precisely, from the recurrence relation (14) and the first five coefficients b n , n = 0, 1, . . . , 5, it can be shown that b n = O(x n/2 /2 )Ai a √ x + c for n even, and b n = O(x n/2 /2 )Ai a √ x + c for odd n, with a := ±3 −5/6 i z, c := 3 −1/3 y. Then, although the asymptotic sequence b ± 2n (x, y, z)/x n 2 is not a genuine Poincaré sequence, its asymptotic order descents in the form of a sawtooth and the asymptotic series (15) provides a good asymptotic approximation at any truncation of the number of terms.  (15) and their asymptotic character. The Table 1 Relative errors in the approximation of (H ) (x, y, z) given in (2) by using (15) with the series truncated after n terms y = z = 1 y = 0.5e i π 9 , z = 0.25 computations have been carried out by using the symbolic manipulation program Wolfram Mathematica 12.2. In particular, the "exac" value of (H ) (x, y, z) is computed by means of numerical integration with the command "NIntegrate." We remark at this point that a general implementation of an algorithm of the modified saddle point method (mixing symbolic and numeric computations) is not possible, since it strongly depends on the phase function of the integrand. The computation of the saddle points and the steepest descent paths is quite handicraft and requires a rigorous previous study.

Numerical experiments
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.