Spiral Caustics of Vortex Beams

: We discuss the nonparaxial focusing of laser light into a three-dimensional (3D) spiral distribution. For calculating the tangential and normal components of the electromagnetic ﬁeld on a preset curved surface we propose an asymptotic method, using which we derive equations for calculating stationary points and asymptotic relations for the electromagnetic ﬁeld components in the form of one-dimensional (1D) integrals over a radial component. The results obtained through the asymptotic approach and the direct calculation of the Kirchhoff integral are identical. For a particular case of focusing into a ring, an analytical relation for stationary points is derived. Based on the electromagnetic theory, we design and numerically model the performance of diffractive optical elements (DOEs) to generate ﬁeld distributions shaped as two-dimensional (2D) and 3D light spirals with the variable angular momentum. We reveal that under certain conditions, there is an effect of splitting the longitudinal electromagnetic ﬁeld component. Experimental results obtained with the use of a spatial light modulator are in good agreement with the modeling results.


Introduction
In the classical geometric optics light is assumed to propagate along light rays, which are straight lines in a uniform medium. One of the amazing phenomena in the geometric optics is the formation of caustics, which are formed as an envelope of optical rays. The geometric optics [1,2] can be looked at as the short-wave limit of the classical optics, with the latter relying on the asymptotic approximation of Maxwell's equations [2]. In the basis of short-wave asymptotics is the fact that the electromagnetic field is described by integrals of rapidly oscillatory functions.
To determine a light ray intensity, a vector Kirchhoff integral needs to be taken. The common assumption of the geometric optics that the intensity at a given point is only determined by an insignificant fragment of the diffractive optical element (DOE) is its pivotal point. In many cases, an increase in the integration region does not lead to a considerable change in the intensity at the observation point. In terms of mathematics, this fact is described by a stationary phase method [2]. However, the diffraction integral not always can be calculated with the stationary phase method as in some cases, its use leads to the appearance of irremovable singularities, with the electromagnetic field intensity at the point tending to infinity. This is what happens in the neighborhood of geometric caustics. Caustics have been studied in numerous works. In the general case, classification of all possible caustics was conducted by V. I. Arnold [3]. The use of the catastrophe theory in the study of caustics is reported in Refs. [4,5]. An approach alternative to that based on the integral representations of rapidly oscillatory functions is discussed in Ref. [6], where the Helmholtz equation is reduced to the solution of a chain of differential equations. The method is analogous to reducing the Helmholtz equation to the solution of an eikonalequation and a transfer equation.Although it has some advantages, this approach only allows one to find the field in the neighborhood of non-singular caustics. Asymptotic methods for analyzing equations and the field near singular points are set forth in Refs. [7][8][9][10]. In this work, we study caustics of spiral vortex beams. The vortex beams are special in that the ray structure of their field depends on the wavelength [11]. Unlike Ref. [11], where classical vortex beams that form axisymmetric caustics were studying, in this work, we look into light fields with a spiral caustic.
Optical elements capable of generating spiral intensity patterns have been treated in a number of works [12][13][14][15][16][17][18][19][20][21][22][23][24]. In particular, E. Abramochkin and V. Volostnikov [12,13] proposed generating plane curves by means of astigmatic transforms. In Ref. [14], a carrier frequency method in combination with several parabolic wavefronts was applied to generating free-form 3D intensity distributions composed of light spots. A similar approach based on the superposition of diffraction-free beams, the spatial spectrum of which is determined on a narrow ring, was discussed in Refs. [15][16][17] for shaping different 3D light curves. A special combination of conical and vortex phases [18,19] makes it possible to form a light spiral in the focal plane of a lens. Note that conventional diffractive spiral axicons [20] and more complex refractive analogs [21] can be used to generate 3D spiral intensity in the near diffraction zone. The simplest optical element for shaping a spiral intensity distribution is a power-exponent phase plate [22,23] or a generalized spiral phase plate [24].
In this work, we propose DOEs with special angular dependence structure to generate 2D and 3D light spirals with the variable angular momentum. The considered approach expands the family of beams with optical vortices, which have found their applications in various fields including optical manipulation and laser structuring [25][26][27][28][29][30][31][32].

Electromagnetic Field on a Curved Surface
When tackling the problem of diffraction by a 3D solid (e.g., a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined/stitched [33][34][35][36][37]. Below, we discuss an asymptotic method for calculating the tangential and normal components of the electromagnetic field on a free-shape curved surface.

Coordinate Systems
Hereinafter, we utilize the following coordinate systems. The Cartesian coordinates (u, v, 0) on the DOE under study and the Cartesian coordinates (x, y, z) in the field computation domain.
The polar coordinates (ρ 0 , θ 0 , 0) on the DOE under study: The cylindrical coordinates (ρ, θ, z) in the field computation domain: Although the cylindrical coordinates are widely used, there are problems where the field needs to be analyzed in the local coordinates on a more complex-shaped surface. One such problem deals with calculating light forces exerted upon a free-shape solid.

2.
2. An Incident Wave in the Coordinates (x, y, z) The electric field of a plane wave may be given by where (α, β) are the spatial spectral coordinates linked with the Cartesian coordinates (x, y, z), E j (α, β) is the j-th component of the spatial spectrum of the electric field.
In these coordinates, the electric field of an incident plane wave on the surface of interest takes the form: where the matrix elements in (7) contain parametric equations that describe the surface in (3) and N j (ξ, η) are the components of the normal vector to the surface. The electric field of a plane wave may also be defined as where k 0 = 2π/λ 0 is the wavenumber for the wavelength of light λ 0 in vacuum and ε, µ are the dielectric permittivity and magnetic permeability of the medium. The functions A(α, β), B(α, β) define the contribution of the TM-and TE-components of the electromagnetic field. These may be replaced with functions E x (α, β), E y (α, β), which are given by In these coordinates, the electric field may also be represented as E(ξ, η) = P(ξ, η)E TME (α, β) A(α, β) B(α, β) × × exp{i[αX(ξ, η) + βY(ξ, η) + γ(α, β)Z(ξ, η)]}, Photonics 2021, 8,24 4 of 20 E TME (α, β) = E TMZ (α, β), E TEZ (α, β) (13) In the Cartesian coordinates, the magnetic field takes the form: In the coordinates linked with the surface of interest, the magnetic field is given by

Computation of the Field: General Case
In the general case of a field composed of superposition of plane waves, the field on the surface is described by the expression: where the functions A(α, β), B(α, β) are expressed using Equation (11) through the transverse components of the spatial spectrum, E x (α, β), E y (α, β). The spatial spectrum is derived from a Fourier transform of the transverse components of the electric field on the DOE: Then, the field on the surface takes the form: In the general case, the relations in (20) may be rewritten as: where K, L are the column-vectors. Equation (21) may be given in a different form: where the matrix G is defined as follows: Integral (23) can be calculated using a stationary phase method: where α s , β s are derived from the solution of a stationary point equation:

Computing the Field in the Neighborhood of a Spiral Caustic Surface
An interesting type of light caustics is a spiral. Light spirals can find uses not only for trapping, guiding, and micromanipulation of microparticles [17,29], but also for laser-aided surface structuring [19,38].
Let us analyze a wavefront that has 3D caustics, whose eikonal function may be generally written as: where the functions r 0 (θ 0 ) and f (θ 0 ) characterize the type of the 3D spiral. Additional vortex phase singularity exp(imθ 0 ) in (27) lets us obtain a vortex spiral wavefront. The vortex component guarantees the intensity null on the optical axis [11,39,40].
The complex amplitude of the electromagnetic field in the initial plane (z = 0) is given by Introducing the cylindrical coordinates in the observation plane: expressions for the Cartesian components of the field may be written in the form: where

Asymptotic Relations for the Diffraction Integral in the Neighborhood of a Spiral Caustic
In this section, we calculate the field in the neighborhood of a caustic surface, which is defined by parametric equations [11,41]: where and S is the solution of a quadratic equation: May the relations for r 0 (θ 0 ) and f (θ 0 ) be linear: A stationary point relative to the angle θ s will be derived from the equation: In this case, the field components in (30) and (31) take the form: where σ s = ±π/4 [2], The expressions in Equations (39) and (40) are 1D integrals. We note that without the use of an asymptotic approach, it is not possible to reduce the Kirchhoff integral to a 1D integral even if both the field and the DOE are radially symmetric. Based on the expressions derived, the number of computing operations can be reduced compared to the straightforward calculation of the Kirchhoff integral in (30), (31).

An Analytical Solution for an Annular Caustic
To verify the above-derived expressions, let us analyze a particular case for which analytical expressions can be deduced. If r 1 = r 0 , f 1 = f 0 , instead of a spiral, we get an annular caustic, for which the stationary point equation in (38) is essentially simplified: After rearrangements, Equation (44) takes the form: The solution to Equation (45) is: where It is worth noting that not all θ s from (46) satisfy the original equation because the solution contains additional roots. Therefore, prior to using a root, we need to check whether it satisfies the original equation in (44).
Next, making use of a stationary phase method, we find that E(ξ, η) contains only a single integral with respect to a variable ρ 0 : where R s and R s are defined in (41) and (42). The availability of analytical expressions for the stationary points in (46) makes the asymptotic analysis and design essentially simpler. The integrals in (48), (49) can be calculated analytically using a technique discussed in Ref. [42].
We note that Equations (48) and (49) may be utilized when calculating a field generated by an arbitrary axisymmetric wavefront at a distance from the optical axis.

Designing DOEs to Generate Spiral Caustics
As a rule, all these optical elements were designed and studied in a paraxial approximation. In the meantime, an optical element that implements the eikonal in (27) presupposes the use of a nonparaxial approach. Expression (27) implies shaping 3D light spirals. Let us analyze the eikonal in more detail, with the spiral longitudinally extended along the optical axis: where ∆ f (θ 0 ) = ( f 1 − f 0 ) 2π θ 0 = αθ 0 . Obviously, the longitudinal length of the formed spiral is determined by two terms in curly brackets in Equation (50), containing linear and nonlinear dependence on ∆ f (θ 0 ), and hence on the angle θ 0 .
Neglecting the constant last term, the DOE phase is approximately described as follows: Expression (51) contains terms nonlinearly dependent on the angle (in the braces) and the last term with the linear angular-dependence. The linear dependence on the angle corresponds to the classic optical vortex of the order p = −kα which in the general case is a fractional value [43]. The presence of this term is not necessary for the longitudinal distribution of the spiral, the nonlinear term is mainly responsible for this. An analogy can be drawn with 'perfect optical vortices' [44] which have a ring intensity distribution of the same radius regardless of the optical vortices present in the beam. However, the change in the value of the vortex singularity can be used to vary the speed of rotation over the ring of trapped particles [45].
The linear part of the angular dependence in an optical element described by Equation (50) corresponds to the classic optical vortex whose order is proportional to the ratio of the longitudinal spiral extent to the incident laser wavelength: Note, the magnitude of the difference ( f 1 − f 0 ) may be several thousand wavelengths (which is normally the case in the experiments), so the order of the vortex in (52) will be of the same magnitude. Too high an order of the extra optical vortex may lead to distortion of not only the on-axis but also the off-axis caustic [11]. We also note that if the order of the optical vortex is too high, even 'perfect optical vortices' that such elements form cease to be 'perfect', with the radius of the intensity ring becoming order-dependent [46]. Besides, in the experimental realization, when an order of the extra optical vortex is too high [47] the limited resolution of the manufacturing technologies or devices used for the generation of the desired vortex light fields (for example, spatial light modulators or digital micromirror devices) lead to distortion of the field. This challenge can be addressed via introducing a compensating optical vortex. The phase of a paraxial vortex lens, which focuses into different points of a light spiral in different planes, may be given by Equation (53) suggests that if m = ( f 1 − f 0 )/λ, the optical element is structurally simplified: Note, despite the compensation of the linear angular dependence, the nonlinear angular dependence, which is responsible for the longitudinal extension of the spiral, remains.
It stands to reason that the corresponding optical vortex may also be imbedded into the nonparaxial spiral DOE in Equation (50): Thus, the technique of shaping a 3D light spiral that uses an optical element implemented through the angular dependence of Equation (37) leads to the generation of both linear and non-linear optical vortices. Linear part of the angular dependence can be varied (both decrease or increase) by introducing an additional vortex phase into the optical element. This provides an extra degree of freedom in variations of the angular momentum of the field and in possible applications.

Results of the Numerical Simulation and the Experiment
Numerical modeling was conducted using a variety of propagation operators, including the asymptotic operators (39) and (40) derived herein and, by way of verification thereof, through the straightforward numerical integration of the vector Kirchhoff integrals (30), (31). Focusing into a plane spiral curve was numerically simulated using Equations (27), (36), (37) for the following parameters: f 1 = f 0 = 100 mm, m = 0, radiation wavelength 532 nm, while the rest of the parameters-including the DOE radius R d -were varied.
In the experiments, the proposed DOEs were optically implemented using a spatial light modulator (SLM) HOLOEYE PLUTO VIS. An experimental optical setup is shown in Figure 1. An incident linearly polarized Gaussian beam was expanded using a set of lenses L1 and L2 and then directed onto the SLM display. The reflected beam was phasemodulated and spatially filtered using a set of lenses L3, L4, and a pupil D. A video-camera mounted on an optical bench recorded intensity distributions at different distances from the DOE plane. To obtain circularly polarized laser beams, an additional quarter-wave plate put behind lens L4 was utilized in some experiments.
Note, despite the compensation of the linear angular dependence, the nonlinear angular dependence, which is responsible for the longitudinal extension of the spiral, remains.
It stands to reason that the corresponding optical vortex may also be imbedded into the nonparaxial spiral DOE in Equation (50): Thus, the technique of shaping a 3D light spiral that uses an optical element implemented through the angular dependence of Equation (37) leads to the generation of both linear and non-linear optical vortices. Linear part of the angular dependence can be varied (both decrease or increase) by introducing an additional vortex phase into the optical element. This provides an extra degree of freedom in variations of the angular momentum of the field and in possible applications.

Results of the Numerical Simulation and the Experiment
Numerical modeling was conducted using a variety of propagation operators, including the asymptotic operators (39) and (40)  , m = 0, radiation wavelength 532 nm, while the rest of the parameters-including the DOE radius d R -were varied.
In the experiments, the proposed DOEs were optically implemented using a spatial light modulator (SLM) HOLOEYE PLUTO VIS. An experimental optical setup is shown in Figure 1. An incident linearly polarized Gaussian beam was expanded using a set of lenses L1 and L2 and then directed onto the SLM display. The reflected beam was phase-modulated and spatially filtered using a set of lenses L3, L4, and a pupil D. A video-camera mounted on an optical bench recorded intensity distributions at different distances from the DOE plane. To obtain circularly polarized laser beams, an additional quarter-wave plate put behind lens L4 was utilized in some experiments. Results of the numerical simulation based on Equations (30) and (31) for the linear x-polarization (meaning that the y-component is zero) and the experiment are shown in Table 1. The inner and outer radii of the spiral were chosen to be r 0 = 1 mm, r 1 = 1.25 mm, and the DOE radius was varied: R d = 1.5 mm (upper row), R d = 0.7 mm (bottom row).
As can be seen (the first column in Table 1), the phase structure of the optical element has a spiral topology similar to the structures considered in [18,24]. It should be noted that the presence of a spiral structure in an optical element, which can be realized also by arrays of metallic nanoparticles [48,49], leads to the appearance of the angular momentum. In contrast to the classical vortex phase singularity, which ensures the formation of the orbital angular momentum (OAM) of the integer order (there is only one angular harmonic) [25,26], the spiral structure leads to the appearance of a set or a sequence of angular harmonics [24,48].
The angular harmonics spectrum can be used to estimate the OAM of a light field [50]: where c n are expansion coefficients of the analyzed field E(r, ϕ) by the angular harmonics exp(inϕ): In practice, the finite version of expression (56) is often used, as well as the OAM estimation by three arbitrary [51] or two maximum [52] coefficients (57). In addition, due to the conservation of the OAM value during propagation in free space [50], it is convenient to calculate it in the input plane, where the field is spatially bounded, and the integral in (57) has a finite radius limit.
Note that a total angular momentum [53], in addition to the OAM, can have a spin angular momentum (SAM) related with the polarization state of a field [54]. Moreover, in various situations, spin-to-orbit angular momentum conversion can occur [55][56][57][58]. However, here we consider only the OAM associated with the phase structure of the beams. Table 1. Numerical modeling and experiment on focusing into a single-turn spiral curve for the linear x-polarization.            Figure 2 shows in more detail the x-components of the fields presented in Table 1 Photonics 2021, 8, x FOR PEER REVIEW 12 of 25  Figure 2 shows in more detail the x-components of the fields presented in Table 1 and the normalized moduli of the expansion coefficients n c of the input fields. As can be  Figure 2 shows in more detail the x-components of the fields presented in Table 1 and the normalized moduli of the expansion coefficients n c of the input fields. As can be  Figure 2 shows in more detail the x-components of the fields presented in Table 1 and the normalized moduli of the expansion coefficients n c of the input fields. As can be  Figure 2 shows in more detail the x-components of the fields presented in Table 1 and the normalized moduli of the expansion coefficients n c of the input fields. As can be Figure 2 shows in more detail the x-components of the fields presented in Table 1 and the normalized moduli of the expansion coefficients |c n | of the input fields. As can be seen, in the first case (the first line in Figure 2) the phase practically does not change along the formed light spiral and the OAM is close to zero (µ = −0.03), and in the second case (the second line in Figure 2) there is a vortex phase singularity along the light spiral and the OAM is close to unity (µ = 0.81). seen, in the first case (the first line in Figure 2) the phase practically does not change along the formed light spiral and the OAM is close to zero (μ = −0.03), and in the second case (the second line in Figure 2) there is a vortex phase singularity along the light spiral and the OAM is close to unity (μ = 0.81). Thus, even for plane/flat curves, the proposed method of the light spiral formation makes it possible to vary the OAM of the light field just by changing the limiting aperture. From the results given in Table 1, the DOE parameters and the polarization state are seen to have an essential effect on the intensity distribution of the longitudinal E-field component (the second column in Table 1). The polarization state is linked with the emergence of a zero intensity along a line perpendicular to the polarization plane. In the case under study, the linear x-polarization leads to the zero-intensity line along the Y-axis.

Optical Element Phase
Another interesting effect observed in the intensity pattern of the longitudinal E-field component  Table 1).
For the optical element in row 1 of Table 1, one can clearly see a spiral-shaped zero-phase line. In this region, the rays pass through the optical element without refraction and in parallel with the optical axis. Thus, in this region, the transverse coordinates of output and input points of the rays are the same, meaning that the difference multipliers are zeroed. Rays from other parts are redirected onto the focal  Table 1: (a) phase of the input field, (b) amplitude and (c) phase of the field in the focal plane, (d) normalized moduli of the expansion coefficients |c n | for the input fields in the first line (red) and in the second line (blue).
Thus, even for plane/flat curves, the proposed method of the light spiral formation makes it possible to vary the OAM of the light field just by changing the limiting aperture.
From the results given in Table 1, the DOE parameters and the polarization state are seen to have an essential effect on the intensity distribution of the longitudinal E-field component (the second column in Table 1). The polarization state is linked with the emergence of a zero intensity along a line perpendicular to the polarization plane. In the case under study, the linear x-polarization leads to the zero-intensity line along the Y-axis.
Another interesting effect observed in the intensity pattern of the longitudinal E-field component |E z (u, v)| 2 is splitting of the spiral curve. This effect is caused by the presence of difference multipliers in the Kirchhoff integral: (ρ cos θ − ρ 0 cos θ s ) and (ρ sin θ − ρ 0 sin θ s ). These multipliers become zero if the dimensions of an optical element are close to those of the light spiral generated (cf. lines 1 and 2 of Table 1).
For the optical element in row 1 of Table 1, one can clearly see a spiral-shaped zerophase line. In this region, the rays pass through the optical element without refraction and in parallel with the optical axis. Thus, in this region, the transverse coordinates of output and input points of the rays are the same, meaning that the difference multipliers are zeroed. Rays from other parts are redirected onto the focal curve at some angles. For the optical element in row 2 of Table 1, all its regions direct radiation to a light spiral at some angles to the optical axis. A comparison of the results in Table 1 shows that an increase in the DOE radius leads to a widened curve thickness. Besides, if the DOE radius is smaller than that of the curve under shaping, the light spiral shaped is increased in scale (see the last column).
Computations based on the straightforward numerical integration of the vector Kirchhoff integrals (30), (31) are essentially time-consuming. This has prompted the use of diverse algorithms for accelerating the computing procedures [59,60]. The asymptotic relationships in Equations (39)- (40) we proposed herein enable the number of computing procedures to be decreased due to reducing the 2D integrals to the 1D ones. Importantly, the asymptotic approximation has little or no effect on the computation accuracy. A comparison of modeling results when focusing into a single-turn plane spiral curve obtained through the straightforward integration of Equations (30) and (31) and the asymptotic relationships (39)-(40) for the circularly polarized incident light are depicted in Table 2.   From (36), the spiral radius is seen to linearly increase from 0 r to 1 r with increasing angle 0 θ from 0 to 2π. Obviously, the variation can easily be made non-linear. Increasing the number of spiral turns is more challenging. Here, two approaches may be used. With    From (36), the spiral radius is seen to linearly increase from 0 r to 1 r with increasing angle 0 θ from 0 to 2π. Obviously, the variation can easily be made non-linear. Increasing the number of spiral turns is more challenging. Here, two approaches may be used. With    Although the asymptotic approach implies one-dimensional integration, additional computation of the stationary point is required (38). Nevertheless, the number of operations to obtain the final result is less than with double integration (30)- (31). This is because the solution of Equation (38)  approach was on average 1.5 times less than with double direct integration in the case of software implementation in Matlab (R2014a) on a personal computer using a processor Intel (R) Core (TM) i5-3570K CPU @ 3.4 GHz.
From (36), the spiral radius is seen to linearly increase from r 0 to r 1 with increasing angle θ 0 from 0 to 2π. Obviously, the variation can easily be made non-linear. Increasing the number of spiral turns is more challenging. Here, two approaches may be used. With the first, a composite DOE is synthesized [61], with its various parts contributing to different turns of the spiral. The second approach relies on superposition (summation) of complex transmission functions of DOEs that form different turns of the spiral [14,62]. Benefits and drawbacks of the both methods are clearly seen from Table 3, which presents the modeling and experimental results for focusing into a two-turn plane spiral curve for circular polarization (intensities of the x and y-components are identical).              The vortex order can be decreased by decreasing the difference 1 0 f f − to several wavelengths. In particular, to get p = 10, we need 1 0 0.005 mm f f − = , i.e. if 0 100 mm f = , then 1 100.005 mm f = (the spiral radius is much greater than its on-axis length). To get a scale-proportionate spiral, we reduce its radius until several microns, then at  The vortex order can be decreased by decreasing the difference 1 0 f f − to several wavelengths. In particular, to get p = 10, we need 1 0 0.005 mm f f − = , i.e. if 0 100 mm f = , then 1 100.005 mm f = (the spiral radius is much greater than its on-axis length). To get a scale-proportionate spiral, we reduce its radius until several microns, then at  The vortex order can be decreased by decreasing the difference 1 0 f f − to several wavelengths. In particular, to get p = 10, we need 1 0 0.005 mm f f − = , i.e. if 0 100 mm f = , then 1 100.005 mm f = (the spiral radius is much greater than its on-axis length). To get a scale-proportionate spiral, we reduce its radius until several microns, then at  The vortex order can be decreased by decreasing the difference 1 0 f f − to several wavelengths. In particular, to get p = 10, we need 1 0 0.005 mm f f − = , i.e. if 0 100 mm f = , then 1 100.005 mm f = (the spiral radius is much greater than its on-axis length). To get a scale-proportionate spiral, we reduce its radius until several microns, then at For sectored elements (rows 2 and 3 of Table 3), the ratio of sector areas needs to be fitted thoroughly to ensure that the thickness and intensity of different turns is the same. For composite elements (rows 2 and 4 in Table 3), there is no such requirement and the curve lines are thinner because radiation from the entire element contributes to each element of the curve. As a drawback, we should mention interference effects leading to the generation of extra segments in the spiral.
Given the above-described parameters and putting f 1 − f 0 = 1 mm, the extra vortex will have a very high order in (52): p = 1 mm/0.000532 mm ≈ 1879. An experimental implementation of vortices with high-order topological charge is feasible [40] but requires synthesizing a multi-level DOE with over 2000 × 2000 pixels.
The vortex order can be decreased by decreasing the difference f 1 − f 0 to several wavelengths. In particular, to get p = 10, we need f 1 − f 0 = 0.005 mm, i.e., if f 0 = 100 mm, then f 1 = 100.005 mm (the spiral radius is much greater than its on-axis length). To get a scale-proportionate spiral, we reduce its radius until several microns, then at f 1 − f 0 = 10 µm ÷ 50 µm, we find that p = 10 ÷ 50. Modeling and experimental results on focusing light into a 3D microspiral are presented in Table 4. Table 4. Focusing into a single-turn 3D micro-spiral curve: modeling vs. experiment.   Table 4. Focusing into a single-turn 3D micro-spiral curve: modeling vs. experiment.  Table 4. Focusing into a single-turn 3D micro-spiral curve: modeling vs. experiment.  Table 4. Focusing into a single-turn 3D micro-spiral curve: modeling vs. experiment.  Table 4. Focusing into a single-turn 3D micro-spiral curve: modeling vs. experiment.       Let us analyze in more detail the OAM of 3D microspirals. Figure 3 shows the results of varying the OAM for the 3D microspiral with p = −10. The initially formed spiral has a linear vortex component of the order p = −10, which is clearly seen from the phase structure in the focal plane (the first line in Figure 3). However, the OAM is somewhat different from this value (µ = −9.92) because it also includes a non-linear dependence on the angle, which ensures the longitudinal extent of the spiral. When using an additional vortex phase, the value of the OAM can be increased or decreased: at m = −5 (middle line in Figure 3) an increase occurs (µ = −14.92), and when m = 10 (lower line in Figure 3) we can see almost complete compensation (µ = 0.08). The residual value of the OAM corresponds to the non-linear part that provides the 3D character of the spiral. We suppose this value is related to the ratio of the spiral gap to the average spiral radius.

Element
The spiral length (a distance over which the spiral remains 'unblurred' and sharp) is clearly seen to increase with increasing f 1 − f 0 . In this case, the optical vortex order p, Equation (52), grows proportionally, as is evident from the structure of the optical elements (column 1 of Table 4).
If the vortex order p is too high, the resulting light caustic gets distorted. This effect can be observed in the bottom row of Table 4. In this case, becoming sufficiently high, the vortex component begins to affect the intensity distribution in the light spiral, 'dispersing' the light energy from the central part closer to the periphery. In particular, from the bottom row of Table 4, the spiral inner radius is seen to increase, getting close to its outer radius.
In addition, the experimental implementation of the phase with large orders of the vortex using SLM is problematic. The minimum zone size based on 3 pixels per zone for the used SLM is 3 × 8 µm = 24 µm.
To avoid these adverse effects, we use nonparaxial spiral DOEs of Equation (55) with additional compensating vortex phase. Table 5 shows modeling results for such DOEs.
somewhat different from this value (μ = −9.92) because it also includes a non-linear dependence on the angle, which ensures the longitudinal extent of the spiral. When using an additional vortex phase, the value of the OAM can be increased or decreased: at m = −5 (middle line in Figure 3) an increase occurs (μ = −14.92), and when m = 10 (lower line in Figure 3) we can see almost complete compensation (μ = 0.08). The residual value of the OAM corresponds to the non-linear part that provides the 3D character of the spiral. We suppose this value is related to the ratio of the spiral gap to the average spiral radius.  The results in Table 5 clearly show that in the initial plane ( f 0 = 100 mm), the initial spiral point located on the radius r 0 = 1 mm is brightest, whereas in the planes f 1 it is the final spiral point which is brightest, located at r 1 = 1.25 mm. Although the DOEs have hardly discernible difference between their phases, the extent of the 3D light spiral from the second DOE is twice as large. Table 5. Focusing into a single-turn 3D spiral curve: modeling vs. experiment. final spiral point which is brightest, located at 1 1.25 mm r = . Although the DOEs have hardly discernible difference between their phases, the extent of the 3D light spiral from the second DOE is twice as large.  Although the DOEs have hardly discernible difference between their phases, the extent of the 3D light spiral from the second DOE is twice as large. spiral point located on the radius 0 1 mm r = is brightest, whereas in the planes 1 f it is the final spiral point which is brightest, located at 1 1.25 mm r = . Although the DOEs have hardly discernible difference between their phases, the extent of the 3D light spiral from the second DOE is twice as large. 100 mm f = spiral point located on the radius 0 1 mm r = is brightest, whereas in the planes 1 f it is the final spiral point which is brightest, located at 1 1.25 mm r = . Although the DOEs have hardly discernible difference between their phases, the extent of the 3D light spiral from the second DOE is twice as large.

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.
The special DOE structure proposed herein because of angular dependence leads to the emergence of the angular momentum of the field. The possibility of the OAM

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.
The special DOE structure proposed herein because of angular dependence leads to the emergence of the angular momentum of the field. The possibility of the OAM

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.
The special DOE structure proposed herein because of angular dependence leads to the emergence of the angular momentum of the field. The possibility of the OAM

Conclusions
In this work, we have deduced asymptotic representations of diffraction integrals that describe behavior of vector electromagnetic fields in the neighborhood of caustics from spiral DOEs.
The distribution of the electric and magnetic fields on a curved surface needs to be known when calculating light forces exerted on a free-form convex conducting microparticle. Besides, when solving a problem of diffraction by a solid (a ball, a cylinder, etc.) bounded by a curved surface, the tangential components of the electromagnetic field need to be joined. The relationships deduced herein can be used for calculating the Poynting vector when evaluating forces exerted on a microparticle trapped in an optical field.
Equations for calculating stationary points and asymptotic relationships that describe components of the electromagnetic field as 1D integrals over a radial variable have been deduced. Without the use of an asymptotic approach, the Kirchhoff integral cannot be reduced to a 1D integral even if both the incident field and the DOE are radially symmetric. One-dimensional asymptotic integrals that describe a vortex DOE to focus into a 3D light spiral have been derived.
The special DOE structure proposed herein because of angular dependence leads to the emergence of the angular momentum of the field. The possibility of the OAM variations due to an additional vortex phase is shown. This provides an extra degree of freedom in controlling the motion of the trapped particle along a 3D spiral curve.