On the Approximate Evaluation of Some Oscillatory Integrals

Abstract

To determine the photon emission or absorption probability for a diatomic system in the context of the semiclassical approximation it is necessary to calculate the characteristic canonical oscillatory integral which has one or more saddle points. Integrals like that appear in a whole range of physical problems, e.g., the atom–atom and atom–surface scattering and various optical phenomena. A uniform approximation of the integral, based on the stationary phase method is proposed, where the integral with several saddle points is replaced by a sum of integrals each having only one or at most two real saddle points and is easily soluble. In this way we formally reduce the codimension in canonical integrals of “elementary catastrophes” with codimensions greater than 1. The validity of the proposed method was tested on examples of integrals with three saddle points (“cusp” catastrophe) and four saddle points (“swallow-tail” catastrophe).

1. Introduction

Many physical processes like the atom scattering on atom, molecule, surface, or a crystal and propagation of electromagnetic waves through a medium, are described generally by multidimensional integrals. The value of these integrals as a function of the parameters $\mathbf{a}=\left\{a_1,a_2,\dots a_n\right\}$ depends on the morphology of the phase $f(\mathbf{a};\mathbf{u})$ of the integrand, where $\mathbf{u}=\left\{u_1,u_2,\dots u_k\right\}$ are the variables of integration.

The main contribution to the integral comes from the neighbourhood of the saddle points ${\mathbf{u}}_s=\left\{u_{1s},u_{2s},\dots u_{ks}\right\}$, where $\frac{\partial f(\mathbf{a};{\mathbf{u}}_{\mathbf{s}})}{\partial u_{1s}}=\frac{\partial f(\mathbf{a};{\mathbf{u}}_{\mathbf{s}})}{\partial u_{2s}}=\dots =\frac{\partial f(\mathbf{a};{\mathbf{u}}_{\mathbf{s}})}{\partial u_{ks}}=0$. For some parameters a it is possible that some higher derivatives are equal to 0 at the saddle points, which causes coalescence of these points, increases the value of the integral that manifests itself as “rainbow and glory” [1]. To analyze the rainbow phenomenon, it is important to know the caustic surface $C({\mathbf{a}}_{\mathbf{c}})$, defined by the relation $\frac{\partial f({\mathbf{a}}_c;u_{\mathbf{s}})}{\partial u_{1s}}=\frac{{\partial }^2\partial f({\mathbf{a}}_c;u_{\mathbf{s}})}{\partial u_{1s}{}^2}=\dots =\frac{\partial \partial f({\mathbf{a}}_c;u_{\mathbf{s}})}{\partial u_{ks}}=\frac{{\partial }^2\partial f({\mathbf{a}}_c;u_{\mathbf{s}})}{\partial u_{ks}{}^2}=0$.

By transforming the variables using a uniform one-to-one mapping, the phases can be transformed into simple forms, which are classified as seven “elementary catastrophes” [2,3]. There are four one-dimensional catastrophes: fold $f(x;u)=u^3+xu$, cusp $f(x,y;u)=u^4+x{u}^2+yu$, swallow-tail $f(x,y,z;u)=u^5+x{u}^3+y{u}^2+zu$, butterfly $f(x,y,z,w;u)=u^6+x{u}^4+y{u}^3+z{u}^2+wu$, and three two dimensional ones: hyperbolic umbilic $f(x,y,z;u,v)=u^3+v^3+xuv-yu-zv$, elliptic umbilic $f(x,y,z;u,v)=u^3-3u{v}^2+x\left(u^2+v^2\right)-yu-zv$, parabolic umbilic $f(x,y,z,w;u,v)=u^{2}v+v^4+x{u}^2+w{v}^2-yu-zv$. The elementary Thom’s catastrophes in the context of the theory of atomic collisions have been discussed in detail by Connor [4].

In the cuspoid case (one integration variable), the oscillatory integrals are usually written in the form:

$$I(\mathbf{a})=\underset{-\infty }{\overset{\infty }{\int }}g(u)e^{if(\mathbf{a};u)}du$$

(1)

where a = {a₁, a₂, …} is a set of parameters. As a varies, as many as K + 1 (real or complex) critical points of the smooth, real-valued phase function f can coalesce in clusters of two or more. The function g has a smooth amplitude. In what follows we denote $\frac{{\partial }^n}{\partial u^n}f(\mathbf{a};u)=f^{(n)}(\mathbf{a};u)$. The critical (stationary) points u_j(a), 1 ≤ j ≤ K + 1, are defined by $f^{(1)}(\mathbf{a};u_j)=0$ [5].

In the case of a single real critical point the integral $I(\mathbf{a})$ is in the leading order approximated by [6]:

$$I(\mathbf{a})\approx I_q(\mathbf{a};u_1)=\sqrt{\frac{2\pi i}{f^{(2)}(\mathbf{a};u_1)}}g(u_1)e^{if(\mathbf{a};u_1)}=\sqrt{\frac{2\pi }{\left|f^{(2)}(\mathbf{a};u_1)\right|}}g(u_1)e^{if(\mathbf{a};u_1)+i{\sigma }_1\frac{\pi }{4}}$$

(2)

where ${\sigma }_1=\mathrm{sgn}(f^{(2)}(\mathbf{a};u_1))$ and the subscript q indicates a quadratic expansion of $f(\mathbf{a};u)$ around $u_1$. The result is easily generalized to the case of j_max ($1\le j_{\max }\le K$) isolated real critical points [7]. The main contribution to the integral comes from the regions around the stationary points $u_j$ where the phase function $f(\mathbf{a};u)$ is slowly varying.

Since the positions of the critical points depend on a, they can move close together and coalesce as a varies. In the uniform asymptotic evaluation of oscillatory integrals the result is expressed in terms of certain canonical integrals [5,7] and their derivatives. Each canonical integral is characterized by a given number of coalescing critical points. One defines a mapping u(a;t) by relating f(a;u) to the normal form of cuspoid catastrophes ${\Phi }_K(\mathbf{b};t)$ in the following way:

$$f(\mathbf{a};u(\mathbf{a};t))=A(\mathbf{a})+{\Phi }_K(\mathbf{b}(\mathbf{a});t)=A(\mathbf{a})+t^{K+2}+{\displaystyle \sum _{m=1}^K{b}_m{t}^m}$$

(3)

with the K + 1 functions A(a) and b(a) determined by the correspondence of K + 1 critical points of f and Φ_K.

In the simplest case of two coalescing critical points (K = 1, fold catastrophe) there is a single point $u_e=u({\mathbf{a}}_e)$ where $f^{(2)}(\mathbf{a};u_e)=0$, i.e., the function $f^{(1)}(\mathbf{a},u)$ has an extremum and there are two stationary points $u_1(\mathbf{a})$ and $u_2(\mathbf{a})$. In some range of the parameter a the stationary points are real and $u_1\le u_e\le u_2$. For $\mathbf{a}={\mathbf{a}}_e$ the two points coalesce and $u_1=u_e=u_2$. For other values of a the stationary points are complex conjugate solutions of Equation (2), i.e., $u_1=u_2^{*}$.

The leading-order uniform approximation in the case of the fold catastrophe is given by [8]

$$I(\mathbf{a})\approx I_F(\mathbf{a})=\pi \sqrt{2}{e}^{iA(\mathbf{a})}\left[\left(\frac{g(u_1)}{\sqrt{{\sigma }_1{f}^{(2)}(\mathbf{a};u_1)}}+\frac{g(u_2)}{\sqrt{-{\sigma }_1{f}^{(2)}(\mathbf{a};u_2))}}\right)\sqrt[4]{\zeta (\mathbf{a})}Ai(-\zeta (\mathbf{a}))-i{\sigma }_1\left(\frac{g(u_1)}{\sqrt{{\sigma }_1{f}^{(2)}(\mathbf{a};u_1)}}-\frac{g(u_2)}{\sqrt{-{\sigma }_1{f}^{(2)}(\mathbf{a};u_2))}}\right)\frac{A{i}^{\prime }(-\zeta (\mathbf{a}))}{\sqrt[4]{\zeta (\mathbf{a})}}\right]$$

(4)

where $A(\mathbf{a})=\frac{1}{2}\left[f(\mathbf{a};u_2)+f(\mathbf{a};u_1)\right]$, $\zeta {(\mathbf{a})}^{3/2}={\sigma }_1\frac{3}{4}\left[f(\mathbf{a};u_2)-f(\mathbf{a};u_1)\right]$ and ${\sigma }_i=\mathrm{sgn}\left(f^{(2)}(\mathbf{a};u_i)\right)$. The branches are chosen so that $\zeta (\mathbf{a})$ is real and positive if the critical points are real, or real and negative if they are complex. (Note that ${\sigma }_i{f}^{(2)}(\mathbf{a};u_i)=\left|f^{(2)}(\mathbf{a};u_i)\right|$ and in the case under consideration ${\sigma }_2=-{\sigma }_1$).

The transitional approximation $I_F^{tr}(\mathbf{a})$ reproduces the uniform approximation $I_F(\mathbf{a})$ on the neighbourhood of $\mathbf{a}\approx {\mathbf{a}}_e$ ($I_F^{tr}({\mathbf{a}}_e)=I_F({\mathbf{a}}_e)$) and enables analytical continuation from the region of real stationary points into the region of complex ones. The transitional approximation is given by [9],

$$I_F^{tr}(\mathbf{a})=2\pi e^{iA(\mathbf{a})}g\left(u_e\right){\left(\frac{2}{\left|f^{(3)}(\mathbf{a};u_e)\right|}\right)}^{1/3}\left\{Ai(-\zeta (\mathbf{a}))+i{\sigma }_1{\left(\frac{2}{\left|f^{(3)}(\mathbf{a};u_e)\right|}\right)}^{1/3}\left[\frac{g^{(1)}(u_e)}{g\left(u_e\right)}-\frac{1}{6}\frac{f^{(4)}(\mathbf{a};u_e)}{f^{(3)}(\mathbf{a};u_e)}\right]A{i}^{\prime }(-\zeta (\mathbf{a}))\right\}$$

(5)

where $\zeta (\mathbf{a})={\sigma }_1{\left(\frac{2}{\left|f^{(3)}(\mathbf{a};u_e)\right|}\right)}^{1/3}{f}^{(1)}(\mathbf{a};u_e)$ and $A(\mathbf{a})=f(\mathbf{a};u_e)+\frac{1}{6}{\left(\frac{f^{(1)}(\mathbf{a};u_e)}{f^{(3)}(\mathbf{a};u_e)}\right)}^2{f}^{(4)}(\mathbf{a};u_e)$.

In order to obtain A(a) and b(a), a set of nonlinear equations has to be solved. These can be solved in principle, but there are practical difficulties in attempting a solution [10]. On the other hand, away from b = 0 the canonical integrals can be approximated in terms of canonical integrals ${\Phi }_J$ corresponding to lower-order catastrophes (i.e., J < K) [7,10,11,12,13,14,15,16].

The motivation to study these types of integrals originates from the investigation of optical spectra of diatomic molecules [9]. For example, in the semiclassical approximation the matrix element of the dipole moment $D(R)$ for the optical transition is proportional to the integral [11]

$$\underset{-\infty }{\overset{\infty }{\int }}dtD(R(t))\exp \left(\frac{i}{ħ}\underset{0}{\overset{t}{\int }}\left(\Delta (R(t^{\prime }))-h\nu \right)d{t}^{\prime }\right).$$

(6)

The radial movement of atoms is described classically, $R=R(t)$. The phase function in the integral (6) is $f(t)=\frac{1}{ħ}\underset{0}{\overset{t}{\int }}\left(\Delta (R(t^{\prime }))-h\nu \right)d{t}^{\prime }$, where $\Delta (R)$ is the energy difference of the upper and lower electronic state energies. The condition $f^{(1)}(t)=0$ gives the saddle points which satisfy the classical Franck–Condon condition $\Delta (R(t_c))=h\nu$. If there are points $t_t$ satisfying the condition $f^{(3)}(t_t)=0$, the method suggested in Section 2 of this paper is a good choice to calculate the integral in Equation (6).

In the following sections we propose a new procedure for the approximate evaluation of oscillatory integrals with several stationary points.

2. A New Procedure for Approximate Evaluation of Oscillatory Integrals

Let there be a point $u_t$ in the integration interval $\left[-\infty ,\infty \right]$ which satisfies the condition $f^{(3)}(\mathbf{a};u_t)=0$. In the neighbourhood of this point one defines a function:

$$F(\mathbf{a};u)=f(\mathbf{a};u_t)+f^{(1)}(\mathbf{a};u_t)\left(u-u_t\right)+\frac{1}{2}{f}^{(2)}(\mathbf{a};u_t){\left(u-u_t\right)}^2+\frac{1}{4!}{f}^{(4)}(\mathbf{a};u_t){\left(u-u_t\right)}^4$$

(7)

The first derivative of this function, $F^{(1)}(\mathbf{a};u)=f^{(1)}(\mathbf{a};u_t)+f^{(2)}(\mathbf{a};u_t)\left(u-u_t\right)+\frac{1}{3!}{f}^{(4)}(\mathbf{a};u_t){\left(u-u_t\right)}^3$, has an inflection at the point $u_t$. If $\mathrm{sgn}\left(f^{(2)}(\mathbf{a};u_t)\right)=\mathrm{sgn}\left(f^{(4)}(\mathbf{a};u_t)\right)$, $F^{(1)}(\mathbf{a};u)$ is monotonic function. In the case when $\mathrm{sgn}\left(f^{(2)}(\mathbf{a};u_t)\right)=-\mathrm{sgn}\left(f^{(4)}(\mathbf{a};u_t)\right)$, the function $F^{(1)}(\mathbf{a};u)$ has two extremes at real points $u_{1,2}=u_t\mp \sqrt{-\frac{2f^{(2)}(\mathbf{a};u_t)}{f^{(4)}(\mathbf{a};u_t)}}$.

If there are m points $u_{p,i}\in \left[-\infty ,\infty \right]$, $i=1,\dots ,m$, satisfying $f^{(3)}(\mathbf{a};u_{p,i})=0$ and ${\sigma }_{p,i}=\mathrm{sgn}\left(f^{(2)}(\mathbf{a};u_{p,i})\right)=-\mathrm{sgn}\left(f^{(4)}(\mathbf{a};u_{p,i})\right)$, these points divide the interval $\left[-\infty ,\infty \right]$ into m + 1 intervals $\left[u_{p,i-1},u_{p,i}\right]$ and the integral $I(a)$ can be written:

$$I(\mathbf{a})=\underset{-\infty }{\overset{\infty }{\int }}g(u)e^{if(\mathbf{a};u)}du=\sum _{i=1}^{m+1}\underset{u_{p,i-1}}{\overset{u_{p,i}}{\int }}g(u)e^{if(\mathbf{a};u)}du$$

(8)

where the end points of the integration $u_{p,0}=-\infty$ and $u_{p,m+1}=\infty$ (${\sigma }_{p,0}=\mathrm{sgn}\left({\lim }_{u\to -\infty }{f}^{(2)}(\mathbf{a};u)\right)$ and ${\sigma }_{p,m+1}=\mathrm{sgn}\left({\lim }_{u\to \infty }{f}^{(2)}(\mathbf{a};u)\right)$ have been introduced. At each interval $\left[u_{p,i-1},u_{p,i}\right]$ the function $f^{(1)}(\mathbf{a};u)$ has a simple property. If ${\sigma }_{p,i-1}={\sigma }_{p,i}$, the function $f^{(1)}(\mathbf{a};u)$ is monotonic on the interval $\left[u_{p,i-1},u_{p,i}\right]$ and has a single real saddle point. In the case ${\sigma }_{p,i-1}=-{\sigma }_{p,i}$ there is a point $u_e\in \left[u_{p,i-1},u_{p,i}\right]$, $f^{(2)}(\mathbf{a};u_e)=0$ and the function $f^{(1)}(\mathbf{a};u)$ has an extreme at $u_e$ and two saddle points.

One defines a function $f_{p,i}(\mathbf{a};u)$ as a series expansion of the phase $f(\mathbf{a};u)$ around $u_{p,i}$ up to the quadratic term: $f_{p,i}(\mathbf{a};u)=f(\mathbf{a};u_{p,i})+f^{(1)}(\mathbf{a};u_{p,i})(u-u_{p,i})+\frac{1}{2!}{f}^{(2)}(\mathbf{a};u_{p,i}){(u-u_{p,i})}^2$. Note that $f_{p,i}(\mathbf{a};u_{p,i})=f(\mathbf{a};u_{p,i})$, $f_{p,i}^{(1)}(\mathbf{a};u_{p,i})=f^{(1)}(\mathbf{a};u_{p,i})$, $f_{p,i}^{(2)}(\mathbf{a};u_{p,i})=f^{(2)}(\mathbf{a};u_{p,i})$, and $f_{p,i}^{(3)}(\mathbf{a};u_{p,i})=f^{(3)}(\mathbf{a};u_{p,i})=0$.

We define the integral $I_{p,i}(\mathbf{a})=\underset{-\infty }{\overset{\infty }{\int }}g(u_{p,i})e^{i{f}_{p,i}(\mathbf{a};u)}du$, which has an exact solution

$$I_{p,i}^{}(\mathbf{a})=\sqrt{\frac{2\pi i}{f^{(2)}(\mathbf{a};u_{p,i})}}g(u_{p,i})e^{i\left(f(\mathbf{a};u_{p,i})-\frac{1}{2}\frac{f^{(1)}{(\mathbf{a};u_{p,i})}^2}{f^{(2)}(\mathbf{a};u_{p,i})}\right)}$$

(9)

Relation (8) can be written as:

$$\begin{array}{ll}I(\mathbf{a})& =\sum _{i=1}^{m+1}\underset{u_{p,i-1}}{\overset{u_{p,i}}{\int }}g(u)e^{if(\mathbf{a};u)}du+\sum _{i=1}^m{I}_{p,i}(\mathbf{a})-\sum _{i=1}^m{I}_{p,i}(\mathbf{a})\\ & =\sum _{i=1}^{m+1}\underset{u_{p,i-1}}{\overset{u_{p,i}}{\int }}g(u)e^{if(\mathbf{a};u)}du+\sum _{i=1}^m\underset{-\infty }{\overset{u_{p,i}}{\int }}g(u_{p,i})e^{i{f}_{p,i}(\mathbf{a};u)}du+\sum _{i=1}^m\underset{u_{p,i}}{\overset{\infty }{\int }}g(u_{p,i})e^{i{f}_{p,i}(\mathbf{a};u)}du-\sum _{i=1}^m{I}_{p,i}(\mathbf{a})\end{array}$$

By combining integrals in the first three sums a simple expression is obtained:

$$I(\mathbf{a})=\sum _{i=1}^{m+1}{I}_i(\mathbf{a})-\sum _{i=1}^m{I}_{p,i}(\mathbf{a})$$

(10)

where integrals $I_i(\mathbf{a})$ are

$$I_i(\mathbf{a})=\underset{-\infty }{\overset{\infty }{\int }}{g}_i(u)e^{i{f}_i(\mathbf{a};u)}du$$

(11)

The functions $g_i(u)$ and $f_i(\mathbf{a};u)$ are shown in

Table 1. Functions $g_i(u)$ and $f_i(\mathbf{a};u)$ for i = 1, i = m+1 and 1 < i < m + 1.

$\begin{array}{ccc}i=1& -\infty \le u<u_{p,1}& u\ge u_{p,1}\\ f_1(\mathbf{a};u)=& f(\mathbf{a};u)& f_{p,1}(\mathbf{a};u)\\ g_1(u)=& g(u)& g(u_{p,1})\end{array}$	$\begin{array}{ccc}i=m+1& u<u_{p,m}& u_{p,m}\le u\le \infty \\ f_{m+1}(\mathbf{a};u)=& f_{p,m}(\mathbf{a};u)& f(\mathbf{a};u)\\ g_{m+1}(u)=& g(u_{p,m})& g(u)\end{array}$
$\begin{array}{cccc}1<i<m+1& u<u_{p,i-1}& u_{p,i-1}\le u<u_{p,i}& u\ge u_{p,i}\\ f_i(\mathbf{a};u)=& f_{p,i-1}(\mathbf{a};u)& f(\mathbf{a};u)& f_{p,i}(\mathbf{a};u)\\ g_i(u)=& g(u_{p,i-1})& g(u)& g(u_{p,i})\end{array}$

.

The relation (10) can be generalized to be valid for the case $m=0$ as well:

$$I(\mathbf{a})={\delta }_{m,0}I(\mathbf{a})+\left(1-{\delta }_{m,0}\right)\left(\sum _{i=1}^{m+1}{I}_i(\mathbf{a})-\sum _{i=1}^m{I}_{p,i}(\mathbf{a})\right)$$

(12)

So far, no approximation has been made. The relation (10) is an identity. An integral of the function, the phase of which has several real stationary points, is divided into the sum of the integrals $I_i(\mathbf{a})$ whose phase functions $f_i(\mathbf{a};u)$ have either one or at most two stationary points.

If the phase function $f_i(\mathbf{a};u)$ has only one real saddle point and its first derivative $f_i^{(1)}(\mathbf{a};u)$ is monotonic, the value of integral $I_i(\mathbf{a})$ can be calculated using Equation (2). In the case when the function $f_i^{(1)}(\mathbf{a};u)$ has a single extreme at the point $u_e$ and two real or complex saddle points, the integral is easily soluble using the approximate methods described in the introduction (Equation (4)). If the phase $f_i(\mathbf{a};u)$ is given by numerical points in the region where a complex pair of saddle points contributes to the integral, the analytical continuation of (5) can be used. The numerical accuracy of this method is determined by the accuracy of the leading-order uniform approximations (2) and (4).

3. Results

The method outlined in Section 2 was tested on three examples that are typical for the spectra of diatomic molecules. For simplicity we use the phase function given by the polynomial phase of the Thom’s elementary catastrophe. The case when $f^{(1)}(t)$ and the difference potential $\Delta (R)$ are both monotonic functions with a single inflection point is illustrated by the analysis of the Pearcey integral $P(x\ge 0,y)$ in Section 3.1.1. In Section 3.1.2 with the Pearcey integral $P(x<0,y)$ we analyze the case when the function $f^{(1)}(t)$ and the difference potential $\Delta (R)$ have two extremes and one inflection point. Finally, in Section 3.2 we illustrate the case when the difference potential has an extreme near the turning point by the analyses of the swallow-tail catastrophe integral $S(x<0,0,z)$. For simplicity, we take $g(u)=1$ in all the examples. The dependence of the integral (6) on the variable transition dipole moment was discussed by Beuc et al. in [9].

3.1. Cusp Catastrophe (K = 2)

Let us consider the Pearcey integral, the canonical integral for the cusp catastrophe ($x,y\in R$):

$$P(x,y)=\underset{-\infty }{\overset{\infty }{\int }}{e}^{i\left(u^4+x{u}^2+yu\right)}du$$

(13)

(Other notations also appear in the literature). The Pearcey integral is symmetrical with respect to variable y: $P(x,y)=P(x,-y)$. For the numerical integration of the Pearcey integral we used the form:

$P(x,y)=2e^{i\frac{\pi }{8}}\underset{0}{\overset{\infty }{\int }}{e}^{-t^4+e^{i\frac{3\pi }{4}}x{t}^2}\cos (e^{i\frac{5\pi }{8}}yt)dt$ [13]. The numerical evaluation of the integral and all other calculations in this paper were done using the Wolfram Mathematica 11.3 computing system.

The phase function in (6) is $f(x,y;u)=u^4+x{u}^2+yu$. There are three saddle points defined by the condition $f^{(1)}(x,y;u)=4u^3+2xu+y=0$ (Figure 1 and Figure 3a):

$$\begin{array}{l}{u}_1(x,y)=\frac{1+i\sqrt{3}}{12}\frac{x}{U(x,y)}-\frac{1-i\sqrt{3}}{2}U(x,y)\\ u_2(x,y)=\frac{1-i\sqrt{3}}{12}\frac{x}{U(x,y)}-\frac{1+i\sqrt{3}}{2}U(x,y)\\ u_3(x,y)=-\frac{1}{6}\frac{x}{U(x,y)}+U(x,y)\end{array}$$

(14)

where $U(x,y)=\frac{1}{2}\sqrt[3]{\sqrt{\delta (x,Y)}-y}$ and $\delta =\frac{8}{27}{x}^3+y^2$. If δ < 0, all saddle points are real and if δ > 0, one saddle point is real and the other two are complex conjugates of each other (Figure 1 and Figure 3a). There are two bifurcation points $u_{e1,e2}=\mp \sqrt{-\frac{1}{6}x}$ defined by the relation $f^{(2)}(x,y;u)=12u^2+2x=0$. There is a single point $u_p=0$ where $f^{(3)}(x,y;u_p)=0$ and $\frac{f^{(2)}(x,y;u_p)}{f^{(4)}(x,y;u_p)}=\frac{x}{12}$, i.e., $\mathrm{sgn}\left(\frac{f^{(2)}(x,y;u_p)}{f^{(4)}(x,y;u_p)}\right)=\mathrm{sgn}(x)$. In the special case $x=0$ one has $u_{e1}=u_{e2}=u_p=0$.

3.1.1. Case $x>0$

If $x\ge 0$, then δ is always positive, the saddle point $u_3(x,y)$ is real and the points $u_1(x,y)$ and $u_2(x,y)$ are complex conjugates (Figure 1).

The function $f^{(1)}(1,y;u)$ is monotonic and, according to Equation (2), the value of the Pearcey integral can be approximated as

$$P(x\ge 0,y)\approx P_q(x\ge 0,y)=\sqrt{\frac{\pi i}{6u_3^2+x}}{e}^{i(u_3^4+x{u}_3^2+y{u}_3)}$$

(15)

From Figure 2 and

Table 2. Comparison values of $P(x,y)$ and approximate function $P_q(x,y)$ for $y=0$ and $x=0$.

x	P(x,0)	Pq(x,0)	y	P(0,y)	Pq(0,y)
1	1.20838 + 0.779288i	1.25331 + 1.25331i	1	1.55093 + 0.427893i	1.09286 + 0.353605i
2	0.924029 + 0.729007i	0.886227 + 0.886227i	2	1.12475 − 0.17608i	0.837873 − 0.359347i
3	0.754294 + 0.657362i	0.723601 + 0.723601i	3	0.384485 − 0.642953i	0.244423 − 0.757992i
4	0.646979 + 0.593695i	0.626657 + 0.626657i	4	−0.385925− 0.545144i	−0.434337 − 0.578749i
5	0.573931 + 0.541858i	0.560499 + 0.560499i	5	−0.670195 + 0.0710806i	−0.66748 + 0.0754601i
6	0.520848 + 0.500054i	0.511663 + 0.511663i	6	−0.235367 + 0.592027i	−0.214718 + 0.594539i
7	0.480234 + 0.465936i	0.473708 + 0.473708i	7	0.430078 + 0.415615i	0.442629 + 0.405754i
8	0.447916 + 0.437618i	0.443113 + 0.443113i	8	0.510179 − 0.260969i	0.506489 − 0.270769i
9	0.421413 + 0.413718i	0.417771 + 0.417771i	9	−0.1039 − 0.542069i	−0.112749 − 0.540578i
10	0.399166 + 0.393242i	0.396333 + 0.396333i	10	−0.532222 − 0.0242518i	−0.532896 − 0.0165838i

, we can freely estimate that the difference of the approximation $P_q(x,y)$ and the exact values of $P(x,y)$ is smaller than few percent if the condition $\sqrt{x^2+y^2}>5$ is satisfied.

3.1.2. Case $x<0$

According to Section 2, when $x<0$ using the relation (10) the Pearcy integral can be written as:

$$\begin{array}{ll}P(x,y)& =I_1(x,y)+I_2(x,y)-I_{p,1}(x,y)\\ & =\underset{-\infty }{\overset{\infty }{\int }}{e}^{i{f}_1(x,y;u)}du+\underset{-\infty }{\overset{\infty }{\int }}{e}^{i{f}_2(x,y;u)}du-\underset{-\infty }{\overset{\infty }{\int }}{e}^{i{f}_{p,1}(x,y;u)}du\end{array}$$

(16)

Here the phase functions have the form $f_{p,1}(x,y;u)=x{u}^2+yu$, $f_1(x,y;u)=\{\begin{array}{cc}f(x,y;u)& u<0\\ f_p(x,y;u)& u\ge 0\end{array}$, $f_2(x,y;u)=\{\begin{array}{cc}{f}_p(x,y;u)& u<0\\ f(x,y;u)& u\ge 0\end{array}$. It is easy to show that $f_2(x,y,u)=f_1(x,-y,-u)$, $I_2(x,y)=I_1(x,-y)$, and the Pearcey integral can be decomposed exactly as

$$P(x,y)=I_1(x,y)+I_1(x,-y)-\sqrt{\frac{\pi }{\left|x\right|}}{e}^{-i\frac{y^2}{4x}+i\frac{\pi }{4}}$$

(17)

The function $f_1(x,y;u)$ has on the interval $u\in \left[-\infty ,\infty \right]$ only one bifurcation point $u_{e1}=-\sqrt{\frac{\left|x\right|}{6}}$, where $f_1^{(2)}(x,y;u_{e1})=0$, and two saddle points ${\tilde{u}}_1(x,y)=u_1(x,y)$, ${\tilde{u}}_2(x,y)=\{\begin{array}{cc}{u}_2(x,y)& y\le 0\\ \frac{y}{2\left|x\right|}& y>0\end{array}$ (Figure 3b).

Using relation (4), the integral $I_1(x,y)$ can be approximated as:

$$I_1(x,y)\approx I_F(x,y)=\pi e^{iA(x,y)}\left[\left(\frac{1}{\sqrt{6u_1^2+x}}+\frac{1}{\sqrt{-6u_2^2\mathsf{\Theta }(-y)-x}}\right)\sqrt[4]{\zeta (x,y)}Ai(-\zeta (x,y))-i\left(\frac{1}{\sqrt{6u_1^2+x}}-\frac{1}{\sqrt{-6u_2^2\mathsf{\Theta }(-y)-x}}\right)\frac{A{i}^{\prime }(-\zeta (x,y))}{\sqrt[4]{\zeta (x,y)}}\right]$$

(18)

where, $A(x,y)=\frac{1}{4}\{\begin{array}{cc}x\left(u_1^2+u_2^2\right)+\frac{3y}{2}\left(u_1+u_2\right)& y\le 0\\ x{u}_1^2+\frac{3y}{2}{u}_1-\frac{y^2}{2x}& y>0\end{array}$, $\zeta {(x,y)}^{3/2}={\sigma }_1\frac{3}{8}\{\begin{array}{cc}x\left(u_1^2-u_2^2\right)+\frac{3y}{2}\left(u_1-u_2\right)& y\le 0\\ x{u}_1^2+\frac{3y}{2}{u}_1+\frac{y^2}{2x}& y>0\end{array}$.

We define the Airy approximation of the Pearcey integral as:

$$P_A(x,y)=I_F(x,y)+I_F(x,-y)-\sqrt{\frac{\pi }{\left|x\right|}}{e}^{-i\frac{y^2}{4x}+i\frac{\pi }{4}}$$

(19)

Paris obtained the asymptotic form of P(x,y) by considering its analytic continuation to arbitrary complex variables x and y [15]. In

Table 3. Values of P(x,y) obtained for large negative values of x when y = 2 and 4 compared to the asymptotic values [15] and the present work.

x	y	P(x,y)	Asymptotic [11]	PA(x,y)
−4	2	1.96341 − 0.73419i	1.97363 − 0.72605i	1.96482 − 0.72731i
−6	2	0.96527 + 0.46413i	0.96537 + 0.46415i	0.96366 + 0.46538i
−8	2	1.00422 − 0.11480i	1.00422 − 0.11480i	1.004077 −0.11392i
−4	4	0.14360 + 0.90244i	−	0.14063 + 0.90013i
−6	4	0.29478 − 0.84373i	0.29399 − 0.84356i	0.29629 − 0.84406i
−8	4	0.75372 − 0.23933i	0.75371 − 0.23933i	0.75379 − 0.23889i

. we compare some values of P(x,y) for large negative values of x when y = 2 and 4 to the asymptotic values [15] and the present work.

Kaminski [16] rewrites (8) as a sum of two contour integrals, one of which has exactly two relevant coalescing saddle points. This allows him to apply a cubic transformation introduced by Chester, Friedman, and Ursell [17] and to construct a uniform asymptotic expansion of (7) as x → −∞ with δ varying in an interval containing 0. The leading-order approximation was already given by Connor [12] and Connor and Farrelly [13]. In

Table 4. Comparison of values of $P(x,y)$ on the caustic with Kaminski [16] and $P_A(x,y)$.

x	y	P(x,y)	Kaminski [10]	PA(x,y)
−1.0	0.544331	2.14158 + 0.0990191i	2.34415 + 0.00118008i	2.1003 + 0.156994i
−2.0	1.5396	0.962205 − 0.450083i	0.926925 − 0.428207i	0.965935 − 0.448303i
−3.0	2.82843	1.13215 + 1.19182i	1.14743 + 1.19594i	1.12358 + 1.19408i
−4.0	4.35465	−0.142478 + 1.20972i	−0.143582 + 1.2217i	−0.146125 + 1.20649i
−5.0	6.08581	−0.888104 + 0.979074i	−0.890885 + 0.983784i	−0.889185 + 0.975844i
−6.0	8.	−1.10157 + 0.582286i	−1.09951 + 0.581515i	−1.1015 + 0.58047i
−7.0	10.0812	−0.249906 − 0.91133i	−0.249866 − 0.914663i	−0.248282 − 0.910954i
−8.0	12.3168	0.321769 − 0.468203i	0.324275 − 0.466919i	0.321939 − 0.467325i
−9.0	14.6969	0.495502 + 0.309572i	0.495034 + 0.311661i	0.494746 + 0.309898i
−10.0	17.2133	−0.704129 + 0.779039i	−0.704954 + 0.779772i	−0.70467 + 0.778148i

the values of $P(x,y)$ are compared to Kaminski’s results [16] and the approximation $P_A(x,y)$ at some points on the caustic $y={\left(-\frac{2}{3}x\right)}^{\frac{3}{2}}$.

From Table 3 and Table 4, and Figure 4, we estimate that the difference of the approximation $P_A(x,y)$ and the exact values of $P(x,y)$ is smaller than a few percent if the condition $\sqrt{x^2+y^2}>4$ is satisfied.

3.2. Swallow-Tail Catastrophe (K = 3)

The swallow-tail canonical integral is defined by:

$$S(x,y,z)=\underset{-\infty }{\overset{\infty }{\int }}{e}^{i\left(u^5+x{u}^3+y{u}^2+zu\right)}du$$

(20)

As a further example, we consider a special case of the swallow-tail integral, i.e., $S(x,0,z)$—the oddoid integral of the order two [18]. For the real x and y the function $S(x,0,z)$ is also real, and for the numerical evaluation the equation $S(x,0,z)=2\underset{0}{\overset{\infty }{\int }}\cos \left(u^5+x{u}^3+zu\right)du$ is used. This integral is of interest in the study of bound-continuum [19] and bound-bound [20] Franck–Condon factors. The analysis is applied to the domain $x<0$.

In that case the phase function is $f(x,0,z;u)=u^5+x{u}^3+zu$ and it is antisymmetric with respect to the variable u: $f(x,0,z;-u)=-f(x,0,z;u)$. There are four saddle points $u_i$ defined by the condition $f^{(1)}(x,0,z;u)=5u^4+3x{u}^2+z=0$:

$$\begin{array}{l}{u}_1(x,z)=-u_4(x,z)=-\sqrt{\frac{3}{10}}\sqrt{\sqrt{x^2-\frac{20}{9}z}-x}\\ u_2(x,z)=-u_3(x,z)=-\sqrt{\frac{3}{10}}\sqrt{-x-\sqrt{x^2-\frac{20}{9}z}}\end{array}$$

(21)

The condition $f^{(2)}(x,0,z;u_e)=20u_e^3+6x{u}_e=0$ defines three real bifurcation points $u_{ei}$ of the “fold” type: $u_{e1}=-u_{e3}=-\sqrt{\frac{3\left|x\right|}{10}}$, $u_{e2}=0$. As there are two real points ($u_{p1}=-\sqrt{\frac{\left|x\right|}{10}}$, $u_{p2}=\sqrt{\frac{\left|x\right|}{10}}$) satisfying the conditions $f^{(3)}(x,z;u_{pi})=0$ and $\mathrm{sgn}\left(f^{(2)}(x,z;u_{pi})\right)=-\mathrm{sgn}\left(f^{(4)}(x,z;u_{pi})\right)$, according to Section 2, the integral $S(x,0,z)$ can be written as:

$$\begin{array}{ll}S(x,0,z)& ={\displaystyle \sum _{i=1}^3{I}_i(x,z)}-{\displaystyle \sum _{i=1}^2{I}_{pi}(x,z)}\\ & ={\displaystyle \sum _{i=1}^3\underset{-\infty }{\overset{\infty }{\int }}{e}^{i{f}_i(x,0,z;u)}du}-{\displaystyle \sum _{i=1}^2\underset{-\infty }{\overset{\infty }{\int }}{e}^{i{f}_{pi}(x,0,z;u)}du}\end{array}$$

(22)

where,

$$f_{p1}(x,0,z;u)=-f_{p2}(x,0,z;-u)=\frac{\sqrt{10}}{250}{\left(-x\right)}^{\frac{5}{2}}+u\left(z+\frac{3}{20}{x}^2\right)+u^2\frac{\sqrt{10}}{5}{\left(-x\right)}^{\frac{3}{2}},$$

$$f_1(x,0,z;u)=\{\begin{array}{cc}f(x,0,z;u)& u\le u_{p,1}\\ f_{p1}(x,0,z;u)& u>u_{p,1}\end{array}, f_2(x,0,z;u)=\{\begin{array}{cc}{f}_{p,1}(x,0,z;u)& u<u_{p,1}\\ f(x,0,z;u)& u_{p,1}\le u\le u_{p,2}\\ f_{p,2}(x,0,z;u)& u>u_{p,2}\end{array}, \mathrm{and}$$

$$f_3(x,0,z;u)=\{\begin{array}{cc}{f}_{p,2}(x,0,z;u)& u<u_{p,2}\\ f(x,0,z;u)& u\ge u_{p,2}\end{array}.$$

Since $f_{p1}(x,0,z;u)=-f_{p2}(x,0,z;-u)$, it follows that $I_{p2}(x,z)=I_{p1}{(x,z)}^{*}$ and ${\displaystyle \sum _{i=1}^2{I}_{pi}(x,z)}=2\mathrm{Re}{I}_{p1}(x,z)$. Also, $f_1(x,0,z;u)=-f_3(x,0,z;-u)$, $I_3(x,z)=I_1{(x,z)}^{*}$, and $I_1(x,z)+I_3(x,z)=2\mathrm{Re}{I}_1(x,z)$. Using Equation (9) to calculate $I_{p1}(x,z)$, one can write Equation (22) in the form:

$$S(x,0,z)=2\mathrm{Re}{I}_1(x,z)+I_2(x,z)-\sqrt[4]{\frac{40{\pi }^2}{{\left|x\right|}^3}}\cos \left(\frac{19x^4-600x^{2}z-2000z^2}{1600\sqrt{10}{\left|x\right|}^{3/2}}+\frac{\pi }{4}\right)$$

(23)

This expression exactly represents the function $S(x,0,z)$. To find an approximate solution of the integral $S(x,0,z)$ one needs to calculate integrals $I_1(x,z)$ and $I_2(x,z)$, using the approximation described by Equation (4).

The function $f_1(x,0,z;u)$ has two saddle points (see Figure 5b): ${\tilde{u}}_1(x,y)=u_1(x,y)$, ${\tilde{u}}_2(x,z)=\{\begin{array}{cc}{u}_2(x,y)& z\ge \frac{x^2}{4}\\ -\frac{5}{\sqrt{40}{\left|x\right|}^{\frac{3}{2}}}\left(z+\frac{3}{20}{x}^2\right)& z<\frac{x^2}{4}\end{array}$. Applying Equation (4) one gets,

$$I_1(x,z)\approx I_{F1}(x,z)=\pi e^{i{A}_1(x,z)}\left[\left(\frac{1}{\sqrt{10u_1^3+3x{u}_1}}+\frac{1}{\sqrt{-10{\tilde{u}}_2^3-3x{\tilde{u}}_2}}\right)\sqrt[4]{{\zeta }_1(x,z)}Ai(-{\zeta }_1(x,z))-i\left(\frac{1}{\sqrt{10u_1^3+3x{u}_1}}-\frac{1}{\sqrt{-10{\tilde{u}}_2^3-3x{\tilde{u}}_2}}\right)\frac{A{i}^{\prime }(-{\zeta }_1(x,z))}{\sqrt[4]{{\zeta }_1(x,z)}}\right]$$

(24)

where $A_1(x,z)=\frac{1}{2}\left[f_1(x,z;{\tilde{u}}_2)+f_1(x,z;u_1)\right]$ and ${\zeta }_1{(x,z)}^{3/2}=\frac{3}{4}\left[f_1(x,z;{\tilde{u}}_2)-f_1(x,z;u_1)\right]$.

The function $f_2(x,0,z;u)$ has two symmetrical saddle points: ${\tilde{u}}_3(x,z)=\{\begin{array}{cc}{u}_2(x,z)& z\le \frac{x^2}{4}\\ -\frac{\sqrt{10}}{40{\left|x\right|}^{\frac{3}{2}}}\left(z+\frac{3}{20}{x}^2\right)& z>\frac{x^2}{4}\end{array}$, ${\tilde{u}}_4(x,z)=-{\tilde{u}}_3(x,z)$ (Figure 5c). Since $f_2(x,z;{\tilde{u}}_3)=-f_2(x,z;{\tilde{u}}_4)$ and $f_2^{(2)}(x,z;{\tilde{u}}_3)=-f_2^{(2)}{}_2(x,z;{\tilde{u}}_4)$, the approximation of the integral $I_2(x,z)$ has a simple form,

$$I_2(x,z)\propto I_{F2}(x,z)=\frac{2\pi }{\sqrt{-10{\tilde{u}}_3^3-3x{\tilde{u}}_3}}\sqrt[4]{{\zeta }_2(x,z)}Ai(-{\zeta }_2(x,z))$$

(25)

where ${\zeta }_2{(x,z)}^{3/2}=-\frac{3}{2}{f}_2(x,z;{\tilde{u}}_3)$.

Finally, we write the approximation of the integral $S(x,0,z)$ as:

$$S_A(x,0,z)=2\mathrm{Re}{I}_{F1}(x,z)+I_{F2}(x,z)-\sqrt[4]{\frac{40{\pi }^2}{{\left|x\right|}^3}}\cos \left(\frac{19x^4-600x^{2}z-2000z^2}{1600\sqrt{10}{\left|x\right|}^{3/2}}+\frac{\pi }{4}\right)$$

(26)

In

Table 5. Comparison of the values of $S(x,0,z)$ and $S_A(x,0,z)$ on the caustics $z=0$ and $z=\frac{9}{20}{x}^2$.

x	z	S(x,0,0)	SA(x,0,0)	z = (9x2)/20	S(x,0,z)	SA(x,0,z)
−1	0.	2.269084	1.913266	0.45	2.043860	1.661334
−2	0.	2.409949	2.394455	1.8	1.663543	1.612423
−3	0.	0.596406	0.604199	4.05	−0.519597	−0.508896
−4	0.	1.292366	1.281152	7.2	−0.915491	−0.908131
−5	0.	0.215598	0.212031	11.25	0.821846	0.816979
−6	0.	0.247304	0.245451	16.2	0.660728	0.660188
−7	0.	0.670358	0.667798	22.05	−0.394882	−0.392698
−8	0.	0.834808	0.83448	28.8	0.204247	0.205210
−9	0.	1.155092	1.154554	36.45	1.055443	1.053662
−10	0.	0.767679	0.767563	45.	−0.891812	−0.890310

we compare the values of the functions $S(x,0,z)$ and $S_A(x,0,z)$ at the caustics i.e., at the points where the function $f^{(1)}(x,0,z;u)$ has an extreme. These comparisons together with the comparison of functions in Figure 6 clearly show that the function $S_A(x,0,z)$ is a good approximation of the function $S(x,0,z)$ if the condition $\sqrt{x^2+y^2}>3$ is satisfied.

4. Discussion and Conclusions

We have shown that the original oscillatory integral can be expressed exactly as a sum of integrals, each having either one or at most two real stationary points. The construction of these integrals introduces new phase functions that are smooth (infinitely differentiable), but only a few first derivatives are continuous in a characteristic point. This means that the proposed method is limited to the leading term only (i.e., the integral is the same as the leading term plus the residue that cannot be further treated as in the standard application of asymptotic analysis and iterated to obtain an asymptotic expansion [6]). However, the method has practical applications, especially in cases where the phase function (or its first derivative!) is tabulated.

The validity of the proposed method was tested on examples of integrals with three saddle points (“cusp” catastrophe) and four saddle points (“swallow-tail” catastrophe). The examples chosen are typical of the spectra of diatomic molecules, but the method described in this paper can be used for numerical computation of canonical integrals occurring in other physical fields as well, e.g., the propagation of electromagnetic, sound or fluid waves, and particularly within the semiclassical theory of atom–atom and atom–surface scattering, chemical reactions etc.