Optical tweezers have played an important role in life sciences and medical applications^[1]. However, conventional tweezers based on tightly focused Gaussian beams rely on high-numerical-aperture optics, which can limit trapping flexibility. Continuous high-intensity exposure may also cause localized heating or photodamage in sensitive biological environments. Optical tweezers based on structured light address these challenges by modulating the amplitude, phase, and polarization of light to engineer optical force fields. Such beams enable multi-particle manipulation, dynamic control, and tailored force profiles^[1–4]. Among these, abruptly autofocusing beams^[5–8], derived from self-accelerating catastrophe beams, exhibit unique propagation characteristics. They maintain low intensity along most of the optical path and sharply focus energy at a focal plane even without a lens or nonlinearity. This sudden concentration of light generates strong intensity gradients for rapid, localized trapping^[9,10]. Notably, previous work^[11] demonstrated that autofocusing beams can achieve higher trapping stiffness than Gaussian beams at the same power, allowing more efficient trapping with less energy and potentially reducing photodamage.

To enhance the capabilities of autofocusing beams, researchers have developed many kinds of autofocusing beams carrying vortex phases, such as circular Airy vortex beams^[12,13], circular Pearcey vortex beams^[14,15], and multi-Airy vortex beams^[16]. A vortex phase makes beams carry orbital angular momentum (OAM), enabling rotational manipulation and potentially lowering photodamage by reducing the illuminated area through a hollow beam profile^[1,12–17]. However, most of these works focus on autofocusing beams generated from low-order catastrophes, leaving the potential of higher-order catastrophe-based vortex beams largely unexplored.

Therefore, based on our previous work^[4], we propose a type of abruptly autofocusing beam: dislocated-superimposed swallowtail vortex beams (DSVBs) via superposing multiple swallowtail catastrophes and vortex phases along dislocated directions and asymmetric azimuths. We investigated the propagation behaviors of DSVBs both experimentally and theoretically, analyzing their propagation dynamics through calculations of power flows, angular momentum densities (AMDs), and optical trapping forces. Our results demonstrate that DSVBs with opposite topological charges (TCs, ±l) exhibit asymmetric rotating and autofocusing behaviors during propagation, resulting from an asymmetric distribution of OAM. Notably, when the superposition numbers and TCs take the same value, multiple solid focuses form during propagation, differing from typical vortex beams, which usually focus to a single hollow spot. These findings advance our understanding of the focusing properties and OAM of autofocusing beams, offering more possibilities for applications in optical trapping and manipulation.

According to catastrophe theory, a caustic field Cn(a→) can be represented by the standard diffraction catastrophe integral^[19], expressed as Cn(a→)=∫−∞+∞exp[iPn(a→,s)]ds. The canonical potential function, denoted as Pn(a→,s)=sn+∑j=1j=n−2ajsj, is governed by dimensionless control parameter a→=(a1,a2,aj) and state parameter s (n represents the degree of s). Here, a→ represents a vector comprising all dimensionless control parameters aj, where j=1,2,…,n−2(n≥3). When the number of aj is less than four, the structurally stable singularities, known as elementary optical catastrophes, can be classified into seven types, ordered by the complexity of their potential functions: fold, cusp, swallowtail, butterfly, elliptic umbilic, hyperbolic umbilic, and parabolic umbilic^[19]. Among these, higher-order swallowtail catastrophe (n=5), denoted as C5(X,Y,Z), can be formally defined through a catastrophe integral^[7,16], C5(X,Y,Z)=∫−∞+∞exp[i(s5+Zs3+Ys2+Xs)]ds. Here, X, Y, and Z represent dimensionless coordinates related to Cartesian coordinates, corresponding to aj.

As illustrated in Fig. 1, we superimpose multiple 2D swallowtail catastrophes along dislocated directions and asymmetric azimuths to generate the DSVBs. After superimposing a vortex phase, the field of DSVBs at the initial plane can be expressed as ψ(x,y,0)=∑m=1NC5(Xm,0,0)×C5(0,Ym,0)[x+iyx2+y2]l,(1)[XmYm]=[cosθmsinθm−sinθmcosθm][x/x0y/y0]+[c0c0].(2)Here, m=1,2,3,…,N denotes the order of the superimposed swallowtail beams, x0 and y0 represent the transverse scale factors, and c0 is a dimensionless constant corresponding to x/x0 and y/y0, which determines the distance of each superposed catastrophe beam from the center [Fig. 1(a)] and affects the dislocation between the main lobes of DVSBs. This spatial dislocation, controlled by c0, modulates the transverse beam structure and tunes the axial positions of the focuses. l is the topological charge for introducing a central vortex phase into DSVBs via [x+iyx2+y2]l, while the angular separation between the superposed components is defined by θm=2π(m−1)/N [Fig. 1(b)]. Notably, the main lobes of different swallowtail components do not converge during propagation, as their transverse acceleration directions are non-intersecting. This behavior stems from the asymmetric intensity distribution of each superposed catastrophe beam governed by control parameters aj, which are set to be Xm or Ym in Eq. (1). Consequently, DSVBs experience a dislocated autofocusing propagation [Fig. 1(c)].

To further investigate beam propagation dynamics, we calculated the Fresnel diffraction integral of DSVBs by $$\begin{gathered} ψDSVB(x,y,z)=exp(ikz)iλz∬Rψ(x′,y′,z=0) \\ ⁢exp[i(x−x′)2+(y−y′)22z]dx′dy′, \end{gathered}$$(3)where λ is the wavelength and k=2π/λ denotes the wavenumber.

Following the method in prior works^[18–21], numerical simulations were performed on the propagation of DSVBs. In the simulation, the wavelength was set to λ=632.8nm, the input power P=1W, x0=y0=46.9μm, and c0=3. Besides, the propagation characteristics could be experimentally investigated using a setup analogous to the off-axis hologram method, similar to Refs. [18–20] [Fig. 2(a)]. The experimental configuration involved transmitting an expanded quasi-plane Gaussian beam. The spatial light modulator (SLM) was programmed with a precomputed complex hologram of DSVBs, generated by simulating the interference pattern between the DSVBs’ wavefront and a reference plane wave. A 4f system (comprising lenses L3 and L4) equipped with a spatial filter was then used to split the first-order diffracted light, modulated by the SLM, transforming the beam into DSVBs. The propagation of the generated DSVBs was captured and analyzed by a charge-coupled device (CCD).

Initially, to observe obvious rotational propagation of DSVBs from the intensity distribution, we began with small superposition numbers, for example, set N=4 to show structure features of individual superimposed catastrophe beams, while propagation dynamics of DSVB with different TCs were examined. Figures 2(b)–2(e) present the corresponding results for l=2 [Figs. 2(b1)–2(b3)], l=0 [Figs. 2(c1)–2(c3)], and l=−2 [Figs. 2(d1)–2(d3)] (insets are the simulation results). At the initial plane, DSVBs exhibit similar initial intensity distributions [Figs. 2(b1)–2(d1)]. However, the focal planes are located at z=133, 183, and 91 mm. Furthermore, as the absolute value of l increases, the central hollow region of the beams expands.

In particular, for l=0, the focal length of the beam is approximately about eight Rayleigh lengths [Fig. 2(c3)], which highlights its robust autofocusing behavior. Additionally, when the beams possess the opposite TCs (l=2 and l=−2), the intensity profiles at the focal plane and their sideview propagation exhibit distinct characteristics [Figs. 2(b2), 2(d2) and 2(b3), 2(d3)]. Their autofocusing and rotating propagation displays asymmetry (Video 1, Video 2, and Video 3), especially the rotational speed. As analyzed in Ref. [4], DSVB with l=0 still exhibits spontaneous autofocusing and forward rotating behavior [Fig. 2(c)], which results in such asymmetry propagation behavior of DSVBs. Moreover, increasing superposition number N will cause the intensity of each superposed beam to decrease, leading to a smaller intrinsic OAM of DVSBs (see Fig. 4) and further weakening of the rotational speed of DSVBs during propagation (see Video 4 and Video 5, which present the propagation of DSVBs with l=0,N=6,8).

Next, propagation of DSVBs was investigated under a special condition where the superposition number N equals the value of TC (l), i.e., N=|l|. Figure 3, Video 4, Video 5, and Video 6 present these results where DSVBs take N=−l=4 [Figs. 3(a1)–3(a3)], N=−l=6 [Figs. 3(b1)–3(b3)], and N=−l=8 [Figs. 3(c1)–3(c3)]. In all cases, DSVBs keep autofocusing and rotating propagation behaviors. As the beams propagate, multiple solid focuses with the largest intensity are generated and located at z=121, 119, and 129 mm. In addition, multiple solid focuses have also been detected in the cases for N=l=4, N=l=6, and N=l=8 (not shown here), although their peak intensity values are notably lower compared to the previous cases.

To explain the propagation dynamics of DSVBs via power flow, we calculated their transverse Poynting vectors S→⊥^[10,16]: S→⊥=i4η0k(ψ∇→⊥ψ*−ψ*∇→⊥ψ),(4)where η0=μ0/ε0 is the impedance of free space, μ0andε0 are the dielectric constant and magnetic permeability in vacuum, and ∇→⊥=∂∂xe→x+∂∂ye→y is the gradient operator.

Figure 4 shows how the power flows of DSVBs work. With TCs, power flows around the center, like a beam rotating forward. At the initial plane, power flows inward [Figs. 4(a1)–4(a3)]. At the focal plane, power flows balance outwards and inwards, with arrows along the hollow middle. Beyond the focal plane, power flows outwards and gets weaker (Video 9, Video 10, and Video 11). In addition, power flows of DSVBs with positive TCs (l>0) rotate clockwise, while those with negative TCs (l<0) rotate counterclockwise. This counterclockwise behavior also appears at l=0, suggesting transverse acceleration directions of lobes are different and N also affects this phenomenon because it tunes DSVBs’ OAM [Fig. 5]. For cases of N=|l|, power flows in the tangential direction become stronger with larger N or l [Figs. 4(b1)–4(b3)]. During propagation, power concentrates inwards and disperses outwards many times (see Video 12, Video 13, and Video 14). Thus, after such many dislocated superpositions, each beam’s intensity will be plus or minus at different specific locations on the propagation axis, also manifesting as multiple solid focuses.

OAM describes the angular momentum generated by the object as it rotates around an axis. Therefore, we can quantitatively calculate the OAM of DSVBs at different TCs and superimposed number N, and find more understanding for the asymmetric propagation and power flow distribution of DSVBs. Generally, we can calculate angular momentum density (AMD) j→ via linear momentum density p→^[4]: j→=r→×p→=r→×(E→×H→)/c2=jxe→x+jye→y+jze→z. r→ is the position vector. c is light speed. E→ is the electric field. H→ is the magnetic field. jx and jy are the transverse AMDs, and jz is the longitudinal AMD. Since OAM is the integration of longitudinal AMD jz over the entire space, i.e., Jz=∬jzdxdy, the average OAM per photon can be expressed as^[26]⟨L⟩=ckJzP=2ck∬jzdxdy∬|ψ|2dxdy.(5)

Figure 5 illustrates the longitudinal AMD (backgrounds) and OAM spectrum of DSVBs at the initial plane. AMDs with a blue background are usually negative; otherwise, the opposite [Figs. 5(a1)–5(f1)]. Generally, OAM follows the change of the modulating TCs. However, the averaged OAM ⟨L⟩ is 1.99ℏ at the initial plane for l=0, while 3.91ℏ and 0.02ℏ for l=2 and l=−2, respectively. Thus, DSVBs with l=0,2 present an enhanced rotation speed (Video 1, Video 2, and Video 3). As reported in Ref. [4], when N>1, the interaction of the transverse linear momentum across the branches of the catastrophe beam generates phase vortices, which manifest as a non-zero OAM, enabling self-rotation even when l=0. Moreover, the inherent OAM creates an OAM discrepancy between catastrophe superposition and vortex phase modulation, leading to asymmetric power flows and propagation characteristics for l=2 and l=−2.

To explain the above asymmetric phenomena, we also analyze the OAM spectrum of DSVBs via the methods introduced in Refs. [23–25], i.e., for any complex amplitude optical field, its OAM spectrum can be obtained by Fourier angular decomposition. Its OAM mode coefficients cl can be calculated as cl=12π∫02πψ(x,y,0)e−ilϕdϕ.(6)

From cl, we can calculate the power (or probability) distribution representing each OAM mode in the light field.

Figures 5(a2)–5(f2) present the OAM spectrum for six cases of DVSBs corresponding to Figs. 2, 3, and 4. In addition to the most weighted OAM mode, a small proportion of other OAM modes coexist [Fig. 5(e2)]. These additional components exhibit an asymmetric distribution around the maximum weight, which explains the asymmetric propagation and power flow distributions. Figure 5(g) shows the effect of TCs, where overall OAM shifts upwards (inset line graph). With both N and l modulation, an increase in N results in a reduction of OAM. In contrast, increasing l leads to an enhancement of OAM. Furthermore, the asymmetry in OAM is commonly occurring under each superposition number N.

DSVBs can be used for optical trapping and manipulation. Their transverse gradient force and scattering force acting on Rayleigh particles can be calculated by^[22]F→g=14ε0εmRe(α)∇|ψ2|e→r,F→s=16πcεm3k4∇|α2|S→⊥,(7)where α=4πR3(εp−εm)(εp+2εm) corresponds to the polarizability, R signifies the radius of the particle, and εp and εm represent the dielectric permittivity and the surrounding medium of the particle, respectively. Thus, we theoretically calculate the trapping force for a polystyrene bead (R=50nm) in the water (εp=2.5, εm=1.7).

Figures 6(a)–6(b) depict the gradient and scattering forces at the focal point. The DSVB with l=0 exhibits the strongest gradient force, while it has the smallest scattering force. In addition, the DSVB with l=−2 has a stronger scattering force, even though its OAM approaches zero. Furthermore, we calculated the trapping force distribution (F→trap=F→g+F→s). At the initial plane, all three cases show similar force distributions, with the force vectors pointing towards the inner parts of the main lobes [Fig. 6(c1)]. For l=0, trapping force at the focal point is nearly twice as large as that at the initial plane, and all the forces point inward [Fig. 6(c3)]. In contrast, in the other two cases, the trapping force at the focal plane is approximately equivalent in magnitude to that at the initial plane [Figs. 2(a3), 2(d3) and 6(c2)–6(c4)]. The force vectors point towards the inner sides of the lobes, where particles will be wrapped and spun forward, possibly reducing photodamage during optical manipulation^[10].

In this study, we propose a type of special autofocusing beam named dislocated-superimposed swallowtail vortex beams by modulating high-order catastrophe beams with vortex phases. Modulation of topological charge (l) and superposition number (N) significantly influences beam propagation and trapping. This is due to two key aspects: non-zero total OAM even when l=0 and asymmetric OAM distribution under modulation of opposite TCs (±l). Moreover, when N=|l|, DSVBs not only preserve non-zero OAM and rational propagation, but also exhibit multiple solid focuses at the center of symmetry during propagation. In this case, compared to conventional OAM modes that arise from vortex phases that generate hollow points, the OAM of DSVBs may offer multiple axial particle trapping at the center and further rotate additional particles around each central spot at the same time. These findings enhance our understanding of higher-order catastrophe beams and offer new possibilities for developing optical manipulation tools in biomedical applications^[27].