Three-Airy Beams, Their Propagation in the Fresnel Zone, the Autofocusing Plane Location, as Well as Generalizing Beams

Abstract

A family of 2D light fields consisting of the product of three Airy functions with linear arguments has been studied theoretically and experimentally. These fields, called three-Airy beams, feature a parameter shift and have a cubic phase and a super-Gaussian circular intensity in the far zone. Transformations of three-Airy beams in the Fresnel zone have been studied using theoretical, numerical, and experimental means. It has been shown that the autofocusing plane of a three-Airy beam is similar to the square root of the shift parameter. We also introduce generalized three-Airy beams containing nine free parameters, and obtain their Fourier transform in a closed form.

1. Introduction

Light fields based on the Airy function usage are intensively studied objects in modern optics. As a solution of the 1D paraxial equation, the Airy function first appeared in 1979, when the non-diffracting Airy beams of infinite energy were found [1]. However, a boom in optical research related to the Airy function began in 2007, when it was noted that a linear exponential factor allows the construction of paraxial Airy beams of finite energy [2]. In general, Airy beams depend on some real parameters and are functionally stable solutions of the paraxial equation. Despite the lack of structural stability, under certain restrictions, there is a propagation zone in which the intensity of Airy beams varies very weakly. In addition, the parabolic shape of the propagation trajectory of Airy beams provides new opportunities for the practical application of such beams since traditional Gaussian beams are characterized by rectilinear propagation.

This discovery was followed by studies and publications that looked at various aspects related to Airy beams. Among the theoretical works, the following results should be noted: the transformation of Airy–Gaussian beams in first-order optical systems [3], the Poynting vector and orbital angular momentum of Airy beams [4], the nonparaxial study of Airy beams [5,6], truncated Airy beams [7,8], the comparison of Airy beams and Hermite–Gaussian modes under astigmatic transformation [7], vortex Airy beams [7,9] and others [10,11,12,13,14]. Among the experimental works, we mention papers related to the study of the propagation of Airy beams in various media [15,16,17], as well as the use of Airy beams for problems related to microparticles [18] (see also the references therein).

As a rule, either 1D Airy beams or their 2D analogs, obtained by multiplying two 1D beams whose arguments are orthogonal Cartesian coordinates in the plane, were considered. Airy beams of more complex shape (for example, circular Airy beams) have been investigated only numerically and experimentally since theoretical study requires the discovery of quite nontrivial integrals in a closed form. However, the theory of integral transformations of Airy functions is not yet sufficiently developed to cope with such challenges.

One of the nontrivial examples when theoretical research turned out to be possible is a family of three-Airy beams [19]. In 2010, these beams were introduced as the product of three 1D Airy patterns, which are shifted from the origin to the vertices of an equilateral triangle and rotated so that the resulting field is invariant under rotation by ${120}^{\circ }$:

$${\mathrm{Ai}}_3\left(\mathbf{r}\right|c)=\mathrm{Ai}\left(\frac{-x\sqrt{3}-y}{2}+c\right)\mathrm{Ai}\left(\frac{x\sqrt{3}-y}{2}+c\right)\mathrm{Ai}(y+c),$$

(1)

where

$$\mathrm{Ai}\left(x\right)=\frac{1}{2\pi }{\int }_{\mathbb{R}}exp\left(\frac{\mathrm{i}{t}^3}{3}+\mathrm{i}xt\right)\mathrm{d}t$$

(2)

is the Airy function definition, $\mathbf{r}=(x,y)$ is a 2D vector, and c is a shift parameter.

The shift parameter determines the shape of the three-Airy beams and their behavior when propagating in the Fresnel zone. As it turned out, the Fourier transform of these beams has a cubic phase and a circular super-Gaussian intensity [19]:

$$\mathcal{F}\left[{\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)\right]\left(\mathbf{r}\right)=\frac{1}{3^{5/6}\pi b^2}exp\left(-\frac{2\mathrm{i}(3x^{2}y-y^3)}{27b^3}\right)\mathrm{Ai}\left(3^{2/3}c+\frac{{2\left|\mathbf{r}\right|}^2}{3^{4/3}{b}^2}\right).$$

(3)

Here, $\mathit{\rho }=(\xi ,\eta )$ is a 2D vector, b is a scaling factor,

$$\mathcal{F}\left[f\left(\mathit{\rho }\right)\right]\left(\mathbf{r}\right)=\frac{1}{2\pi }{\int \int }_{{\mathbb{R}}^2}{\mathrm{e}}^{-\mathrm{i}\langle \mathbf{r},\mathit{\rho }\rangle }f\left(\mathit{\rho }\right){\mathrm{d}}^2\mathit{\rho }$$

(4)

is 2D Fourier transform, and $\langle \mathbf{r},\mathit{\rho }\rangle =x\xi +y\eta$ is the standard scalar product of 2D vectors. As is seen in (3), only the Airy factor depends on the shift parameter. A simple proof of the relation (3) is presented in Appendix A.

The most interesting cases are when, at the central point, $x=y=0$, the Airy function reaches a maximum (i.e., $3^{2/3}c=a_n^{\prime }$) or goes to zero (i.e., $3^{2/3}c=a_n$). Here $a_n$ and $a_n^{\prime }$ are zeros of the functions $\mathrm{Ai}\left(x\right)$ and ${\mathrm{Ai}}^{\prime }\left(x\right)$, respectively. All of them are real and located on the negative semi-axis. Asymptotic expansions are known for both families of zeros [20],

$$a_n\approx -{\left(\frac{3\pi }{8}(4n-1)\right)}^{2/3},a_n^{\prime }\approx -{\left(\frac{3\pi }{8}(4n-3)\right)}^{2/3}.$$

(5)

Even for small n, the initial terms of these expansions give a very good approximation to the exact values of the zeros. For example,

$$\begin{array}{c}{a}_1=-2.338\dots \approx -2.320,a_1^{\prime }=-1.018\dots \approx -1.115,\\ a_2=-4.087\dots \approx -4.081,a_2^{\prime }=-3.248\dots \approx -3.261,\\ a_3=-5.520\dots \approx -5.517,a_3^{\prime }=-4.820\dots \approx -4.826,\\ a_4=-6.786\dots \approx -6.784,a_4^{\prime }=-6.163\dots \approx -6.167,\end{array}$$

(6)

where the first number is the exact value of zero, and the second is due to Equation (5).

In subsequent years, three-Airy beams were constructed and studied in various ways. Torre in [21] found their intensity moments of second order. In [22], the three-Airy beam for the case $3^{2/3}c=a_1^{\prime }$ was implemented on the basis of the Fourier transform method. Its propagation in free space and a nonlinear medium (SBN crystal) was studied experimentally. In [23], Desyatnikov et al. also used the Fourier transform method, but the initial field was formed by the holographic method. In that article, two beams with parameters $3^{2/3}c=a_3$ and $3^{2/3}c=a_3^{\prime }$ were investigated, and autofocusing of three-Airy beams was discussed. In addition, Desyatnikov considered the propelling behavior of such beams when an optical singularity is embedded, as well as the prospects for their use for transporting particles. In [24], vortex three-Airy beams with different topological charges have been investigated.

In this work, we continue to study three-Airy beams by theoretical and numerical means, as well as in optical experiments. In Section 2, we discuss the ways to evaluate the Fresnel transform of these beams. Here, the experimental implementation of three-Airy beams is also considered. In Section 3, we investigate the autofocusing behavior of these beams. In Section 4, we generalize the relation (3) and provide some corollaries concerning the beams of finite and infinite energy. In Conclusion, we discuss possible applications of the results obtained.

We consider the study of three-Airy beams at various values of the shift parameter, the characteristics of their propagation in the Fresnel zone, and the determination of the autofocusing plane depending on the shift parameter as a necessary basis for the subsequent transition to the study of multiple beams of Pearcey, swallowtail, and other diffraction catastrophes, whose non-Gaussian nature also leads to the appearance of an autofocusing plane when propagating in the Fresnel zone.

2. Three-Airy Beams in the Fresnel Diffraction Zone

2.1. Theory

If a coherent light field of a wavelength $\lambda$ is specified by its complex amplitude $F_0\left(\mathbf{r}\right)$ in the initial plane $z=0$, then the field propagation in free space is described in the paraxial approximation by the equation $\{{\partial }_x^2+{\partial }_y^2+2\mathrm{i}k{\partial }_z\}F(\mathbf{r},z)=0$, where $F(\mathbf{r},0)=F_0\left(\mathbf{r}\right)$. If $F_0\left(\mathbf{r}\right)$ vanishes sufficiently rapidly as $\left|\mathbf{r}\right|\to \infty$, then $F(\mathbf{r},z)$ related to $F_0\left(\mathbf{r}\right)$ by the Fresnel transform [25]:

$$\begin{array}{cc}\hfill F(\mathbf{r},z)={\mathbf{FR}}_z\left[F_0\left(\mathit{\rho }\right)\right]\left(\mathbf{r}\right)& =\frac{k}{2\pi \mathrm{i}z}{\int \int }_{{\mathbb{R}}^2}exp\left(\frac{\mathrm{i}k}{2z}{|\mathbf{r}-\mathit{\rho }|}^2\right)F_0\left(\mathit{\rho }\right){\mathrm{d}}^2\mathit{\rho }.\hfill \end{array}$$

(7)

where $k=2\pi /\lambda$ is the wavenumber of light.

The numerical evaluation of the double integral (7) can be carried out by expanding the initial field into a series of Hermite–Gaussian (HG) modes (see Appendix B).

When choosing $F_0\left(\mathit{\rho }\right)={\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)$, the question arises: what Gaussian beam width w is optimal for approximating a three-Airy beam by the superposition of HG modes (A11) if the double series is replaced by a truncated finite sum?

$$\begin{array}{c}\hfill {\mathbf{funcFR}}_{z,w}\left[{\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)\right]\left(\mathbf{r}\right)\cong \sum _{\begin{array}{c}0\le n\le N\\ 0\le m\le M\end{array}}{c}_{n,m}{\mathrm{e}}^{-\mathrm{i}(n+m)arg\sigma }{\mathcal{H}}_{n,m}\left(\frac{\mathbf{r}}{w}\right).\end{array}$$

(8)

This problem can be solved numerically using the least squares method (see [19] for detail), selecting the shift parameter $c=3^{-2/3}{a}_1^{\prime }$, for which the Fourier image of the three-Airy beam has a circular intensity very similar to Gaussian. The solution is $w=1.027/b$, and this value remains optimal even if $3^{2/3}c\ne a_1^{\prime }$.

Figure 1 depicts the intensity and phase distributions of the beam (8) for the shift parameters $c=3^{-2/3}{a}_3^{\prime }$ and $c=3^{-2/3}{a}_3$. Other parameters are $b=1$, $w=1$ and $z\in [0,\pi /2]$. As can be seen, when the three-Airy beam propagates, there is a plane in which the intensity has a trefoil shape with a strong central peak. We consider this plane to be the autofocusing (AF) plane.

For the initial beam $F_0\left(\mathit{\rho }\right)={\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)$, we were unable to calculate analytically the double integral (7), while we reduce it to a 1D integral:

$$\begin{array}{cc}\hfill & {\mathbf{FR}}_z\left[{\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)\right]\left(\mathbf{r}\right)\hfill \\ \hfill & =\frac{1}{\sqrt{\pi }}exp\left(\frac{9\mathrm{i}c{b}^{2}z}{4k}+\frac{27\mathrm{i}{b}^6{z}^3}{64k^3}\right){\int }_{\mathbb{R}}{\mathrm{e}}^{-t^2}{\mathrm{Ai}}_3(b\mathbf{r}\left|c+\frac{9b^4{z}^2}{16k^2}+\sqrt{\frac{b^{2}z}{k}}·{\mathrm{e}}^{-\pi \mathrm{i}/4}t\right)\mathrm{d}t.\hfill \end{array}$$

(9)

This integral representation is useful in the theoretical study of the elliptic umbilic catastrophe [26], but it is not applicable for numerical simulation related to the propagation of three-Airy beams in the Fresnel zone. The reason is that the integrand in (9) is the product of a rapidly decreasing Gaussian function, ${\mathrm{e}}^{-t^2}$, and a rapidly increasing Airy function of the complex argument, ${\mathrm{Ai}}_3$∼${\mathrm{e}}^{{\left|t\right|}^{3/2}}$ as $t\to -\infty$. Of course, the integral (9) converges absolutely. However, to obtain acceptable calculation accuracy, we are forced to make the integration domain too large, which leads to an excessive increase in calculation time. As a result, for numerical simulation, the use of the series expansion in terms of HG modes (8) is preferable.

2.2. Experiment

The experimental part of our investigation of three-Airy beams has been conducted as follows. Figure 2 shows the optical set-up. We use solid-state laser L with a wavelength of 532 nm. Lenses $L_1$ and $L_2$ ensure beam expansion to the desired size. Lens $L_3$ forms intensity distribution on the sensitive area of the video camera. To reduce the formation distance, a lens with a focal length of 1000 mm was used.

The amplitude–phase distribution presented in Figure 3a,b was encoded into a pure phase one by partial encoding method [27]. The partial encoding technique makes it possible to vary the error of the field generation and the diffraction efficiency in a wide range, choosing the best ratio for a specific task. A threshold value $\alpha \in [0,1]$ is proportional to the accuracy of the formation of a given field and inversely proportional to efficiency. The simulation showed a slight difference in the formed patterns calculated for $\alpha =1$ and $\alpha =0.5$, so in the experiments, we used an optical element for $\alpha =0.5$ (see Figure 3d), because in this case, higher diffraction efficiency is ensured.

The DOE shown in Figure 3d was fabricated using the method photolithography on a quartz plate with a resolution of 1 $\mathsf{\mu }$m, and the microrelief depth was 530 nm. Figure 4 shows the shape of the surface microrelief section DOE. The diameter of the DOE is 4 mm.

The three-Airy beam propagation simulation and corresponding experimentally recorded intensity patterns obtained for phase element (Figure 3d), illuminated by a converging beam (with a focus at a distance of 1000 mm from the input plane), are shown in Figure 5.

3. Autofocusing Plane of Three-Airy Beams

During free-space propagation, many non-Gaussian paraxial light fields demonstrate an autofocusing (AF) behavior. It means that a light beam propagating in the Fresnel zone and with a wide intensity profile at some point suddenly changes sharply and takes the form of a single, highly localized peak of large magnitude.

The most famous beams with the AF property are circular Airy beams ${\mathrm{e}}^{a(r_0-r)}\mathrm{Ai}(r_0-r)$, where $r_0$ is the radius of the Airy ring on the input plane and a is an exponential truncation factor [28,29] (see also [30] and the references therein). Due to the radially symmetrical shape of such fields, the study of the location of their AF plane when propagating in the Fresnel zone is reduced to a 1D problem and allows the application of asymptotical methods based on the idea of the stationary phase.

For three-Airy beams, ${\mathrm{Ai}}_3\left(\mathbf{r}\right|c)$, even the very definition of the concept of “autofocusing plane” becomes a problem since the transverse intensity distribution of such beams is not radially symmetric and, as a consequence, a 3D profile of its intensity is not bell-shaped. However, in the simplest version, when the shift parameter is defined by the relation $3^{2/3}c=a_1^{\prime }\approx -1.018$, the intensity of the three-Airy beam in the Fresnel zone changes weakly, transforming from triangular to circular shape. In this case, the intensity reaches a maximum value on the axis of the beam propagation.

Thus, there are two ways to find the location of the AF plane of the three-Airy beam. First, one can write an equation for its intensity extrema points located on the optical axis (i.e., at $\mathbf{r}=\mathbf{0}$):

$$\begin{array}{c}\hfill {\partial }_z{\left|{\mathrm{Ai}}_3(\mathbf{0},z|c)\right|}^2=0\end{array}$$

(10)

and try to solve it within the framework of some additional asymptotic assumptions. Second, one can represent the Fresnel transform of the three-Airy beam as a certain integral of the pure phase exponential and apply the stationary phase method to find the caustic. We consider both of them.

3.1. Three-Airy Beam Intensity at the Points of the Optical Axis

It is known that the Fresnel transform (7) of a field $F_0\left(\mathit{\rho }\right)$ can be written in terms of the Fourier transform of the field. Specifically, since

$$F_0\left(\mathit{\rho }\right)={\mathcal{F}}^{-1}\left[\mathcal{F}\left[F_0\right]\left(\mathbf{r}\right)\right]\left(\mathit{\rho }\right)=\frac{1}{2\pi }{\int \int }_{{\mathbb{R}}^2}{\mathrm{e}}^{\mathrm{i}\langle \mathbf{r},\mathit{\rho }\rangle }\mathcal{F}\left[F_0\right]\left(\mathbf{r}\right){\mathrm{d}}^2\mathbf{r},$$

(11)

then

$${\mathbf{FR}}_z\left[F_0\left(\mathit{\rho }\right)\right]\left(\mathbf{r}\right)={\mathcal{F}}^{-1}\left[exp\left(-\frac{\mathrm{i}z}{2k}{\left|\mathit{\rho }\right|}^2\right)\mathcal{F}\left[F_0\right]\left(\mathit{\rho }\right)\right]\left(\mathbf{r}\right).$$

(12)

Substituting the three-Airy beam $F_0\left(\mathit{\rho }\right)={\mathrm{Ai}}_3\left(b\mathit{\rho }\right|c)$ and its Fourier image (3), one obtains

$${\mathrm{Ai}}_3(b\mathbf{r},z|c)=\frac{1}{3^{5/6}\pi b^2}{\mathcal{F}}^{-1}\left[exp\left(-\frac{\mathrm{i}z}{2k}{\left|\mathit{\rho }\right|}^2-\frac{2\mathrm{i}(3{\xi }^2\eta -{\eta }^3)}{27b^3}\right)\mathrm{Ai}\left(3^{2/3}c+\frac{{2\left|\mathit{\rho }\right|}^2}{3^{4/3}{b}^2}\right)\right]\left(\mathbf{r}\right).$$

(13)

Let us simplify this integral for $\mathbf{r}=\mathbf{0}$, passing to polar coordinates, $\xi +\mathrm{i}\eta =\rho {\mathrm{e}}^{\mathrm{i}\varphi }$, integrating over the polar angle and introducing the variable $t={\rho }^2/9b^2$:

$$\begin{array}{cc}\hfill {\mathrm{Ai}}_3(\mathbf{0},z|c)& =\frac{3^{7/6}}{2\pi }{\int }_0^{\infty }exp\left(-\frac{9\mathrm{i}z{b}^2}{2k}t\right)J_0\left(2t^{3/2}\right)\mathrm{Ai}\left(3^{2/3}[c+2t]\right)\mathrm{d}t.\hfill \end{array}$$

(14)

Then the equation ${\partial }_z{\left|{\mathrm{Ai}}_3(\mathbf{0},z|c)\right|}^2=0$ takes the form

$$\begin{array}{c}\hfill \mathrm{Im}\left\{G_1^{*}·G_0\right\}=0,\mathrm{where}{G}_m={\int }_0^{\infty }exp\left(\frac{9\mathrm{i}z{b}^2}{2k}t\right)J_0\left(2t^{3/2}\right)\mathrm{Ai}\left(3^{2/3}[c+2t]\right)t^m\mathrm{d}t,\end{array}$$

(15)

and the asterisk means complex conjugation.

Some intensity curves, $|{\mathrm{Ai}}_3{(\mathbf{0},z|c)|}^2$, are depicted in Figure 6 and Figure 7. Here, we use the dimensionless variable $z:=2z{b}^2/k$ and select the shift parameter c so that the Fourier image of the three-Airy beam reaches its maximum value at a point lying on the optical axis, $c=c_n=3^{-2/3}{a}_n^{\prime }$. Thus, the location of the AF plane $z_s\left(c_n\right)$ depends on n. A numerical study of the discrete function $z_s\left(c_n\right)$, where $20\le n\le 150$, revealed that $z_s\left(c_n\right)\approx \sqrt{n}+5$ for large n (in particular, $z_s\left(c_{60}\right)-\sqrt{60}=4.812$), while, of course, the numerical results cannot be considered to be a true asymptotic behavior of $z_s\left(c\right)$ as $c\to \infty$.

3.2. Caustic of Three-Airy Beam upon Propagation in the Fresnel Zone

Finding the caustic is one of the traditional problems that arise when the properties of diffraction catastrophe integrals and light fields associated with them are studied. Each diffraction catastrophe is connected to a definite caustic due to the equations on the stationary phase and the Hessian. A good introduction to this field (catastrophe optics) is the paper [31] (see also [32] and the references therein). Recently, caustical investigation has been applied to aberration laser beams [33], vortex beams [34], and Ince–Gaussian beams [35].

Let us consider the problem of finding the AF plane of the beam ${\mathrm{Ai}}_3(\mathbf{r},z|c)$ as the process of the appearance of a caustic point of the beam when it propagates in the Fresnel zone. Replacing the radial Airy function in Equation (13) with the definition (2) and changing the variable $t\to 3^{-2/3}t$, we obtain

$$\begin{array}{cc}\hfill {\mathrm{Ai}}_3(\mathbf{r},z|c)& =\frac{3^{-5/6}}{4{\pi }^3}{\int \int \int }_{{\mathbb{R}}^3}exp\left(\mathrm{i}P(\xi ,\eta ,t)\right)\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}t,\hfill \end{array}$$

(16)

where $P(\xi ,\eta ,t)=(x\xi +y\eta )-\frac{1}{4}z({\xi }^2+{\eta }^2)-\frac{2}{27}(3{\xi }^2\eta -{\eta }^3)+\frac{1}{27}{t}^3+t\left[c+\frac{2}{9}({\xi }^2+{\eta }^2)\right]$.

By making the substitution of the output variables $(x,y,z,c)=\frac{4}{9}\left(X,Y,2Z,-A\right)$ in the function $P(\xi ,\eta ,t)$,

$$\begin{array}{c}\hfill P(\xi ,\eta ,t)=\frac{4}{9}(X\xi +Y\eta -At)+\frac{2}{9}(t-Z)({\xi }^2+{\eta }^2)+\frac{1}{27}(t^3-6{\xi }^2\eta +2{\eta }^3),\end{array}$$

(17)

and equating its gradient to zero, ${\partial }_{\xi }P={\partial }_{\eta }P={\partial }_{t}P=0$, we obtain a system for finding the stationary points:

$$\begin{array}{c}\xi \eta +(Z-t)\xi =X,\frac{1}{2}({\xi }^2-{\eta }^2)+(Z-t)\eta =Y,({\xi }^2+{\eta }^2)+\frac{1}{2}{t}^2=2A.\end{array}$$

(18)

From the last equation, it is seen that a real solution $(\xi ,\eta ,t)$ exists for $A>0$ only.

Equating the Hessian of $P(\xi ,\eta ,t)$ to zero, we have a fourth relationship between the variables:

$$(t^3-6{\xi }^2\eta +2{\eta }^3)+2Z({\xi }^2+{\eta }^2-t^2)+t(Z^2-3[{\xi }^2+{\eta }^2])=0.$$

(19)

This equation, when jointly solved with the system (18), allows us to find the location and shape of the caustic curve.

Ideally, our further plan of action is as follows. We find an explicit solution of the system (18) in the form of $\xi =\xi (X,Y,Z,A)$, $\eta =\eta (X,Y,Z,A)$, and $t=t(X,Y,Z,A)$. Then, we substitute it into (19) and obtain a caustic equation, $\Phi (X,Y,Z,A)=0$, which for the values $X=Y=0$ leads to AF plane: $\Phi (0,0,Z,A)=0\Rightarrow Z=Z\left(A\right)$.

Carrying out this plan and passing from Cartesian coordinates to polar ones, $(\xi ,\eta )=(\rho cos\varphi ,\rho sin\varphi )$, we obtain

$$\begin{array}{c}X=\frac{1}{2}{\rho }^{2}sin2\varphi +(Z-t)\rho cos\varphi ,Y=\frac{1}{2}{\rho }^{2}cos2\varphi +(Z-t)\rho sin\varphi ,{\rho }^2+\frac{1}{2}{t}^2=2A,\\ (t^3-2{\rho }^{3}sin3\varphi )+2Z({\rho }^2-t^2)+t(Z^2-3{\rho }^2)=0.\end{array}$$

Later, we use the equation ${\rho }^2+\frac{1}{2}{t}^2=2A$ to replace $\rho =\sqrt{2A}cos\theta$, $t=2\sqrt{A}sin\theta$, where $-\pi /2\le \theta \le \pi /2$. In fact, it makes the transition from Cartesian to spherical coordinates $(\xi ,\eta ,t)\to (\rho ,\varphi ,\theta )$.

Making another change in variables, $(X,Y,Z)=(A\tilde{X},A\tilde{Y},\sqrt{A}\tilde{Z})$, we obtain the caustic equations without the explicit presence of the variable A (thus, the use of the stationary phase method at $A\gg 1$ becomes justified):

$$\begin{array}{c}\tilde{X}={cos}^2\theta sin2\varphi +\sqrt{2}(\tilde{Z}-2sin\theta )cos\theta cos\varphi ,\end{array}$$

(20)

$$\begin{array}{c}\tilde{Y}={cos}^2\theta cos2\varphi +\sqrt{2}(\tilde{Z}-2sin\theta )cos\theta sin\varphi ,\end{array}$$

(21)

$$\begin{array}{c}2\sqrt{2}{cos}^3\theta sin3\varphi =sin\theta {(\tilde{Z}-2sin\theta )}^2+2{cos}^2\theta (\tilde{Z}-3sin\theta ).\end{array}$$

(22)

From the last equation, we find the solution $\varphi =\varphi (\tilde{Z},\theta )$. Substituting it into Equations (20) and (21), we obtain the caustic curve in parametric form.

The case $cos\theta =0$ is the most interesting. If $\theta =-\pi /2$, then there is no solution due to the constraint $z\ge 0$. If $\theta =\pi /2$, then we obtain the singular solution $(\tilde{X},\tilde{Y},\tilde{Z})=(0,0,2)$. This is a caustic point. It determines the location of the AF plane. Returning to the original variables $(x,y,z)$, we obtain $x=y=0$ and

$$z_s=\frac{8}{3}\sqrt{-c}.$$

(23)

Thus, if the shift parameter c, while remaining negative, increases in absolute value, then the AF plane moves away from the initial plane.

It is interesting to note that the same result is obtained if we replace the Bessel function in Equation (15) with its asymptotic expression, $J_0\left(x\right)$∼$\sqrt{2/\pi x}cos(x-\pi /4)$. However, this way requires much more effort.

In the numerical experiments above, we considered the case $c=c_n=3^{-2/3}{a}_n^{\prime }$. Then,

$$z_s\left(c_n\right)=\frac{8}{3}\sqrt{-3^{-2/3}{a}_n^{\prime }}\stackrel{\left(5\right)}{\approx }\frac{4}{3}{\left[\pi (4n-3)\right]}^{1/3}.$$

(24)

The exponent of n in the last formula does not coincide with the analogous indicator in the dependence $z_s\left(c_n\right)\approx \sqrt{n}+5$, found empirically for $n\le 150$ in numerical experiments. Nevertheless, already at $n=100$, we obtain a quite good correspondence with an exact numerical result (see Figure 6): $z_s\left(c_{100}\right)\stackrel{\left(24\right)}{=}14.3522\approx 15.012$. Taking into account the asymptotic nature of the formula (23), we expect that its accuracy will increase with increasing n, while the accuracy of the empirical formula will decrease.

Fixing $\tilde{Z}=2$ (the AF plane), let us investigate the question: what is the shape of the caustic curve in the plane of variables $(\tilde{X},\tilde{Y})$? We omit the case $cos\theta =0$, leading to the caustic point. Then,

$$\begin{array}{c}\hfill sin3\varphi =\frac{sin\theta {(2-2sin\theta )}^2+2{cos}^2\theta (2-3sin\theta )}{2\sqrt{2}{cos}^3\theta }.\end{array}$$

(25)

The curve of the RHS of (25) versus $\theta \in [-\pi /2,\pi /2]$ is shown in Figure 8. As is seen, there are two intervals, $\theta \in [-0.651,-0.418]$ and $[0.339,1.123]$, for which the RHS value is placed in $[-1,1]$, i.e., for which a solution ${\varphi }_0=\varphi \left(\theta \right)$ exists. The general solution $\varphi \left(\theta \right)$ is a set of the values $\left\{{\varphi }_0,{\varphi }_0+\frac{2\pi }{3},{\varphi }_0+\frac{4\pi }{3}\right\}$ and $\left\{\frac{\pi }{3}-{\varphi }_0,\pi -{\varphi }_0,\frac{5\pi }{3}-{\varphi }_0\right\}$. If we use these six values $\varphi \left(\theta \right)$ to construct the parametric curve given by Equations (20) and (21), then we obtain a caustic, consisting of two parts:

• outside—a triangular-shaped curve (six values for $\theta \in [-0.651,-0.418]$ are used),
• inside—a 3-bladed curve with sharp vertices (six values for $\theta \in [0.339,1.123]$ are used).

The shape of the caustic (shown on the right in Figure 8) is very similar to what is obtained in numerical and optical experiments.

4. Fourier Transform of Generalized Three-Airy Beams

In this section, we consider a 2D light field generalizing the three-Airy beam (1). Let us define the field as the product of three Airy functions with linear arguments of general form,

$${\mathrm{Ai}}_3\left(\mathbf{r}\right|\mathbf{M})=\prod _{1\le n\le 3}\mathrm{Ai}(A_{n}x+B_{n}y+C_n),$$

(26)

where $A_n$, $B_n$, $C_n$ are real parameters combined into a matrix

$$\mathbf{M}=\left(\begin{array}{ccc}{A}_1& A_2& A_3\\ B_1& B_2& B_3\\ C_1& C_2& C_3\end{array}\right).$$

(27)

We do not discuss what values these parameters can take. The only restriction is that each of the Airy multipliers is nontrivial (i.e., we discard cases like $A_1=B_1=0$).

The main result of this section is that the Fourier image of the field (26) can be found in a closed form:

$$\begin{array}{c}\hfill \mathcal{F}\left[{\mathrm{Ai}}_3\left(\mathit{\rho }\right|\mathbf{M})\right]\left(\mathbf{r}\right)=\frac{1}{{2\pi |\Delta |}^{1/3}}exp(\mathrm{i}{P}_1+\mathrm{i}{P}_3)\mathrm{Ai}(P_0+P_2),\end{array}$$

(28)

where $\Delta$ is a constant and $P_n=P_n\left(\mathbf{r}\right)$ are homogeneous real-valued polynomials of degree n. All of them are quite cumbersome since the LHS of (28) contains 9 free parameters. We write $\Delta$ and $P_n$ using auxiliary quantities. Specifically, let

$${\Delta }_1=\left|\begin{array}{cc}{A}_2& A_3\\ B_2& B_3\end{array}\right|,{\Delta }_2=\left|\begin{array}{cc}{A}_3& A_1\\ B_3& B_1\end{array}\right|,{\Delta }_3=\left|\begin{array}{cc}{A}_1& A_2\\ B_1& B_2\end{array}\right|,$$

(29)

that means that $det\mathbf{M}=C_1{\Delta }_1+C_2{\Delta }_2+C_3{\Delta }_3$. Let also $R_n=B_{n}x-A_{n}y$. Then,

$$\begin{array}{c}\Delta ={\Delta }_1^3+{\Delta }_2^3+{\Delta }_3^3,P_0=\frac{det\mathbf{M}}{{\Delta }^{1/3}},P_2=\frac{{\Delta }_1{\Delta }_2{\Delta }_3}{{\Delta }^{4/3}}\sum _{1\le n\le 3}\frac{R_n^2}{{\Delta }_n},\\ P_1=\frac{1}{\Delta }\left|\begin{array}{ccc}{R}_1& R_2& R_3\\ {\Delta }_1^2& {\Delta }_2^2& {\Delta }_3^2\\ C_1& C_2& C_3\end{array}\right|,P_3=\frac{1}{3{\Delta }^2}\left\{\left|\begin{array}{ccc}{R}_1^3& R_2^3& R_3^3\\ {\Delta }_1^3& {\Delta }_2^3& {\Delta }_3^3\\ 1& 1& 1\end{array}\right|+\prod _{1\le m<n\le 3}(R_n{\Delta }_n-R_m{\Delta }_m)\right\}.\end{array}$$

(30)

Formula (28) can be proved at least in two ways. The first is based on the 1D Fourier image of the product of two Airy functions ([36] Section 3.5.3): If

$$F\left(x\right)={\int }_{\mathbb{R}}{\mathrm{e}}^{-\mathrm{i}xt}\mathrm{Ai}(A_{1}t+B_1)\mathrm{Ai}(A_{2}t+B_2)\mathrm{d}t,$$

(31)

then

$$\begin{array}{cc}\hfill F\left(x\right)& =\frac{1}{|A_1^3-A_2^3{|}^{1/3}}exp\left(\frac{\mathrm{i}(A_1^2{B}_1-A_2^2{B}_2)x}{A_1^3-A_2^3}+\frac{\mathrm{i}(A_1^3+A_2^3)x^3}{3{(A_1^3-A_2^3)}^2}\right)\times \hfill \\ \hfill & \times \mathrm{Ai}\left(\frac{1}{{(A_1^3-A_2^3)}^{1/3}}\left[A_1{B}_2-A_2{B}_1-\frac{A_1{A}_2{x}^2}{A_1^3-A_2^3}\right]\right),A_1\ne A_2\hfill \\ \hfill F(x\ne 0)& =\frac{1}{2\sqrt{\pi \left|Ax\right|}}exp\left(\frac{\mathrm{i}{x}^3}{12A^3}+\frac{\mathrm{i}(B_2+B_1)x}{2A}-\frac{\mathrm{i}A{(B_2-B_1)}^2}{4x}+\frac{\pi \mathrm{i}}{4}\mathrm{sgn}\left(Ax\right)\right),A_1=A_2=A\ne 0\hfill \\ \hfill F\left(0\right)& =\frac{1}{\left|A\right|}\delta (B_2-B_1).A_1=A_2=A\ne 0\hfill \end{array}$$

(32)

The second way is similar to the proof of Equation (3) in Appendix A but lengthy. Nevertheless, both ways are straightforward. Moreover, a desire to find the representations of the polynomials $P_n$ in (28) which demonstrate their invariance under permutations of the triples $\{A_n,B_n,C_n\}$, requires quite tedious algebra.

The formula (28) has numerous and important corollaries. Some of them, interesting for mathematics, are mentioned in [37]. For optics, the light field ${\mathrm{Ai}}_3\left(\mathbf{r}\right|\mathbf{M})$ is a field of finite or infinite energy provided $\Delta >0$ or $\Delta <0$, respectively. The condition $\Delta <0$ can be expressed in terms of the relative position of the vectors $(A_n,B_n)$, $1\le n\le 3$, on the plane: If there is a straight line passing through the origin such that all three vectors are located along one side from it, then $\Delta <0$. And vice versa, if such a line does not exist, then $\Delta >0$.

Consider an example: Let the product of Airy functions in Equation (26) be as follows:

$${\mathrm{Ai}}_3\left(\mathbf{r}\right|\mathbf{M})=\mathrm{Ai}(-xsin\alpha +ycos\alpha +C_1)\mathrm{Ai}(xsin\alpha +ycos\alpha +C_2)\mathrm{Ai}(y+C_3).$$

(33)

Then

$$\begin{array}{c}(A_1,B_1)=(-sin\alpha ,cos\alpha ),(A_2,B_2)=(sin\alpha ,cos\alpha ),(A_3,B_3)=(0,1),\\ {\Delta }_1=-2cos\alpha sin\alpha ,{\Delta }_2={\Delta }_3=sin\alpha ,\Delta =2{sin}^3\alpha (1-4{cos}^3\alpha ),\\ R_1=xcos\alpha +ysin\alpha ,R_2=xcos\alpha -ysin\alpha ,R_3=x,\end{array}$$

(34)

and the polynomials $P_n$ are reduced to simpler expressions:

$$\begin{array}{cc}\hfill P_0& =\frac{C_1+C_2-2C_{3}cos\alpha }{2^{1/3}{(1-4{cos}^3\alpha )}^{1/3}},\hfill \\ \hfill P_1& =\frac{C_2-C_1}{2sin\alpha }·x+\frac{C_3-2(C_1+C_2){cos}^2\alpha }{1-4{cos}^3\alpha }·y,\hfill \\ \hfill P_2& =\frac{1}{2^{1/3}{(1-4{cos}^3\alpha )}^{1/3}}\left(\frac{1}{2{sin}^2\alpha }·x^2-\frac{2cos\alpha }{1-4{cos}^3\alpha }·y^2\right),\hfill \\ \hfill P_3& =-\frac{{cos}^2\alpha }{{sin}^2\alpha (1-4{cos}^3\alpha )}·x^{2}y+\frac{1+4{cos}^3\alpha }{3{(1-4{cos}^3\alpha )}^2}·y^3.\hfill \end{array}$$

(35)

One of the most interesting cases is when the RHS of Equation (28) contains a circular Airy factor. For this, the polynomial $P_2$ should be radially symmetric, which leads to the equation

$$\frac{1}{2{sin}^2\alpha }=-\frac{2cos\alpha }{1-4{cos}^3\alpha }\Rightarrow 8{cos}^3\alpha -4cos\alpha -1=0.$$

(36)

Introducing the variable $t=2cos\alpha$, one obtains the cubic equation $t^3-2t-1=0$. Its roots are $t=\{-1,(1\pm \sqrt{5})/2\}$ that corresponds to $\alpha =\{2\pi /3,\pi /5,3\pi /5\}$. We discuss these three cases separately.

(i)

If $\alpha =2\pi /3$, then $\Delta =3^{5/2}/2^3$, the field (33) is of finite energy, and

$$\begin{array}{cc}\hfill \mathcal{F}\left[{\mathrm{Ai}}_3\left(\mathit{\rho }\right|\mathbf{M})\right]\left(\mathbf{r}\right)& =\frac{1}{3^{5/6}\pi }exp(\frac{\mathrm{i}}{3}\left\{\sqrt{3}(C_2-C_1)x-(C_1+C_2-2C_3)y\right\}\hfill \end{array}$$

$$\begin{array}{c}\hfill -\frac{2\mathrm{i}}{27}(3x^{2}y-y^3))\mathrm{Ai}\left(\frac{C_1+C_2+C_3}{3^{1/3}}+\frac{2}{3^{4/3}}(x^2+y^2)\right).\end{array}$$

(37)

When $C_1=C_2=C_3=c$, this formula coincides with Equation (3).

(ii)

If $\alpha =\pi /5$, then

$$\begin{array}{c}cos\alpha =\frac{1+\sqrt{5}}{4},sin\alpha =\frac{1+\sqrt{5}}{4}\sqrt{5-2\sqrt{5}},\Delta =-\sqrt{5}{sin}^3\alpha <0,\end{array}$$

(38)

and the field (33) is of infinite energy. Its Fourier image is

$$\begin{array}{c}\mathcal{F}\left[{\mathrm{Ai}}_3\left(\mathit{\rho }\right|\mathbf{M})\right]\left(\mathbf{r}\right)=\frac{\sqrt{10+2\sqrt{5}}}{5^{2/3}2\pi }exp(\mathrm{i}{P}_1+\mathrm{i}{P}_3)\mathrm{Ai}(P_0+P_2),\\ \begin{array}{c}{P}_0=-\frac{1}{5^{1/6}}\left(C_1+C_2-\frac{1+\sqrt{5}}{2}{C}_3\right),\\ P_1=\sqrt{\frac{5+\sqrt{5}}{10}}(C_2-C_1)·x+\frac{(3+\sqrt{5})(C_1+C_2)-4C_3}{2\sqrt{5}}·y,\\ P_2=-\frac{1+\sqrt{5}}{5^{2/3}}(x^2+y^2),P_3=\frac{2(2+\sqrt{5})}{5}·x^{2}y+\frac{2(4+\sqrt{5})}{15}·y^3.\end{array}\end{array}$$

(39)

Since $P_2$ is negative, the field (39) has an infinite number of zero circles. Figure 9 depicts an example when $C_1=C_2=c$, $C_3=\frac{1}{2}(3+\sqrt{5})c$, and $5^{1/3}c=a_1^{\prime }$. Then, the polynomial $P_1$ vanishes, and the intensity becomes maximum at the origin.

(iii)

If $\alpha =3\pi /5$, then

$$\begin{array}{c}cos\alpha =-\frac{\sqrt{5}-1}{4},sin\alpha =\frac{\sqrt{5}-1}{4}\sqrt{5+2\sqrt{5}},\Delta =\sqrt{5}{sin}^3\alpha >0,\end{array}$$

(40)

and the field (33) is of finite energy. Its Fourier image is

$$\begin{array}{c}\mathcal{F}\left[{\mathrm{Ai}}_3\left(\mathit{\rho }\right|\mathbf{M})\right]\left(\mathbf{r}\right)=\frac{\sqrt{10-2\sqrt{5}}}{5^{2/3}2\pi }exp(\mathrm{i}{P}_1+\mathrm{i}{P}_3)\mathrm{Ai}(P_0+P_2),\\ \begin{array}{c}{P}_0=\frac{1}{5^{1/6}}\left(C_1+C_2+\frac{\sqrt{5}-1}{2}{C}_3\right),\\ P_1=\sqrt{\frac{5-\sqrt{5}}{10}}(C_2-C_1)·x-\frac{(3-\sqrt{5})(C_1+C_2)-4C_3}{2\sqrt{5}}·y,\\ P_2=\frac{\sqrt{5}-1}{5^{2/3}}(x^2+y^2),P_3=-\frac{2(\sqrt{5}-2)}{5}·x^{2}y+\frac{2(4-\sqrt{5})}{15}·y^3.\end{array}\end{array}$$

(41)

The polynomial $P_2$ is positive, so the field (41) has a finite number of zero circles. As before, we choose the coefficients $C_n$ in such a way as to vanish the polynomial $P_1$: $C_1=C_2=c$, $C_3=\frac{1}{2}(3-\sqrt{5})c$. When $5^{1/3}c=a_1^{\prime },a_1,a_4^{\prime },a_4$, the initial three-Airy beams, and their Fourier images are depicted in Figure 10.

Finally, let us consider the field (33) for the limiting case $\alpha \to +0$. Then $sin\alpha \cong \alpha$, $cos\alpha \cong 1$, $\Delta \cong -6{\alpha }^3$, and

$$\begin{array}{c}{P}_0=\frac{2C_3-C_1-C_2}{6^{1/3}},P_1=\frac{C_2-C_1}{2\alpha }·x+\frac{2(C_1+C_2)-C_3}{3}·y,\\ P_2=-\frac{1}{6^{4/3}}\left(\frac{3x^2}{{\alpha }^2}+4y^2\right),P_3=\frac{x^{2}y}{3{\alpha }^2}+\frac{5y^3}{27}.\end{array}$$

(42)

Substituting these expressions into Equation (28) and integrating both sides over $x\in \mathbb{R}$ (at his stage, the integral over $\xi$ reduces to Dirac’s delta function and $\alpha$ vanishes), we have

$$\begin{array}{c}\hfill 2\pi {\int }_{\mathbb{R}}{\mathrm{e}}^{-\mathrm{i}y\eta }\prod _{1\le n\le 3}\mathrm{Ai}(\eta +C_n)\mathrm{d}\eta =\frac{1}{6^{1/3}}exp\left[\frac{\mathrm{i}(2C_1+2C_2-C_3)y}{3}+\frac{5\mathrm{i}{y}^3}{27}\right]\\ \hfill \times {\int }_{\mathbb{R}}exp\left[\frac{\mathrm{i}(C_2-C_1)x}{2}+\frac{\mathrm{i}{x}^{2}y}{3}\right]\mathrm{Ai}\left(\frac{2C_3-C_1-C_2}{6^{1/3}}-\frac{3x^2+4y^2}{6^{4/3}}\right)\mathrm{d}x.\end{array}$$

(43)

This quite unusual identity is applicable in the theory of hypergeometric functions and helps to investigate the elliptic umbilic diffraction catastrophe [26], and has various corollaries. We provide only one of them, the integral representation of the product of both Airy functions (see also Section 3.6.3 in [36]):

$$\mathrm{Ai}\left(c\right)\mathrm{Bi}\left(c\right)=3^{1/6}{\int }_{\mathbb{R}}{\mathrm{Ai}}^2\left(\eta \right)\mathrm{Ai}(\eta +3^{1/3}c)\mathrm{d}\eta .$$

(44)

To prove this, we substitute $y=0$ into (43), then the integral over x is known due to Berry ([31], Equation (B17)). The formula (44) helps to study the propagation of certain Gaussian and non-Gaussian light fields of finite and infinite energy containing the factor $\mathrm{AiBi}$.

5. Conclusions

The problems studied above related to the propagation of three-Airy beams in the Fresnel zone have close connections with hot topics in optics. We believe that generalized three-Airy beams will be useful for various applications of modern optics and photonics, such as laser structuring of photosensitive material surfaces and laser ablation [38,39], optical trapping and manipulation [40,41].

For the theory, the results obtained can be used as a basis for studying other light fields. For example, when talking about Airy laser beams, we usually mean the beams constructed for the function $\mathrm{Ai}\left(x\right)$, ignoring the function $\mathrm{Bi}\left(x\right)$. The only exception known to the authors is Ref. [42], in which the study of the initial light field $F_0\left(x\right)=exp(-x^2/w^2)\mathrm{Bi}(Ax+B)$ was conducted by Torre. It is quite easy to transfer all the results obtained for Airy–Gaussian beams in [3] (here, Airy means Ai) to the case of Bi. One of the possible ways is to apply the fact that both functions, $\mathrm{Ai}\left(x\right)$ and $\mathrm{Bi}\left(x\right)$, are solutions of the differential equation, $y^{\prime \prime }=xy$. The second way is to use algebraic relations ([43], Equations (10.4.6) and (10.4.9)),

$$2\omega \mathrm{Ai}\left(\omega x\right)=-\mathrm{Ai}\left(x\right)+\mathrm{i}·\mathrm{Bi}\left(x\right),2\omega \mathrm{Bi}\left(\omega x\right)=-\mathrm{Bi}\left(x\right)+3\mathrm{i}·\mathrm{Ai}\left(x\right),$$

(45)

where $\omega ={\mathrm{e}}^{2\pi \mathrm{i}/3}$ is a cubic root of unity.

For the initial field $F_0\left(x\right)=exp(-x^2/w^2)\{\mathrm{Ai}+\mathrm{i}·\mathrm{Bi}\}(Ax+B)$, it is also easy to construct a solution of the paraxial equation. It can be seen that the resulting beam, $F(x,z)$, exhibits the property of “bending the trajectory” as it propagates.

In addition, the connection between Airy functions with Bessel functions of orders $\pm 1/3$ makes it possible to construct fractional Bessel–Gaussian (BG) beams of these orders. These beams are quite different from the BG beams of integer and half-integer orders known before [44,45,46]. It would be especially interesting to study three-Airy beams of infinite energy, i.e., fractional-order Bessel beams, to reveal the effect of larger depth of focus [47]. Moreover, replacing Airy functions in the initial fields above with the Pearcey and Swallowtail functions leads to light fields related to Bessel functions of the order of $\pm 1/4$ and $\pm 1/5$.

Finally, the study of the propagation of generalized three-Airy beams ${\mathrm{Ai}}_3\left(\mathbf{r}\right|\mathbf{M})$ in the Fresnel zone requires the use of diffraction catastrophe integrals of order higher than 3. For example, considering Equation (28) for the general case of $P_2$ instead of the radially symmetric case $P_2\left(\mathbf{r}\right)=c{r}^2$, one can obtain

$$\begin{array}{cc}\hfill {\int }_{\mathbb{R}}{\mathrm{e}}^{\mathrm{i}x\xi }\mathrm{Ai}\left(2^{2/3}\xi \right)\mathrm{Ai}(\xi +C_1)\mathrm{Ai}(\xi +C_2)\mathrm{d}\xi & =\frac{{\left|x\right|}^{-1/4}}{2^{5/3}{\pi }^{3/2}}exp\left(\frac{\mathrm{i}x{X}^2}{4}-\frac{\mathrm{i}{x}^3}{12}-\frac{\pi \mathrm{i}}{4}\mathrm{sgn}x\right)\hfill \\ \hfill & \times \left\{\begin{array}{cc}\mathrm{Pe}\left[(C_2-C_1){\left|x\right|}^{1/4},-X{\left|x\right|}^{1/4}\right],& (x>0),\\ {\mathrm{Pe}}^{*}\left[-(C_2-C_1){\left|x\right|}^{1/4},X{\left|x\right|}^{1/4}\right],& (x<0),\end{array}\right.\hfill \end{array}$$

(46)

where $X=x/2-(C_1+C_2)/x$, the asterisk means complex conjugation, and

$$\mathrm{Pe}(x,y)={\int }_{\mathbb{R}}exp\left(\mathrm{i}{t}^4+\mathrm{i}y{t}^2+\mathrm{i}xt\right)\mathrm{d}t$$

(47)

is the Pearcey function [48,49], i.e., the diffraction catastrophe of order 4. In particular,

$${\int }_{\mathbb{R}}{\mathrm{e}}^{\mathrm{i}x\xi }\mathrm{Ai}\left(2^{2/3}\xi \right){\mathrm{Ai}}^2\left(\xi +\frac{x^2}{4}\right)\mathrm{d}\xi =\frac{{\Gamma (5/4)\left|x\right|}^{-1/4}}{2^{2/3}{\pi }^{3/2}}exp\left(-\frac{\mathrm{i}{x}^3}{12}-\frac{\pi \mathrm{i}}{8}\mathrm{sgn}x\right),$$

(48)

since $\mathrm{Pe}(0,0)=2\Gamma (5/4){\mathrm{e}}^{\pi \mathrm{i}/8}$. This result confirms a very important statement given in [50]: Fourier and Fresnel transforms of light fields constructed using diffraction catastrophes of order n sometimes require the use of diffraction catastrophes whose order exceeds n.