Composite Diffraction-Free Beam Formation Based on Iteratively Calculated Primitives

To form a diffraction-free beam with a complex structure, we propose to use a set of primitives calculated iteratively for the ring spatial spectrum. We also optimized the complex transmission function of the diffractive optical elements (DOEs), which form some primitive diffraction-free distributions (for example, a square or/and a triangle). The superposition of such DOEs supplemented with deflecting phases (a multi-order optical element) provides to generate a diffraction-free beam with a more complex transverse intensity distribution corresponding to the composition of these primitives. The proposed approach has two advantages. The first is the rapid (for the first few iterations) achievements of an acceptable error in the calculation of an optical element that forms a primitive distribution compared to a complex one. The second advantage is the convenience of reconfiguration. Since a complex distribution is assembled from primitive parts, it can be reconfigured quickly or dynamically by using a spatial light modulator (SLM) by moving and rotating these components. Numerical results were confirmed experimentally.


Introduction
The term "diffraction-free beam" was introduced to denote a laser beam propagating along the optical axis without changing the transverse distribution, i.e., without the influence of diffraction effects. The most famous among diffraction-free beams are the Bessel modes [1][2][3], which are the solution of the Helmholtz equation in cylindrical coordinates. In addition, Mathieu beams [4,5] are famous for the elliptic coordinate system and parabolic beams [6,7] for the parabolic coordinate system, as well as various generalized beams [8][9][10]. A general property of classical diffraction-free beams is the concentration of the spatial spectrum on a narrow ring. Such a property is often used to generate diffraction-free beams [2,11]. There are also other beams with diffraction-free properties whose spatial spectrum differs significantly from a narrow ring. They include Airy beams, Olver beams, and their modifications [12][13][14][15].
Among the effective methods of forming diffraction-free beams are the applications of axicons [16][17][18], diffraction optical elements (DOEs) [19][20][21], or spatial light modulators (SLMs) [22][23][24]. With these approaches, in contrast to the focusing of a narrow ring [2,11], a significant part of the energy of the incident beam goes to the formation of a diffraction-free beam. However, in some cases, additional coding of the calculated complex amplitudes into a phase-only mask may be required [25,26]. Another simple method of energetically efficient formation of various diffraction-free beams is based on a partial diaphragm of the annular light distribution [27,28], formed, for example, by a tandem of an axicon and a lens [29,30], or a toroidal lens [31,32].
Diffraction-free beams play a special role in holographic optical tweezers for controlled laser manipulation of an ensemble of microparticles. Structured laser beams with a complex intensity and phase distribution make it possible to hold and move nano-and micro-sized trapped objects in space. In this case, it is especially important to have the ability not only to form beams with a given transverse field distribution but also to quickly rearrange the structure of the beam. Usually, in this case, SLMs are used for the generation and dynamic movement of a set of focused light points [49][50][51]. Simple optical elements shaping various contour beams [52][53][54] and autofocusing beams [55][56][57][58] are also known. However, this does not ensure the diffraction-free nature of the beams.
The design of optical elements to generate complex diffraction-free patterns, as a rule, requires the use of iterative algorithms [59][60][61][62][63]. Although iterative methods seem to be a universal tool, they do not have global convergence, i.e., only a certain local minimum is guaranteed (one of the possible solutions with some accuracy). The disadvantage of an iterative approach is also the fact of repeated use of the direct and inverse operator at each iteration, which leads to significant time and computational costs.
In this paper, we propose a compromise approach, when a complex image is assembled from a set of diffraction-free primitive pictures (for example, a square or/and a triangle), for which the corresponding DOEs are pre-calculated and optimized using an iterative algorithm based on the narrow ring spatial spectrum. Based on the superposition of such DOEs, taking into account the rotation and displacement of primitives, it is possible to dynamically (for example, using SLM) generate various more complex diffraction-free patterns.

Theoretical Background
It is known that a light field bounded by a narrow annular diaphragm is transformed into a diffraction-free beam when focused by a lens. The length of the diffraction-free region in the longitudinal direction is determined by the following parameters [2]: where f is the focal length of the lens, d is the middle radius of the annular diaphragm, and the radius R 0 is the smaller value of the lens radius or the effective radius of the diffraction pattern [64]: where λ is the radiation wavelength, ∆d is the width of the ring slit. Thus, the length of the diffraction-free region defined by Equation (1) decreases both with increasing the radius and the width of the ring slit. Thereby, the diffraction-free beam is formed precisely by a very thin ring (∆d → 0). In an ideal case, when the spatial spectrum is the annular delta function, an arbitrary diffraction-free light field can be described by the Whittaker integral [65][66][67]: where α t < 1 is the constant corresponding to the radius of the annular spatial spectrum, and A(ϕ) is an arbitrary angular function on the annular spatial spectrum. Note, the pattern of the angular spectrum may be obtained as the result of an iterative solution to the inverse problem [68].
However, in the case of an infinitely thin spectral ring, energy losses are inevitable, and the convergence of the iterative process of the calculation deteriorates due to the small number of degrees of freedom [59,69,70]. That is why, in the given work, we consider in fact pseudo-diffraction-free beams [61,71], which have all the properties of diffraction-free beams at a certain propagation distance defined by Equation (1) related to the size of the bounding aperture.

Iterative Algorithm
We use an iterative algorithm that was previously successfully used to calculate a DOE shaping a diffraction-free beam of the "light sheet" type [72] suitable for fluorescence microscopy [73,74]. The iterative algorithm is based on the modeling of the action of the lens using the Fourier transform and superimposing the annular diaphragm in the spectrum plane ( Figure 1). The iterative algorithm presented in Figure 1 consists of the following main blocks: truncation of the complex field in the spectral plane by an annular diaphragm, modeling, action of the lens using the Fourier transform and replacement of the intensity by given one in the focal plane. The stopping criterion for the iterative algorithm is determined by the root mean square error (RMSE) of the generated intensity distribution from the given (primitive) picture. When the RMSE falls below the specified threshold MSE 0 , the calculation is finished . Figure 2 presents a block diagram of an iterative algorithm for the formation of a diffraction-free structured laser beam with a transverse intensity distribution in the form of a given primitive (in particular, "triangle"). The iterative process is described in more detail below.
At the zero iteration n = 0, the complex distribution is given in the form of a constant amplitude and a random phase with values in the range from 0 to 2π ("input field" in Figure 1, block "begin" in Figure 2): Then the input complex field is truncated by an annular diaphragm (ψ n (x, y) in Figure 1, block "1" in Figure 2):ψ n (x, y) = ψ n (x, y) · T(x, y), where T(x,y) is the ring-function: T(x, y) = 1, r s ≤ x 2 + y 2 ≤ r e , 0, else where r s = d − ∆d/2, r e = d + ∆d/2. The annular spectral distribution defined by Equation (5) is Fourier transformed simulating the lens action ( {ψ n (x, y)} in Figure 1, block "2" in Figure 2): The RMSE between the calculated intensity |Ψ n (u, v)| 2 and the given distribution (primitive) D(u, v) is estimated. If it does not exceed a given value of MSE 0 , then the exit criterion (block "3" in Figure 2) is successfully fulfilled, otherwise, we go to the next stage of the algorithm.
At the next stage, the complex distribution Ψ n (u, v) is replaced corresponding to the primitive distribution ("primitive" in Figure 1 and block "4" in Figure 2): Next, the Fourier transform is performed again ( {Ψ n (u, v)} in Figure 1, block "5" in Figure 2): The jump to the next iteration is carried out either at this step: or after reducing the input field to a pure phase distribution: In the second case, the convergence will be worse due to the loss of part of the degrees of freedom. Therefore, we used the first option in accordance with expression (10). However, in this case, at experimental implementation, additional encoding of amplitude information into phase information will be required [24][25][26]75].
Note that the proposed algorithm is very simple and will work effectively for fairly elementary distributions D(u, v), i.e., precisely for "primitives". We propose to form a more complex image from a set of various primitive pictures (taking into account rotation and displacement), for which the corresponding generating DOEs have already been calculated. This approach is convenient when using SLMs for the dynamic reconstruction of diffraction-free pictures, for example, in optical trapping [33][34][35][36] and laser material processing [37][38][39].

Simple Primitives
This section presents the results of calculating the complex transmission function ψ(x, y) using the iterative algorithm described in Section 2.1 with the aim of forming various primitive diffraction-free distributions D(u,v), for example, a triangle and a square ( Figure 3).  Table 1 shows the transverse distribution of the intensity |Ψ n (u, v)| 2 of the calculated field with a "triangle"-type primitive for different numbers of iterations n. It can be seen that the iterative algorithm after some number of iterations comes to a stagnation state, i.e., further iterations do not lead to significant changes. The complex transmission function ψ n (x, y) and the corresponding intensity distribution |Ψ n (u, v)| 2 at the n = 90 iteration are presented in Figure 4. To optimize the characteristics of the formed field, variations of the average radius of the circular spectrum d were considered (Table 2) at a fixed width ∆d. The calculation results for the "triangle" primitive ( Figure 5) show that increasing the radius of the ring d allows improving the convergence of the algorithm and minimizing the achieved RMSE for a given number of iterations. This is obviously due to an increase in the number of points inside the ring, which leads to an increase in the number of degrees of freedom. Note, the large radius of the spectral ring allows us to take into account high frequencies and clearly describe small details, for example, the borders of the primitive. Table 2. DOE phase ψ n=90 (x, y) with fixed maximal outer ring radius R 0 and variable width of the ring slit ∆d to form a diffraction-free "triangle" primitive.  Obviously, the larger the slit width (area), the higher the energy efficiency of the beam, and more pixels (degrees of freedom) can be used to improve the convergence of the algorithm. As can be seen from the numerical calculations (Table 2), an increase in the width of the ring slit ∆d leads to a decrease in RMSE. However, the diffraction-free properties of the beam will deteriorate. This clearly follows from Equations (1) and (2), since the length of the diffraction-free region z max decreases with the increase of the ring slit width ∆d. Recall that an ideal diffraction-free beam is formed at an infinitely narrow ring (∆d → 0) [65][66][67].
Therefore, in order to preserve the diffraction-free properties, we chose another way to increase the area of the ring, namely by increasing the middle ring radius d, keeping ∆d with a small fixed size. As the results of the calculations show (Table 3), increasing the middle ring radius d with a fixed ring slit width ∆d also leads to a decrease in RMSE, although not so significantly. This is an inevitable compromise between providing diffraction-free properties and forming a given image with acceptable accuracy. In a similar way, we calculated the complex transmission function ψ n=90 (x, y) for the "square"-type primitive D(u,v) ( Figure 6). From the comparison of the RMSE graph for a "square" (Figure 6d) and a "triangle" (Figure 5), it can be seen that error increases with the increase in the area of the primitive. Next, we target to form complex diffraction-free beams based on the superposition of already calculated primitives. To show the advantage of this approach, let us first consider the iterative calculation of a complex picture, for example, the "sign" picture ( Figure 7a), consisting of four triangles and one square. The calculation results are shown in Figure 7. Obviously, the RMSE for complex pictures is greater than the RMSE for individual primitives (Figures 5 and 6d). A more detailed comparison is shown in Figure 8.

Complex Patterns
For the formation of a diffraction-free beam with a complex predetermined transverse intensity distribution in the form of a grayscale image, we use the same iterative algorithm. Figure 9 presents a complex pattern D(u,v) with a grayscale "bird" image; amplitude and phase of the calculated DOE ψ n=90 (x, y), and the corresponding transverse intensity distribution in the focal plane |Ψ n=90 (u, v)| 2 . The results of the calculation of a DOE shaping the light field with the complex pattern "bird" (Figure 9a) with different numbers of iterations are shown in Figure 10. The resulting intensity distribution |Ψ n=10 (u, v)| 2 at n = 10 iteration has RMSE equal to 0.1147, which is insufficient for the exact reproducibility of the structured beam. The algorithm comes to stagnation approximately at the n = 50 iteration with RMSE = 0.0785, so at n = 90, the RMSE decreased slightly to 0.0723. Thus, for complex patterns, RMSE minimization requires significantly more iterations with RMSE values remaining significant (Figure 9) compared to RMSE for primitive patterns (Figure 8).

Superposition of Simple Primitives
In this section, we consider the construction of a diffraction-free beam with a complex predetermined transverse intensity distribution based on the superposition of previously calculated primitives. This method is effective for the optical trapping and manipulation of particles [33][34][35][36], as it provides dynamic restructuring of the configuration of optical traps [35,[49][50][51].
To design a multi-order DOE matched with several primitives (such as in Section 3.1), we use the method of spatial carrier frequencies [76]. It is known from the theory of diffraction and the properties of the Fourier transform that if the complex transmission function of DOE ψ(x, y) is multiplied by the carrier function exp[i(ax + by)], then in the focal plane of the element the intensity distribution: will be shifted by a distance proportional to a along the x-axis and by a distance proportional to b along the y-axis. The rotation of the complex function ψ(x, y) by ϕ 0 degrees corresponds to the rotation of the intensity distribution |Ψ(u, v)| 2 by ϕ 0 degrees also. We consider the multi-order DOE as a superposition of DOEs calculated for different primitives ψ (p) (x, y) taking into account both the rotation angles ϕ 0p and the displacements (a p , b p ): Figure 11 shows the results of the action of the multi-channel DOEs τ(x, y) calculated by Equation (12) and matched with the "triangle" primitives. Obviously, we can compose the multi-channel DOE τ(x, y) calculated by Equation (12) from different types of primitives. Figure 12 shows the multi-channel DOE matched with one "square"-type primitive at the central part and four a "triangle"-type primitives rotated by angles ϕ 0p = πp/2 and shifted from the central part. Figure 12d presents the synthesis scheme for such multi-channel DOE in detail.

Diffraction-Free Properties
To visualize the diffraction-free properties of a beam generated by the multi-channel DOE composed of the superposition of DOEs calculated for different primitives (Equation (13)), we will calculate the transverse intensity distribution at different distances from the focal plane of the lens.
Such intensity distribution can be calculated using the Fourier transform with the addition of the defocusing function as follows: where f is the focal length and z is the defocusing distance. Table 4 presents a comparison of numerical modeling results for a primitive "triangle" and for the composed diffraction-free beam (superposition of simple primitives) defocused from the focal plane (f = 100 mm) on different distances.
It can be seen from Table 4 that a diffraction-free beam based on the "triangle" primitive with a slit width ∆d = 0.2 R 0 (1st line) retains its structure better than for ∆d = 0.4 R 0 (2nd line): the deviation between the intensity distributions in the focal plane (∆z = 0) and at a distance ∆z = −75 mm to the focus is RMSEz = 0.0072 (for ∆d = 0.2 R 0 ) and RMSEz = 0.009 (for ∆d = 0.40 R 0 ). Table 4 shows that the composed diffraction-free beam (3 rd line) keeps well enough its structure up to 75 mm from the focal plane: the deviation between the intensity distributions in the focal plane (∆z = 0) and at a distance ∆z = −75 mm before the focus does not exceed 0.07. We expect that the generated diffraction-free beam retains its structure also after the focal plane. The RMSE between the intensity distribution at the focal plane and the intensity distribution defocused at 75 mm distance (after focus) does not exceed 0.06. The calculation results for the generated diffraction-free beam at defocus distance z in the range from −75 mm to 75 mm are presented in Figure 13. We expect that the generated diffraction-free beam retains its structure also after the focal plane. The RMSE between the intensity distribution at the focal plane and the intensity distribution defocused at 75 mm distance (after focus) does not exceed 0.06. The calculation results for the generated diffraction-free beam at defocus distance z in the range from −75 mm to 75 mm are presented in Figure 13.  Table 5 presents a comparison of numerical modeling results for a complex pattern of the "bird" image defocused from the focal plane (f = 100 mm) on different distances. It can be seen that a diffraction-free beam based on the complex pattern of the "bird" with a slit width Δd = 0.2 R0 (1st line) in the focal plane has a larger deviation from the given pattern D(u,v), than for Δd = 0.4 R0 (2nd line): RMSE = 0.07 and RMSE = 0.05, respectively. However, in the first case, when defocusing, the structure is better preserved: the deviation between the intensity distributions in the focal plane and at a distance Δz =−50 mm to the focus is RMSEz = 0.16 and RMSEz = 0.19, respectively.   Table 5 presents a comparison of numerical modeling results for a complex pattern of the "bird" image defocused from the focal plane (f = 100 mm) on different distances. It can be seen that a diffraction-free beam based on the complex pattern of the "bird" with a slit width ∆d = 0.2 R 0 (1st line) in the focal plane has a larger deviation from the given pattern D(u,v), than for ∆d = 0.4 R 0 (2nd line): RMSE = 0.07 and RMSE = 0.05, respectively. However, in the first case, when defocusing, the structure is better preserved: the deviation between the intensity distributions in the focal plane and at a distance ∆z = −50 mm to the focus is RMSEz = 0.16 and RMSEz = 0.19, respectively.

Experimental Results
For the experimental investigation of the designed DOEs, the optical scheme based on a reflective spatial light modulator (SLM) HOLOEYE PLUTO VIS (1920×1080 pixels, pixel size of 8 µm) was used (see Figure 14). Figure 14. The optical scheme for investigation of the designed DOEs forming diffraction-free laser beams: Laser is a solid-state laser (λ = 532 nm); L1, L2, L3, and L4 are spherical lenses (f 1 =25 mm, f 2 =150 mm, f 3 =500 mm, f 4 =125 mm); M1 and M2 are mirrors, SLM is a reflective spatial light modulator (HOLOEYE PLUTO VIS); D is a circular diaphragm, CAM is a video camera (TOUPCAM UHCCD00800KPA, 3264 × 2448 pixels, 1.67 µm pixel size) mounted on an optical rail.
A Gaussian beam from a solid-state laser (λ = 532 nm, P out = 20 mW) was extended and collimated with a combination of two lenses L1 and L2 with focal lengths of 25 and 150 mm. The collimated laser beam was directed onto the SLM with the help of mirrors M1 and M2. The optical system consisting of two lenses L3 and L4 with focal lengths of 500 and 125 mm as well as a circular diaphragm D was used for spatial filtering of the modulated laser beam reflected from the SLM. Then a video camera CAM (TOUPCAM UHCCD00800KPA; 1600 × 1200 pixels, with a pixel size of 3.34 µm) mounted on an optical rail was used to record the intensity distributions of the formed laser beam at different distances. In the experiments, we used a pure-phase DOE resulting from the encoding of the amplitude-phase transmission function presented in Figure 12.
For the encoding, we used the partial encoding method [75]. Similar to other techniques of amplitude encoding, the used method is sensitive to the inhomogeneity and irregularity of an encoded amplitude distribution. Of course, this introduces the encoding error. In fact, in our work, we used the amplitude encoding only to encode zero amplitude values in the central region of the elements (which have a localized ring structure) to direct light from this region to higher diffraction orders. A limitation of this encoding method is the pixel resolution of the SLM or the manufacturing resolution of the fabricated DOE used to implement the encoded phase mask because the number of pixels increases during coding.
The intensity distributions formed by the investigated pure-phase DOE are shown in Figure 15. The experimentally obtained results are in good agreement with numerically obtained ones. The experimental results qualitatively confirmed numerical calculations and demonstrated the possibility of implementing the calculated ring distribution (when the informative part of the field is concentrated in a narrow region) using SLM, a device with a fixed resolution and number of pixels.

Discussion
Let us compare the RMSE for the "sign"-pattern (referred to as "pattern") and the superposition of primitives composed as "sign"-picture (referred to as "composition") when we used the iterative algorithm with N = 90 iterations ( Figure 16). Figure 16. Comparison of the RMSE for the "sign"-pattern (referred to as "pattern") and the superposition of primitives composed as "sign"-picture (referred to as "composition") (a) and difference in percentage (b) for N = 90 iterations with a fixed middle ring radius d = 0.9R 0 .
As seen from Figure 16, the use of the multi-order DOE of Equation (13) composed of primitives in comparison with a DOE calculated for the complex pattern makes it possible to reduce the RMSE already at the first iteration by almost 30% (from 0.044 to 0.034). Moreover, the advantage in RMSE remains during 50 iterations. After the 50 th iteration, the situation changes since the iterative algorithm for primitives enters stagnation (after about the 30 th iteration, Figure 16a). Thus, the use of a composition of primitives is convenient due to the rapid (for the first few iterations) achievements of an acceptable error (less than 0.02). Changing the structure of the composition does not require additional iterative calculation, as it is achieved by simple operations of rotation and displacement of individual parts of the composition.
From a practical point of view, the results obtained in the article will make it possible to carry out optimization structured light traps, introducing dynamic corrections to the formed intensity distribution with diffraction-free properties.

Conclusions
In this work, an iterative algorithm based on a ring spectrum and a lens was developed for the formation of a complex structure diffraction-free beam. DOEs were calculated and studied for various values of the parameters responsible for the number of iterations and the radius of the circular spectrum. It was shown that, for primitive patterns, the iterative algorithm provides the rapid (for the first few iterations) achievement of an acceptable RMSE (<0.02) and then comes to a stagnation state when further iterations do not lead to significant changes in the error. The calculated RMSE does not exceed 0.006 for a triangletype primitive, 0.015 for a square-type primitive, 0.06 for a superposition of primitives, and 0.07 for a complex grayscale image.
To calculate the structured laser beam given from the superposition of primitives, a multi-order DOE based on the method of spatial carrier frequencies was developed. A diffraction-free beam with a more complex transverse intensity distribution was formed based on the superposition of the DOEs pre-calculated for different primitives (taking into account rotation and shifting). The proposed method of synthesis of multi-order DOEs is convenient due to the rapid or dynamic (when using SLM) reconfiguration of the composition structure without additional iterative calculation, as it is achieved by simple operations of rotation and displacement of individual parts of the composition.
To visualize the diffraction-free properties of the calculated beam, the transverse distribution of the intensity at different distances from the focal plane of the lens was calculated. We estimated that the RMSE between the focal intensity distribution (for a lens with a focal length of 100 mm) and the intensity at a distance of 75 mm from the focal plane does not exceed 0.07. To test the diffraction-free properties of the calculated beam, an optical experiment was carried out and the intensity distributions were recorded before and after passing through the lens. The experimentally obtained results are in good agreement with the numerical ones.