## Abstract

Perturbing the external control parameters of nonlinear systems leads to dramatic changes in their bifurcations. A branch of singular theory, the catastrophe theory, analyzes the generating function that depends on state and control parameters. It predicts the formation of bifurcations as geometrically stable structures and categorizes them hierarchically. We apply the approach to evaluate the catastrophe diffraction integral with respect to two-dimensional cross sections through the control parameter space and thus transfer these bifurcations to optics, where they manifest as caustics in transverse light fields. For all optical catastrophes that depend on a single state parameter (cuspoids), we analytically derive a universal expression for the propagation of all corresponding cuspoid beams. We show that the dynamics of the resulting cuspoids can be expressed by higher-order optical catastrophes with dynamically changing control parameters. We show analytically and experimentally that particular swallowtail beams develop caustics with geometrical structures corresponding to higher-order butterfly catastrophes during propagation, whereas differently tailored swallowtail beams decay to lower-order cusp catastrophes.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. INTRODUCTION

Singularities of the generating function of a nonlinear system describe how the structures of its bifurcations are related to external control parameters. Small perturbations of the control parameters may result in dramatic changes. This is what is called a catastrophe in the context of the singularity and the catastrophe theory [1,2]. A classification of the most fundamental potentials exists, which distinguishes different local effects near catastrophes and the system’s qualitative behavior [1,2]. Each order of a catastrophe forms a specific geometrically stable structure in the control parameter space. According to the unique geometric structure, a hierarchical categorization with increasing dimensionality of the control parameter space is given by fold, cusp, swallowtail, hyperbolic umbilic, elliptic umbilic, butterfly, and parabolic umbilic catastrophes.

In general, any catastrophe has a wide range of different areas of physics that it explains [3–6], where individual catastrophes are present for miscellaneous systems. Exemplary, the fold catastrophe manifests in rainbows [7], the cusp catastrophe emerges in models of social sciences [8], and the swallowtail catastrophe occurs in orbits of hydrogen [9], to name only some exemplary orders of catastrophes.

Manifestation of catastrophes in light phenomena as they occur in nature, where high-intensity caustics emerge in defined geometries according to their respective order, is well understood [10,11]. Caustics are formed due to reflections on smoothly curved surfaces and emerge, e.g., as a cusp in a cup, or arise due to refraction at spatially modulated boundaries, e.g., when forming ramified networks at the ground of shallow waters. Beneath their potential to optically visualize the complex dynamics of nonlinear systems, it is their sharp high-intensity boundary with unique propagation paths on curved trajectories [12–14], which makes caustics in light highly attractive for a broad range of applications [14–17]. However, tailored creation of caustics in light beyond their natural occurrence is necessary to control their properties and dynamics in all facets.

Over the past years, application of spatial light modulators has facilitated the artificial creation of catastrophes in light fields; however, until now only the most fundamental orders were realized as paraxial Airy (fold) and Pearcey (cusp) beams [12,13]. Artificial higher-order optical swallowtail and butterfly catastrophes were transferred to caustics in transverse light fields by mapping the cross sections of the higher-dimensional control parameter space to the two-dimensional (2D) transverse plane [18]. The approach presented allows realizing fundamentally different geometrical caustic structures in the initial transverse plane without propagating the light fields by choosing corresponding cross sections in the control parameter space.

Only recently was a study proposed addressing the propagation properties of Airy and Pearcey caustics, as well as a specific swallowtail caustic, giving a first glance into the dynamics of higher-order optical catastrophes [19]. Their general dynamics, however, are unknown. Of particular interest is the stability of higher-order caustic structures during propagation. They are expected to decay to lower-order caustics, as it is typical for high-dimensional singularities, since the initial caustic structures represent mappings of the related higher-order catastrophes to the transverse plane. Furthermore, the propagation of higher-order caustic beams is expected to provide similar breakthrough potential to their already well-known mates of lower orders.

Thus, with this work we contribute to the analysis of the dynamics of optical swallowtail beams, a particular solution of the paraxial Helmholtz equation. We analytically calculate and prove experimentally how higher-order swallowtail beams evolve in space. We show that the dynamics of swallowtail beams can again be described in terms of swallowtail beams or contain fingerprints of optical butterfly catastrophes, where the control parameters depend on $z$. For a different set of parameters, we demonstrate how the corresponding initial swallowtail caustic decays into a lower-order cusp during propagation.

## 2. CALCULATING THE PROPAGATION OF OPTICAL CATASTROPHES

Caustics emerging as optical catastrophes that depend on a single state parameter $s$ are from co-rank 1 and are called cuspoids [20]. We consider the complete class of optical cuspoids representing a particular solution of the paraxial Helmholtz equation and emerging as a solution of the diffraction catastrophe integral ${C}_{n}(\mathbf{a})$ [11,20]:

All resulting light structures ${C}_{n}(\mathbf{a})$ exhibit individual caustic profiles, which are characterized by degenerate critical points of the potential function ${P}_{n}(\mathbf{a},s)$. Here, the vector $\mathbf{a}$ consists of all dimensionless control parameters ${a}_{j}$, with $j=1,\dots ,n-2$, and spans the *control parameter space* with co-dimension $n-2$. We identify two of the $n-2$ control parameters with dimensionless transverse spatial coordinates $X$, $Y$ [11] and keep the remaining $n-4$ control parameters constant. Thereby, the real space transverse coordinates $X=x/{x}_{0},Y=y/{y}_{0}$ are normalized by introducing characteristic structure sizes ${x}_{0},{y}_{0}$. We restrict the following discussion to optical cuspoids with $n\ge 4$.

We emphasize that the Airy beam $Ai(X)={C}_{3}(\mathbf{a})$ results as a one-dimensional (1D) structure identified with the spatial coordinate $\mathbf{a}={a}_{1}=X$. Thus, our approach is also valid for $n=3$ by reducing all transverse 2D considerations to 1D. However, the Pearcey beam $\mathrm{Pe}(X,Y)={C}_{4}(\mathbf{a})$ leads *per se* to a 2D distribution, by identifying $\mathbf{a}={({a}_{1},{a}_{2})}^{T}={(X,Y)}^{T}$. Consequently, for higher-order optical cuspoids, such as, e.g., the swallowtail beams $\mathrm{Sw}(\mathbf{a})={C}_{5}(\mathbf{a})$, we have to choose a single control parameter from three existing ones $\mathbf{a}={({a}_{1},{a}_{2},{a}_{3})}^{T}$ to be constant in order to map the characteristics of the corresponding catastrophe to the 2D transverse light field by identifying the two remaining control parameters with spatial coordinates ${(X,Y)}^{T}$. Thus, for the optical swallowtail catastrophe, three generic swallowtail beams arise as orthogonal cross sections through the control parameter space, and, correspondingly, more for the butterfly beam $Bu(\mathbf{a})={C}_{6}(\mathbf{a})$, which maps the characteristics of the butterfly catastrophe. This approach is thoroughly presented in [18].

Since caustics represent geometrically stable structures determined by the singular mapping of the potential function ${P}_{n}$ on the $n-2$ dimensional plane $\mathbf{a}$, we calculate degenerate fixed points ${s}_{i}$ as where the first and second derivatives vanish. The submanifold ${s}_{i}$ then subsequently determines the geometric structure of the catastrophe:

By transferring these potentials to optics via Eq. (1), their singularities correspond to geometrically stable catastrophes and manifest as caustics in paraxial wave structures.

Starting with this formalism to create the class of paraxial cuspoid beams including the well-known Airy and Pearcey beams, as well as the less known swallowtail and butterfly beams, we derive analytical expressions for the propagation of these intriguing light structures and, furthermore, analyze the $z$-dependent evolution of their caustics in the regime of paraxial light by referring to Eq. (2).

In order to analytically calculate the propagation of a scalar transverse light field, we use the Huygens–Fresnel integral [21] and derive $z$-dependent expressions for cuspoid beams [Eq. (1)] in the paraxial regime. Especially since the Schrödinger equation describing the dynamics of quantum particles shows formal similarities to the paraxial Helmholtz equation providing an expression for the propagation of electromagnetic waves, the equivalence of time $t$ and longitudinal propagation distance $z$ in suited systems is present [22], which justifies discussing the propagation of beams in terms of their dynamics where applicable.

Our approach links two control parameters ${a}_{\alpha},{a}_{\beta}$ with the transverse coordinates $X$, $Y$, where $\alpha ,\beta \in j=1,\dots ,n-2$, and treats the remaining control parameters ${\mathbf{a}}^{\prime}$ as constant. Consequently, the $z$ dependence of the cuspoid beams is derived to be

Here, $k=2\pi /\lambda =\sqrt{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}$ is the wave number and $\lambda $ is the wavelength. Similarly to the Rayleigh length of Gaussian beams, we can define characteristic Rayleigh lengths as ${z}_{ex}=2k{x}_{0}^{2}$ and ${z}_{ey}=2k{y}_{0}^{2}$. They play an important role in the dynamics of the beams and depend on the wavelength $\lambda =2\pi /k$ and the corresponding structure sizes ${x}_{0},{y}_{0}$.

Using Eq. (3) with proper choice of parameters, one can easily derive the well-known propagation expressions for the Airy [23] and the Pearcey beam [13].

Equation (3) reveals fundamentally new insights in the dynamics of caustics in light: in the case of Airy and Pearcey beams, their propagation is again expressed by Airy and Pearcey functions (displacement and form-invariant scaling, respectively) [12,13]. Starting with the order $n\ge 5$ with co-dimension $n-2$, the propagation of the corresponding cuspoid beam can be described in terms of cuspoid beams up to order $2(n-2)$, if identifying appropriate control parameters with transverse spatial coordinates. This relation is manifested in the propagator ${T}_{\alpha ,\beta}(s,z)=\mathrm{exp}[-i\frac{z}{{z}_{ex}}{s}^{2\alpha}]\mathrm{exp}[-i\frac{z}{{z}_{ey}}{s}^{2\beta}]$ of Eq. (3), where one of the exponents $2\alpha $ or $2\beta $ can always be chosen to be larger than the leading exponent $n$ of the potential function ${P}_{n}$, and contributes only if $z\ne 0$. Thus, by choosing $\alpha $ or $\beta $ to be the coefficient, which corresponds to the state variable of the highest degree (i.e., $n-2$), the propagation of this $n$th-order cuspoid beam is given by a static cuspoid beam with a $2(n-2)$-dimensional control parameter space and $z$-dependent control parameters.

In the following, we exemplary demonstrate our concept by applying Eq. (3) for two initial transverse swallowtail beams to show their fundamentally different dynamics. First, we calculate and experimentally obtain the dynamics of a $\mathrm{Sw}(X,Y,{a}_{3})$ beam. We show that its propagation can be described by an optical catastrophe of the same order with varying control parameters. Second, analytical and experimental investigation of a $\mathrm{Sw}(X,{a}_{2},Y)$ beam reveals that its dynamics are represented by a higher-order optical butterfly catastrophe whose control parameters are functions of $X,Y,z,{a}_{2}$.

## 3. EXPERIMENTAL DETAILS

In order to obtain the spatial intensity distribution experimentally, we have used the setup shown in Fig. 1. As the light source serves a frequency-doubled $\mathrm{Nd}\text{:}{\mathrm{YVO}}_{4}$ laser with a wavelength of $\lambda =532\text{\hspace{0.17em}}\mathrm{nm}$. The linearly polarized beam is expanded and collimated. As a plane wave, it illuminates a HOLOEYE HEO 1080P reflective LCOS phase-only spatial light modulator with full-HD resolution. We modulate both the amplitude and phase of the beam by encoding both information in a phase-only pattern and applying an appropriate Fourier filter [24]. The beam is then imaged by a camera after passing through a microscope objective. Both the camera and the microscope objective are mounted on a $z$-direction movable stage capable of scanning the propagation of the light fields by obtaining transverse intensity patterns.

## 4. DYNAMICS OF THE $\mathsf{Sw}\mathbf{(}\mathsf{X}\mathsf{,}\mathsf{Y}\mathsf{,}{\mathsf{a}}_{\mathsf{3}}\mathbf{)}$ BEAM

As pointed out in our analytic description, particular emphasis will be on the comparison between the dynamics of a $\mathrm{Sw}(X,Y,{a}_{3})$ beam with those of the $\mathrm{Sw}(X,{a}_{2},Y)$ or $\mathrm{Sw}({a}_{1},X,Y)$ beams, since the propagation of the $\mathrm{Sw}(X,Y,{a}_{3})$ beam is a mapping from the three-dimensional (3D) control parameter space onto itself, whereas the evolution of the two other beams is described in the higher four-dimensional (4D) control parameter space. In order to investigate these dynamics, this section substantiates our analytical investigations with experimental realizations of the $\mathrm{Sw}(X,Y,{a}_{3})$ beam and the arbitrarily chosen $\mathrm{Sw}(X,{a}_{2},Y)$ beam. The dynamics of the $\mathrm{Sw}({a}_{1},X,Y)$ beam are similar to that of the $\mathrm{Sw}(X,{a}_{2},Y)$ beam, as stated in Eq. (3), and are not shown here explicitly. Its dynamics can be deduced easily by considering Eq. (3).

For the $\mathrm{Sw}(X,Y,{a}_{3})$ beam with $n=5$, we chose $\alpha =1$ and $\beta =2$. Equation (3) then leads to a noncanonical fifth-order potential function in the exponential, which can be brought into canonical form by substituting $s\to u+z/5{z}_{ey}$ to suppress the next-to-leading order term. The dynamics of the $\mathrm{Sw}(X,Y,{a}_{3},Z)$ beam, where $Z=z/{z}_{ex}$ is the normalized longitudinal coordinate, can then be described in terms of again an optical swallowtail catastrophe $\mathrm{Sw}({C}_{1},{C}_{2},{C}_{3})$, whose control parameters are functions of $X,Y,{a}_{3},Z$ in the form of

wherePlotting Eq. (4) for different $z$ positions allows visualizing the dynamics of a $\mathrm{Sw}(X,Y,{a}_{3})$ beam. We set ${a}_{3}=0$. The analytically calculated propagation is depicted at the top in Fig. 2. Furthermore, we experimentally obtain the volume intensity distribution of the beam for different $z$ positions and have shown them at the bottom in Fig. 2. We chose transverse feature sizes of ${x}_{0}={y}_{0}=8\text{\hspace{0.17em}}\mathrm{\mu m}$ and a propagation distance of 20 mm, which covers the most interesting effects in the sketched intensity volume.

Analysis of Eq. (3) reveals that at $z=0\text{\hspace{0.17em}}\mathrm{mm}$ the $\mathrm{Sw}(X,Y,{a}_{3})$ beam shows mirror symmetry with respect to $Y$:

where the asterisk denotes a complex conjugation. Due to the symmetry, the beam’s propagation for negative $z$ values is not shown, but has also been observed numerically and experimentally.A unique feature of the $\mathrm{Sw}(X,Y,{a}_{3})$ beam becomes apparent: during propagation, the transverse field at the origin ($z=0\text{\hspace{0.17em}}\mathrm{mm}$) turns out to consist of a *fast diffracting contribution*, which quickly vanishes due to its transverse momentum (e.g., visible at $z\approx 10\text{\hspace{0.17em}}\mathrm{mm}$), subsequently revealing a field contribution (remaining field at $z\approx 20\text{\hspace{0.17em}}\mathrm{mm}$) that resembles a Pearcey beam in the field distribution as well as in the propagation.

This similarity is additionally manifested in the distribution of Fourier components since they are located on a parabola for both the Pearcey and $\mathrm{Sw}(X,Y,{a}_{3})$ beams in the form of $\tilde{\mathrm{P}}\mathrm{e}({k}_{x},{k}_{y})=\delta ({({x}_{0}{k}_{x})}^{2}-{y}_{0}{k}_{y})\mathrm{exp}[i{({x}_{0}{k}_{x})}^{4}]$ and $\tilde{\mathrm{S}}\mathrm{w}({k}_{x},{k}_{y},0)=\delta ({({x}_{0}{k}_{x})}^{2}-{y}_{0}{k}_{y})\mathrm{exp}[i{({x}_{0}{k}_{x})}^{5}]$, respectively [13,18]:

Due to this spectral similarity, we express the $\mathrm{Sw}(X,Y,{a}_{3}=0)$ beam as a convolution of a Pearcey beam with a dimensionally reduced swallowtail beam:

Here, $\delta (\dots )$ is the delta function. Note that the initial beam profile of the artificially designed $\mathrm{Sw}(X,Y,{a}_{3})$ beam shows a swallowtail caustic as the mapping of a cross section through the higher-dimensional 3D control parameter space to the lower-dimensional 2D transverse plane. Due to this, the swallowtail structure in Fig. 2 becomes unstable during propagation with respect to the otherwise geometrically stable swallowtail catastrophe in nonlinear systems, and decays to a lower-dimensional optical cusp catastrophe. We will discuss the dynamics of the caustic in more detail in Section 6.

## 5. DYNAMICS OF THE $\mathsf{Sw}\mathbf{(}\mathsf{X}\mathbf{,}{\mathsf{a}}_{\mathbf{2}}\mathbf{,}\mathsf{Y}\mathbf{)}$ BEAM EXPRESSED IN TERMS OF A STATIC BUTTERFLY BEAM

The dynamics of diverse optical swallowtail beams and their caustics are fundamentally different if compared with each other. In order to demonstrate that caustic beams of a certain order $n$ may show structural similarities with higher-order caustics during propagation, we exemplary investigate a $\mathrm{Sw}(X,{a}_{2},Y)$ beam, and, thus, set ${a}_{2}=\text{const.}$ and choose $\alpha =1$, $\beta =3$. Applying Eq. (3) with these parameters leads to a butterfly integral. Its reduction to the canonical butterfly form requires appropriate substitutions for two distinct cases $z<0$ or $z>0$.

Starting with $z<0$, we substitute $s\to u\xb7{(-z/{z}_{ey})}^{(-1/6)}+{z}_{ey}/6z$ and express the propagation of the $\mathrm{Sw}(X,{a}_{2},Y,Z)$ beam in terms of the higher-order static butterfly beam $\mathrm{Bu}({B}_{1},{B}_{2},{B}_{3},{B}_{4})$ with an overall phase factor $\psi $ in dependence of new control parameters $\mathbf{B}=\{{B}_{1},{B}_{2},{B}_{3},{B}_{4}\}$ as functions of $X,Y,Z,{a}_{2}$:

For $z>0$, an overall conjugation procedure is required, subsequently followed by the substitution $s\to u\xb7{(z/{z}_{ey})}^{(-1/6)}+{z}_{ey}/6z$. Defining $q={z}_{ey}/(6{z}_{ex})$ and $\gamma ={(Z/(6q))}^{(1/6)}$ for this case, one can obtain

Equation (9) is plotted for different $z$ values and illustrated at the top in Fig. 3. We obtained the dynamics as well in the experiment, which is shown below. The transverse dimensions are ${x}_{0}={y}_{0}=8\text{\hspace{0.17em}}\mathrm{\mu m}$ and the propagation distance is 20 mm. Due to the point symmetry of the $\mathrm{Sw}(X,{a}_{2},Y)$ beam [cf. Eq. (3)],

and for illustrative reasons, we restrict the evaluation to $z\ge 0$.The plotted analytical description of the propagation of a $\mathrm{Sw}(X,{a}_{2},Y)$ beam is in high agreement with the corresponding experimental realization. The transverse plane at $z=0\text{\hspace{0.17em}}\mathrm{mm}$ shows no mirror or point symmetry, as shown in Fig. 3, and can be proved by considering Eq. (1). However, during propagation the structure separates into two individual structures, showing point symmetry. This suggests that we can assume the $\mathrm{Sw}(X,{a}_{2},Y)$ beam at the origin to consist of two major contributions: one fast diffracting part, which breaks the symmetry in the initial plane, and a second one, which remains during propagation and exhibits point symmetry. After a certain propagation distance ($z\approx 10\text{\hspace{0.17em}}\mathrm{mm}$), the emerging point symmetric structure can easily be considered as a doubled Pearcey-beam-like structure. Nevertheless, the contribution rather occurs due to the similarity of the $\mathrm{Sw}(X,{a}_{2},Y)$ beam to a $\mathrm{Bu}(X,{a}_{2},Y,{a}_{4})$ beam, whose $z=0\text{\hspace{0.17em}}\mathrm{mm}$ real space appearance [18] resembles the field distribution of the $\mathrm{Sw}(X,{a}_{2},Y)$ beam after a certain propagation distance. Their spectra both correspond to cubic functions, ${y}_{0}{k}_{y}={({x}_{0}{k}_{x})}^{3}$, and, thus,

This leads to the real space expression

## 6. DYNAMICS OF THE SWALLOWTAIL CAUSTIC: DECAY TO LOWER-ORDER CUSP

The propagation of the $\mathrm{Sw}(X,Y,{a}_{3})$ beam shown in Fig. 2 gives rise to further study of the dynamics of the corresponding caustic. By parameterizing Eq. (4) and analyzing the structure of the caustics according to Eq. (2), we find the $z$-dependent dynamics of the caustic in the transverse plane.

Figure 4 illustrates the parameterized surface of the caustic in real space ($x$, $y$, $z$). The plane at $z=0$, where the green line highlights the form of the caustic, equals the front plane in Fig. 2.

The surface appears to consist of two intertwined slices, merging point symmetrically at the origin. Each of them has a cuspoid form, which changes with $z$ position, becomes flattened, and is strongly bent. This can be recognized as the *fast diffracting contribution* described above. The remaining part is the cusp of the other slice. According to Eq. (5), we find that the previously discussed properties of the beam resemble the dynamics of the caustic. Since the lower-dimensional initial light field contains the fingerprint of the optical swallowtail caustic due to the mapping of a cross section from the higher-dimensional control parameter space, the expected decay of the swallowtail catastrophe to a lower-order cusp is apparent. Although the structural stability of the swallowtail catastrophe with a 3D control parameter space is lost in 2D, the demonstrated swallowtail light fields show caustic structures, which exhibit novel and unique propagation properties like, e.g., high intensities on curved paths.

## 7. CONCLUSION

Perturbing the external control parameters of a nonlinear system described by a potential function leads to so-called catastrophes at locations where bifurcations suddenly shift [1,2,25]. These singularities of the gradient maps of potential functions manifest as caustics in light, and were studied extensively as natural phenomena in the late seventies and eighties [3,10,11,26,27]. Mapping catastrophes to light fields via the paraxial catastrophe integral of Eq. (1) is a well-known approach. However, it took until 2007 to use spatial light modulators for creating Airy beams [12]. This allows for the first time controlled mapping of caustics to paraxial light, and starts the renaissance of designing artificial caustic beams.

In our contribution, we considered the complete class of cuspoid beams, and, thus, optical catastrophes depending on a single state parameter $s$, by mapping higher-order catastrophes to the lower-dimensional initial transverse plane of paraxial beams and analytically derived a general equation for their propagation. Our approach connects two of the control parameters, $\mathbf{a}$ with the transverse spatial coordinates ${(x,y)}^{T}$, and, thus, the cuspoid beams are realized as cross sections through the higher-dimensional control parameter space. We showed that, depending on the control parameters identified with the spatial coordinates, the propagation of cuspoid beams of order $n$ can be calculated in terms of higher-order static cuspoid beams of order $2(n-2)$.

We demonstrated this by analytically calculating and experimentally obtaining the dynamics of the $\mathrm{Sw}(X,Y,{a}_{3})$ and $\mathrm{Sw}(X,{a}_{2},Y)$ beams. We proved that the evolution of the latter beam is linked to a static butterfly beam. Furthermore, we analyzed the dynamics of the swallowtail caustic and showed how its surface evolves in the 3D real space. We thereby demonstrated that the swallowtail catastrophe at an initial plane continuously decays to a lower-order cusp catastrophe during propagation.

The demonstrated optical catastrophes are highly attractive for microscopy and super-resolution applications. The propagation of the high-intensity rims near caustics capable of forming tailored structures is unique for each order of catastrophe and parameter set, and pave the way to advanced micromachining on tailored curves [28] and the realization of waveguides with a rich diversity of light-guiding paths [14].

## Funding

Westfälische Wilhelms-Universität Münster (Open Access Publication Fund).

## REFERENCES

**1. **R. Thom, *Structural Stability and Morphogenesis* (W. A. Benjamin, 1972).

**2. **V. I. Arnol’d and G. S. Wassermann, *Catastrophe Theory*, 3rd ed. (Springer, 2003).

**3. **M. V. Berry, “Waves and Thom’s theorem,” Adv. Phys. **25**, 1–26 (1976). [CrossRef]

**4. **D. Gifford, C. Miller, and N. Kern, “A systematic analysis of caustic methods for galaxy cluster masses,” Astrophys. J. **773**, 116 (2013). [CrossRef]

**5. **R. Höhmann, U. Kuhl, H.-J. Stöckmann, L. Kaplan, and E. J. Heller, “Freak waves in the linear regime: a microwave study,” Phys. Rev. Lett. **104**, 093901 (2010). [CrossRef]

**6. **A. Mathis, L. Froehly, S. Toenger, F. Dias, G. Genty, and J. M. Dudley, “Caustics and rogue waves in an optical sea,” Sci. Rep. **5**, 12822 (2015). [CrossRef]

**7. **H. Trinkhaus and F. Drepper, “On the analysis of diffraction catastrophes,” J. Phys. A **10**, L11–L16 (1977). [CrossRef]

**8. **J. E. Sheridan and M. A. Abelson, “Cusp catastrophe model of employee turnover,” Acad. Manage. J. **26**, 418–436 (1983). [CrossRef]

**9. **R. C. Méndez and F. F. Alfaro, “Fold, cusp, and swallowtail catastrophes in the evolution of planar closed orbits of hydrogen in crossed electric and magnetic fields,” Revista Colombiana de Fisica **45**, 28–44 (2013).

**10. **M. V. Berry, J. F. Nye, and F. J. Wright, “The elliptic umbilic diffraction catastrophe,” Philos. Trans. R. Soc. London A **291**, 453–484 (1979). [CrossRef]

**11. **M. V. Berry and C. Upstill, “Catastrophe optics: morphologies of caustics and their diffraction patterns,” Prog. Opt. **28**, 257–346 (1980). [CrossRef]

**12. **G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating Airy beams,” Phys. Rev. Lett. **99**, 213901 (2007). [CrossRef]

**13. **J. D. Ring, J. Lindberg, A. Mourka, M. Mazilu, K. Dholakia, and M. R. Dennis, “Auto-focusing and self-healing of Pearcey beams,” Opt. Express **20**, 18955–18966 (2012). [CrossRef]

**14. **A. Zannotti, M. Rüschenbaum, and C. Denz, “Pearcey solitons in curved nonlinear photonic caustic lattices,” J. Opt. **19**, 094001 (2017). [CrossRef]

**15. **F. Diebel, B. M. Bokic, D. V. Timotijevic, D. M. Jovic Savic, and C. Denz, “Soliton formation by decelerating interacting Airy beams,” Opt. Express **23**, 24351–24361 (2015). [CrossRef]

**16. **J. Baumgartl, M. Mazilu, and K. Dholakia, “Optically mediated particle clearing using Airy wavepackets,” Nat. Photonics **2**, 675–678 (2008). [CrossRef]

**17. **M. Manousidaki, D. G. Papazoglou, M. Farsari, and S. Tzortzakis, “Abruptly autofocusing beams enable advanced multiscale photo-polymerization,” Optica **3**, 525–530 (2016). [CrossRef]

**18. **A. Zannotti, F. Diebel, M. Boguslawski, and C. Denz, “Optical catastrophes of the swallowtail and butterfly beams,” New J. Phys. **19**, 053004 (2017). [CrossRef]

**19. **M. V. Berry, “Stable and unstable Airy-related caustics and beams,” J. Opt. **19**, 055601 (2017). [CrossRef]

**20. **J. F. Nye, *Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations* (IOP Publishing, 1999).

**21. **A. E. Siegman, *Lasers* (University Science Books, 1990).

**22. **A. Szameit and S. Nolte, “Discrete optics in femtosecond-laser-written photonic structures,” J. Phys. B **43**, 163001 (2010). [CrossRef]

**23. **M. V. Berry and N. L. Balazs, “Nonspreading wave packets,” Am. J. Phys. **47**, 264–267 (1979). [CrossRef]

**24. **J. A. Davis, D. M. Cottrell, J. Campos, M. J. Yzuel, and I. Moreno, “Encoding amplitude information onto phase-only filters,” Appl. Opt. **38**, 5004–5013 (1999). [CrossRef]

**25. **E. C. Zeeman, “Catastrophe theory,” in *Structural Stability in Physics* (Springer, 1979), pp. 12–22.

**26. **M. V. Berry, “Cusped rainbows and incoherence effects in the rippling-mirror model for particle scattering from surfaces,” J. Phys. A **8**, 566–584 (1975). [CrossRef]

**27. **Y. A. Kravtsov and Y. I. Orlov, “Caustics, catastrophes, and wave fields,” Usp. Fiz. Nauk **141**, 591–627 (1983). [CrossRef]

**28. **A. Mathis, F. Courvoisier, L. Froehly, L. Furfaro, M. Jacquot, P. A. Lacourt, and J. M. Dudley, “Micromachining along a curve: femtosecond laser micromachining of curved profiles in diamond and silicon using accelerating beams,” Appl. Phys. Lett. **101**, 071110 (2012). [CrossRef]