Dynamics of the optical swallowtail catastrophe

Perturbing the external control parameters of nonlinear systems leads to dramatic changes of its bifurcations. A branch of singular theory, the catastrophe theory, analyses 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 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, we analytically derive a universal expression for the propagation of all corresponding caustic beams. We reveal that the dynamics of the resulting caustics can be expressed by higher-order optical catastrophes. We show analytically and experimentally that particular swallowtail beams dynamically transform to higher-order butterfly caustics, whereas other swallowtail beams decay to lower-order cusp catastrophes.


Introduction
Singularities of the generating function of a nonlinear system describe how the structures of its bifurcations are related to external control parameters. Small per-turbations of the control parameters may result in dramatic changes. This is what is called a catastrophe in the context of singularity and catastrophe theory [1,2]. A classification of the most fundamental potentials distinguishes different local effects near catastrophes and defines the system's qualitative behaviour [1,2]. Each order of a catastrophe forms a specific geometrically stable structure in control parameter space. According to the unique geometric structure, a hierarchical categorization with increasing dimensionality of the control parameter space is given by the fold, cusp, swallowtail, hyperbolic umbilic, elliptic umbilic, butterfly, and parabolic umbilic catastrophe.
In general, any catastrophe has a wide range of different areas of physics that it explains [3,4,5,6], where individual catastrophes stand for miscellaneous systems. Exemplary, in meteorology the fold catastrophe manifests in rainbows [7], the cusp catastrophe emerges in models of social sciences describing an employee turnover [8], and the swallowtail catastrophe occurs in orbits of hydrogen in atom physics [9], to name only some exemplary orders of catastrophes.
The manifestation of catastrophes in natural light phenomena is given by high-intensity caustics that appear in defined geometries and different orders [10,11]. Caustics form 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 highintensity boundary with unique propagation paths on curved trajectories [12,13,14] that makes caustics in light highly attractive for a broad range of applications [15,16,17,14]. However, the tailored creation of caustics in light beyond their natural occurrence is necessary to control their properties and dynamics in all facets.
Only recently, the application of spatial light modulators has facilitated the artificial creation of catastrophes in light fields, but 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 cross-sections of the higher-dimensional control parameter space to the two-dimensional (2D) transverse plane [18]. The presented approach 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. Their 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 higherorder caustic beams is expected to provide similar breakthrough potential as their already well-known mates of lower order.
Thus, with this work we provide for the first time to our knowledge the general analysis of the propagation of optical catastrophes mapped to transverse light fields by evaluating the catastrophe diffraction integral for a potential that depends on a single state parameter. We prove our universal findings by demonstrating optical swallowtail beams, a particular solution of the paraxial Helmholtz equation. We analytically calculate and prove experimentally how the higher-order swallowtail beams evolve in space. We show that the dynamics of the swallowtail beams can again be described in terms of swallowtail beams or contain fingerprints of the 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.

The propagation of optical catastrophes
We consider the complete class of caustic beams emerging as optical catastrophes that depend on a single state parameter s. It represents a particular solution of the paraxial Helmholtz equation and emerges as solution of the diffraction catastrophe integral C n (a) [11,19] C n (a) = R e iPn(a,s) ds, P n (a, s) = s n + n−2 j=1 a j a 0j s j .
(1) All resulting light structures C n (a) exhibit individual caustic profiles that are characterized by degenerate critical points of the potential function P n (a, s).
Here, the vector a consists of all control parameters a j with j = 1, ..., n − 2 and spans the control parameter space with co-dimension n − 2. We identify two of the n − 2 control parameters with transverse spatial coordinates [11] and keep the remaining n − 4 control parameters constant. Dimensionless parameters are ensured by introducing characteristic structure sizes a 0j . We restrict the following discussion to optical catastrophes with n ≥ 4.
We emphasize that the Airy beam Ai(x) = C 3 (a) results as 1D structure identified with the spatial coordinate 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 Pe(x, y) = C 4 (a) leads per se to a 2D distribution, by identifying a = (a 1 , a 2 ) T = (x, y) T . Consequently, for higher-order optical catastrophes as e.g. the swallowtail beams Sw(a) = C 5 (a), we have to chose a single control parameter out of three existing ones a = (a 1 , a 2 , a 3 ) T to be constant in order to map 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(a) = C 6 (a) that maps 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 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 (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 caustic 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 (2).
In order to analytically calculate the propagation of a scalar transverse light field, we apply the angular spectrum method [20] and derive z-dependent expressions for caustic beams (1) in the paraxial regime. That is, the propagation of a known transverse electric field distribution C n (a) can be determined by first calculating the 2D Fourier transform C (k α , k β , a ) with respect to two control parameters a α , a β identified with x, y, where α, β ∈ j = 1, ..., n−2 and k α , k β are corresponding Fourier frequencies. We denote the remaining (constant) control parameters with a . The general Fourier transform of (1) can be found in the supplementary material (SM).
is the wave number and λ the wavelength. Consequently, the z-dependence of the caustic beams is derived to be Similar to the Rayleigh length of Gaussian beams, we can define characteristic Rayleigh lengths z ex = 2kx 2 0 and z ey = 2ky 2 0 that play an important role for the dynamics of the beams and depend on the wavelength λ = 2π/k and the corresponding the structure sizes x 0 , y 0 . In general, z e = (z eα , z eβ ) T = 2k(a 2 0α , a 2 0β ) T . Using Eq. (3) with the proper choice of parameters, one can easily derive the well known propagation expressions for the Airy [21] and Pearcey beam [13], respectively.
Eq. (3) reveals fundamentally new insights in the dynamics of caustics in light: In the case of the Airy and the Pearcey beams, their propagation is again expressed by Airy and Pearcey functions (displacement or form-invariant scaling, respectively) [12,13]. Starting with the order n ≥ 5 with co-dimension n − 2, the propagation of the corresponding caustic beam can be described in terms of caustic beams up to order 2(n − 2), if identifying the appropriate control parameters with transverse spatial coordinates. This relation is manifested in the propa- , where one of the exponents 2α or 2β can always be chosen to be larger then the leading exponent n of the potential function P n , and contributes only if z = 0. Thus, by choosing α or β to be the coefficient that corresponds to the state variable of highest degree (i.e. n − 2), the propagation of this n th -order caustic beam is given by a static caustic beam with 2(n − 2)-dimensional control parameter space and zdepending control parameters.
In the following, we exemplary demonstrate our concept by applying (3) for two initial transverse swallowtail beams to show their fundamentally different dynamics. First, we calculate and experimentally obtain the dynamics of a 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, the analytical and experimental investigation of a 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 .
In order to obtain the spatial intensity distribution experimentally we have used the setup shown in Fig. 1. As light source serves a frequency-doubled Nd:YVO 4 laser with a wavelength of λ = 532 nm. The linearly polarized beam is expanded and collimated. As plane wave, it illuminates a HOLOEYE HEO 1080P reflective LCOS phase-only spatial light modulator (SLM) with full HD resolution. We modulate both, amplitude and phase of the beam by encoding both informations in a phase-only pattern and applying an appropriate Fourier filter [22]. The beam is then imaged by a camera after passing a microscope objective. Both, camera and microscope objective are mounted on a in z-direction movable stage capable to scan the propagation of the light fields by obtaining transverse intensity patterns.

Dynamics of a swallowtail beam
As pointed out in our analytic description, particular emphasis will be on the difference between the dynamics of a Sw (x, y, a 3 ) beam in comparison to these of the Sw (x, a 2 , y) or Sw (a 1 , x, y) beams, since the propagation of the 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 fourdimensional (4D) control parameter space. In order to investigate these dynamics, this section substantiates our analytical investigations by experimental realizations of the Sw (x, y, a 3 ) beam and, arbitrarily chosen, the Sw (x, a 2 , y) beam. The dynamics of the Sw (a 1 , x, y) beam are similar to those of the Sw (x, a 2 , y) beam, as stated in (3), and are not shown here explicitly. The analytically calculated propagation is given in the SM. For the Sw (x, y, a 3 ) beam with n = 5 we chose α = 1 and β = 2. (3) then leads to a noncanonical 5 th -order potential function in the exponential which can be brought in canonical form by applying the Tschirnhaus transform [23]. The dynamics of the Sw (x, y, a 3 , z) beam can then be described in terms of again an optical swallowtail catastrophe Sw (A 1 , A 2 , A 3 ), whose control parameters are functions of x, y, a 3 , z in the form of: Sw(x, y, a3, z) = exp [iψ] · Sw (A1, A2, A3) , where(4) A1 = A1(x, y, a3, z), A2 = A2(y, a3, z), A3 = A3(a3, z), ψ = ψ(x, y, a3, z).
The complete analytical expression for the dynamics are given in the SM. Plotting (2) for different z-positions allows visualizing the dynamics of a Sw (x, y, a 3 ) beam. We set a 3 = 0. The analytically calculated propagation is depicted in Fig. 2, top. Furthermore, we obtain the volume intensity distribution of the beam experimentally for different z-positions and show them in Fig. 2, bottom. We chose transverse feature sizes of x 0 = y 0 = 8 µm and a propagation distance of 20 mm, which covers the most interesting effects in the sketched intensity volume. Analysis of (2) reveals that at z = 0 mm the Sw (x, y, a 3 ) beam shows mirror symmetry with respect to y: where the asterisk denotes the 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 Sw (x, y, a 3 ) beam becomes apparent: During propagation, the transverse field at the origin (z = 0 mm) turns out to consist of a fast diffracting contribution that quickly vanishes due to its transverse momentum (e.g. visible at z ≈ 10 mm), subsequently revealing a field contribution (remaining field at z ≈ 20 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 Sw (x, y, a 3 ) beam [13,18]: Due to this spectral similarity, we express the Sw (x, y, a 3 = 0) beam as a convolution of a Pearcey beam with a dimensionally reduced swallowtail beam: Here, δ(...) is the delta-function. Note that the initial beam profile of the artificially designed 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 5.

Dynamic transformation of swallowtail to butterfly caustics
The dynamics of the 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 oder n may show structural similarities with higher-order caustics during propagation, we exemplary investigate a Sw (x, a 2 , y) beam, thus set a 2 = const. and chose α = 1, β = 3. Applying (3) and subsequently performing a Tschirnhaus transform with appropriate substitutions, we express the propagation of the Sw (x, a 2 , y, z) beam in terms of the higher-order static butterfly beam Bu (B 1 , B 2 , B 3  The explicit functions for the new control parameters are given in the SM. The butterfly beam is defined as Bu (a 1 , a 2 , a 3 , a 4 ) = C 6 (a 1 , a 2 , a 3 , a 4 ). We have condensed γ = (−z/z ey ) 1/6 , and used characteristic lengths z e as defined previously. We note that the propagation of a Sw (a 1 , x, y) beam can be expressed in a similar way in terms of static butterfly beams by following (3).
(3) is plotted for different z-values and illustrated in Fig. 3, top. We obtained the dynamics as well in the experiment, which is shown below. The transverse dimensions are x 0 = y 0 = 8 µm and the propagation distance is 20 mm. Due to the point-symmetry of the Sw (x, a 2 , y) beam (cf. (3)) Sw (x, a 2 , y, z) = −Sw * (x, −a 2 , y, −z) , and for illustrative reasons, we restrict the evaluation to z ≤ 0. The asterisk denotes the complex conjugation.
The plotted analytical description of the propagation of a Sw (x, a 2 , y) beam is in high agreement with the corresponding experimental realization. The transverse plane at z = 0 mm shows no mirror nor point symmetry, which can be seen in Fig. 3 and be proved by considering (1). However, during propagation the structure separates into two individual structures, showing point symmetry. This suggests that we can assume the Sw (x, a 2 , y) beam at the origin to consist of two major contributions: One fast diffracting part that breaks the symmetry in the initial plane, and a second one that remains during propagation and exhibits point symmetry. After a certain propagation distance (z ≈ 10 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 Sw (x, a 2 , y) beam with a Bu (x, a 2 , y, a 4 ) beam, whose z = 0 mm real space appearance [18] resembles the field distribution of the Sw (x, a 2 , y) beam after a certain propagation distance. Their spectra both correspond to cubic functions, thus Sw (x, a 2 , y) =Bu (x, a 2 , y, a 4 )·e −i[(x0kx) 6 −(x0kx) 5 +(x0kx) 4 ] .
(10) This leads to the real space expression which indicates an inverse convolution. Properties of the static Sw (x, a 2 , y) beam at the origin are mainly determined by the characteristics of the Bu (x, a 2 , y, a 4 ) beam, therefore revealing similar field distributions during propagation as the Bu (x, a 2 , y, a 4 ) beam (cf. [18]).

Dynamic decay of the swallowtail caustic to a lowerorder cusp
The propagation of the Sw (x, y, a 3 ) beam shown in Fig. 2 gives rise to further study the dynamics of the corresponding caustic. By parametrizing (2) and analyzing the structure of the caustics according to (2), we find the z-dependent dynamics of the caustic in the transverse plane. Fig. 4 illustrates the parametrized 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 at Fig. 2.
The surface appears to consist of two intertwined slices, merging point symmetrically at the origin. Each of them has a cuspoid form that changes with the 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 respective other slice. According to (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 higherdimensional control parameter space, the expected decay of the swallowtail catastrophe to a lower-order cusp is apparent. Though 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 that exhibit novel and unique propagation properties like e.g. high-intensities on curved paths.

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,24]. These singularities of the gradient map of potential functions manifest as caustics in light, and where extensively studied as natural phenomena in the late seventies and eighties [3,10,11,25,26]. Mapping catastrophes to light fields via the paraxial catastrophe integral of (1) is a well-known approach. However, it took until 2007 that spatial light modulators were used to create Airy beams [12]. This allows for the first time the 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 caustic beams depending on a single state parameter s to map higher-order catastrophes to the lower-dimensional initial transverse plane of paraxial beams and derived analytically a general equation for their Fourier spectra (SM) and propagation. Our approach connects two of the control parameters a with the transverse spatial coordinates (x, y) T , thus the caustic 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 caustic beams of order n can be calculated in terms of higher-order static caustic beams of order 2(n − 2).
We demonstrated this by analytically calculating and experimentally obtaining the dynamics of the Sw (x, y, a 3 ) and 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 their high-intensity rims near the caustics capable to form tailored structures are unique for each order of catastrophe and parameter set and pave the way towards advanced micromachining on tailored curves [27] and the realization of waveguides with a rich diversity of light guiding paths [14].
We perform a 2D Fourier transform with respect to two control parameters a α , a β . The n − 2 dimensional control parameter space a is spanned by the coordinates a j , where j = 1, ..., n − 2. We introduce k α , k β as Fourier frequencies of a α , a β , where α, β ∈ [1, ..., n − 2]. With a we denote all remaining control parameters a except a α , a β . We neglected 2π scaling factors for reasons of clarity.
Thus, the Kronecker deltas and delta functions in (1) determine which quadrants of the Fourier space are occupied by Fourier components and which distribution they obey, may it be polynomial or an expression with rational exponent.

Analytical expressions for the dynamics of the optical swallowtail catastrophes
With Eq. (3) in the main document, it is a straight forward task to calculate analytically the propagation of any order of optical catastrophe that was mapped to the 2D initial transverse field via Eq. (1). Since we