Generating Bessel-Gaussian Beams with Controlled Axial Intensity Distribution

: This paper investigated the di ﬀ raction of a Gaussian laser beam on a binary mask and a refractive axicon. The principles of the formation of a zero-order Bessel beam with sharp drops of the axial ﬁeld intensity edges were discussed. A laser optical system based on an axicon for the formation of a Bessel beam with quasi-uniform distribution of axial ﬁeld intensity was proposed. In the laser optical system, the inﬂuence of the axicon apex did not a ﬀ ect the output beam. The results of theoretical and experimental studies are presented. It is expected that the research results will have practical application in optical tweezers, imaging systems, as well as laser technologies using high-power radiation.


Introduction
Research on the methods of focusing light beams associated with the names of such scientists as Ernst Abbe, George Airy, and John Strutt have a long history and remain relevant to this day. Tight focusing of laser radiation, which allows one to obtain high power densities and directivity, is in demand in many problems of photonics and laser physics [1][2][3][4][5][6].
In practice, Gaussian laser beams are widely used. When they are focused using an optical system with spherical (or parabolic) surfaces, it is possible to concentrate laser radiation in the region corresponding to the beam waist diameter. However, in Gaussian beams, the wavefront is flat only in the waist region, the length of which is equal to the confocal parameter. With a decrease in the waist diameter of a Gaussian beam, its confocal parameter decreases [7][8][9].
By choosing the type and shape of the surfaces of the laser optical system elements, it is possible to correct the geometry of the wavefront. Correction methods can also be applied to increase the focusing length of the formed beam. In the end, when the lens takes the form of a cone, the phase front of the transmitted wave becomes conical. Such a conical lens is called an axicon [10]. The axicon forms a beam with a radial distribution of field intensity, which is described by the Bessel function of the first kind of the zero-order. Such beams are often used for scientific and applied problems [11,12], in which they have advantages over Gaussian beams [13,14]. The main ones are: (i) the transverse distribution of field intensity remains practically unchanged over a significant portion along the optical axis, (ii) the distribution of the field is restored behind obstacles, and (iii) excellent depth of field.
Axicon is a reliable optical element for the formation of a zero-order Bessel beam, including schemes with high-power laser radiation. However, it has disadvantages such as the dependence of the parameters of the generated beam on a specific element, the dependence of the intensity distribution of the output beam on the quality of element manufacturing (in particular, on the quality of the apex), and uneven distribution of axial intensity [15].
There is a huge number of different tools used for the formation of Bessel beams [16,17] which allow controlling axial intensity distribution. Diffraction [18,19] and holographic methods [20] are outside of the scope of this study due to the features of their implementation. Methods capable of realizing the required Bessel beam intensity profiles are widely used in microscopy [21], including biological imaging [22][23][24], as well as in applications of ultrashort laser pulses [25,26]. In almost all of the abovementioned applications, a Bessel beam with a controlled quasi-uniform axial intensity distribution is formed.
The only way to change the axial intensity is to apodize the beam incident on the element (axicon or lens). In research, phase masks with a gradient, step, or both dependences of the phase shift value on the radius vector in combination with an axicon or a lens are actively used [27,28]. There are also sophisticated ways to achieve the same goal [29].
Research in the field of parametrization of the quasi-uniform axial intensity distribution of the Bessel beam is relevant. The purpose of this work was to develop a laser optical system for forming a Bessel beam with a controlled quasi-uniform axial distribution of field intensity on a standard element base with a minimum number of optical elements without significant restrictions on their radiation resistance. Research is being carried out on the possibility of realizing such a laser optical system by apodizing the input aperture of the axicon.

Modeling
The fundamental TEM 00 mode formed by laser resonators of a stable configuration corresponds to a Gaussian beam with a radially symmetric field distribution [30]: where u 00 is the distribution of the complex amplitude of the Gaussian beam in the section with the z coordinate from the waist; r = x 2 + y 2 denotes the distance from the observation point in the beam to the z axis; x, y are the transverse coordinates; A w denotes the amplitude constant or the field amplitude on the axis (r = 0) in the section of the beam waist (z = 0); h(z) and R Φ (z) are the dependences of the caustic and radius of wavefront curvature of the beam in section z; 2h w denotes the waist diameter of the Gaussian beam; z c = πh 2 w M 2 λ is Rayleigh length of the Gaussian beam; λ is laser wavelength; M 2 is the parameter quality of a Gaussian beam.
Let us consider the transformation of a Gaussian beam with field distribution Equation (1) by a binary mask and a refractive axicon with a refractive index n and apex angle 2α 0 ( Figure 1). In wave optics, the effect of an axicon is described by the complex amplitude transmittance [31]: is phase transmittance, p = sin β is the numerical aperture of the axicon, β = arcsin(n cos α 0 ) + α 0 − π 2 is the angle of wavefront inclination to the optical axis after passing through the axicon in the case of the plane homogeneous wave.
If the effect of the axicon has only a phase nature, then the generated Bessel beam of the zero-order is characterized by the length of the Bessel zone z B and the diameter of the core D B (diameter of the central maximum) of the field intensity distribution of the Bessel beam. They are defined by the following expressions [32,33]: Such a Bessel beam also has an uneven distribution of axial field intensity [33,34]: where z AP is the analysis plane after the axicon and P denotes the power of laser radiation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 9 w 2.405 , tan 2 sin Such a Bessel beam also has an uneven distribution of axial field intensity [33,34]: where AP z is the analysis plane after the axicon and P denotes the power of laser radiation. Let us investigate the influence of binary mask parameters and the input Gaussian beam on forming a quasi-rectangular axial intensity distribution of the zero-order Bessel beam. Amplitude transmittance is determined by the finite axicon aperture and the nature of its transmission within the clear aperture. The calculation of the transverse distribution of field intensity in the analysis plane AP z after the axicon is performed using a certain transfer operator. The choice of the transfer operator depends on the distance between the original plane and the analysis plane. The nonparaxial scalar model based on Rayleigh-Sommerfeld theory allows to obtain correct results at sufficiently close distances from the aperture of the optical element at which the diffraction occurs. Using the integral Rayleigh-Sommerfeld transform of the first type, we find the amplitude-phase distribution of the beam field in the plane of analysis after the axicon [34,35]: where u′ is distribution of the complex field amplitude before the axicon, ( ) 2 2 2 2 AP 2 cos L r z r = ρ + + − ρ ϕ− θ denotes oblique length equal to the distance between the point at Let us investigate the influence of binary mask parameters and the input Gaussian beam on forming a quasi-rectangular axial intensity distribution of the zero-order Bessel beam. Amplitude transmittance τ A (r) is determined by the finite axicon aperture and the nature of its transmission within the clear aperture. The calculation of the transverse distribution of field intensity in the analysis plane z AP after the axicon is performed using a certain transfer operator. The choice of the transfer operator depends on the distance between the original plane and the analysis plane. The nonparaxial scalar model based on Rayleigh-Sommerfeld theory allows to obtain correct results at sufficiently close distances from the aperture of the optical element at which the diffraction occurs. Using the integral Rayleigh-Sommerfeld transform of the first type, we find the amplitude-phase distribution of the beam field in the plane of analysis after the axicon [34,35]: where u is distribution of the complex field amplitude before the axicon, L 2 = ρ 2 + r 2 + z 2 AP − 2ρr cos(ϕ − θ) denotes oblique length equal to the distance between the point at which the initial distribution of the complex amplitude of the field u is specified and a point in the analysis plane after the axicon. To calculate the axial intensity of the output field, we use the following expression for L: Appl. Sci. 2020, 10, 7911 4 of 9 When calculating amplitude-phase field distributions using the integral Rayleigh-Sommerfeld transform within the nonparaxial scalar model, it is necessary to consider the efficiency and error of numerical methods when working with rapidly oscillating functions. One of such practical methods is Levin's method [36], which allows one to calculate integrals with complex integrands of strongly oscillating functions. According to Levin's method, the one-dimensional oscillating integral is transformed into an ordinary differential equation which can then be solved by the collocation method [37].
The parameters of the axicon directly determine the distribution of the beam field in the analysis plane after the axicon, but do not affect the general nature of the effect under study. The presence of a binary mask in front of the axicon leads to the fact that the longitudinal dependence of axial intensity of the Bessel beam does not obey Equation (3). The distribution of the axial intensity of the Bessel zone is modified in such a way that the width of the central core and concentric rings remains unchanged in a certain interval, while the length of the Bessel zone decreases. With different parameters of the binary mask transmission function, it is possible to achieve a sharp change in the edge of the axial intensity distribution of the Bessel zone (Figure 2a-d).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 9 which the initial distribution of the complex amplitude of the field u′ is specified and a point in the analysis plane after the axicon. To calculate the axial intensity of the output field, we use the following expression for L : 2 When calculating amplitude-phase field distributions using the integral Rayleigh-Sommerfeld transform within the nonparaxial scalar model, it is necessary to consider the efficiency and error of numerical methods when working with rapidly oscillating functions. One of such practical methods is Levin's method [36], which allows one to calculate integrals with complex integrands of strongly oscillating functions. According to Levin's method, the one-dimensional oscillating integral is transformed into an ordinary differential equation which can then be solved by the collocation method [37].
The parameters of the axicon directly determine the distribution of the beam field in the analysis plane after the axicon, but do not affect the general nature of the effect under study. The presence of a binary mask in front of the axicon leads to the fact that the longitudinal dependence of axial intensity of the Bessel beam does not obey Equation (3). The distribution of the axial intensity of the Bessel zone is modified in such a way that the width of the central core and concentric rings remains unchanged in a certain interval, while the length of the Bessel zone decreases. With different parameters of the binary mask transmission function, it is possible to achieve a sharp change in the edge of the axial intensity distribution of the Bessel zone (Figure 2a-d). The analysis of the axial intensity distributions of the Bessel zone was carried out using the gradient of the scalar function, which shows the direction and rate of its most rapid increase/decrease [38,39]. In  The analysis of the axial intensity distributions of the Bessel zone was carried out using the gradient of the scalar function, which shows the direction and rate of its most rapid increase/decrease [38,39]. In Figure 2, orange plots represent the value of the function gradient for each point in the analysis plane after the z AP axicon. Based on the data, it follows that there was an abrupt change in the speed and direction of the gradient of the function.
The position of the beam focusing plane after the axicon can be calculated using Equation (2)  theory correlate well with the ray model. It is important to note that the length and position of the Bessel zone directly depend on the parameters of the amplitude mask: the left boundary along the propagation direction is determined by the size of the inner aperture of the mask, and the right boundary is determined by the outer aperture of the mask. The effective length of the Bessel zone is defined as the difference between the positions corresponding to the outer and inner boundaries of the amplitude mask.

Experimental Setup
The possibility of implementing the previously proposed method has been confirmed experimentally. The source was a diode-pumped CW single-frequency laser (Cobolt AB, Solna, Sweden) with a wavelength of λ = 659.6 ± 0.3 nm, which also had a spatial single-mode transverse field structure (TEM 00 ) with a quality factor M 2 < 1.1. An afocal optical system was used to expand the laser beam to a diameter of 4 mm (at 1/e 2 width). Then, a quasi-parallel laser beam illuminated the mask and the axicon. The axicon AX2510-A (Thorlabs, Newton, NJ, USA) had the following parameters: aperture diameter 25.4 mm, apex angle 2α 0 = 160 • . The mask (BMSTU, Moscow, Russia) was installed directly, close to the front surface of the axicon. The mask we used was glass plane-parallel plates with a circular aperture 25.4 mm in diameter, on one of the surfaces of which a sufficiently thick layer of chromium (τ < 0.01) was applied. Figure 1 shows the masks used in the research. The parameters of the mask (clear and inner apertures) were designed the same as in modeling and mentioned in the caption of Figure 2. To register the axial intensity distribution after the axicon, an afocal optical system with 3 x magnification was placed in front of a sCMOS camera CS2100M-USB (Thorlabs, Newton, NJ, USA) on the optical axis.

Results and Analysis
The axial field intensity distributions were acquired with a step of 1 mm along the optical axis. The effect of a diaphragmed Gaussian beam on the axial intensity and length of the Bessel zone is shown in Figure 3a-d. The dots denote the experimental data, and the solid line denotes the approximated experimental data. It can be seen that the experimental results are in good agreement with the theoretical calculations presented in Section 2.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 9 positions of the focusing plane according to the ray model. The calculations based on the wave theory correlate well with the ray model. It is important to note that the length and position of the Bessel zone directly depend on the parameters of the amplitude mask: the left boundary along the propagation direction is determined by the size of the inner aperture of the mask, and the right boundary is determined by the outer aperture of the mask. The effective length of the Bessel zone is defined as the difference between the positions corresponding to the outer and inner boundaries of the amplitude mask.

Experimental Setup
The possibility of implementing the previously proposed method has been confirmed experimentally. The source was a diode-pumped CW single-frequency laser (Cobolt AB, Solna, Sweden) with a wavelength of 659.6 0.3 λ = ± nm, which also had a spatial single-mode transverse field structure ( 00 TEM ) with a quality factor was installed directly, close to the front surface of the axicon. The mask we used was glass planeparallel plates with a circular aperture 25.4 mm in diameter, on one of the surfaces of which a sufficiently thick layer of chromium ( 0.01 τ < ) was applied. Figure 1 shows the masks used in the research. The parameters of the mask (clear and inner apertures) were designed the same as in modeling and mentioned in the caption of Figure 2. To register the axial intensity distribution after the axicon, an afocal optical system with 3 х magnification was placed in front of a sCMOS camera CS2100M-USB (Thorlabs, Newton, NJ, USA) on the optical axis.

Results and Analysis
The axial field intensity distributions were acquired with a step of 1 mm along the optical axis. The effect of a diaphragmed Gaussian beam on the axial intensity and length of the Bessel zone is shown in Figure 3a-d. The dots denote the experimental data, and the solid line denotes the approximated experimental data. It can be seen that the experimental results are in good agreement with the theoretical calculations presented in Section 2.   Figure 4. The results of calculating the gradient of the function are presented in yellow lines in Figure 3. The most interesting area from the entire dependence of the axial intensity distribution of the Bessel zone was the leading and trailing edges. When a Gaussian beam hit an axicon without a mask, the gradient of the function had equal rise and fall times, Figure 3a. In the case when the outer or inner part of the Gaussian beam was diaphragmed, only one edge of the Bessel beam field intensity changed, as described in Section 2. With the external diaphragm, the leading edge in the function remained unchanged, and with the internal diaphragm, the trailing edge remained unchanged. In this case, the rate of change of the intensity distribution profile in the studied areas increased by more than 4 times. In Figure 3b and Figure 3c, one can see these jumps in the gradient speed. After installing the mask with annular aperture, two sharp changes in axial intensity were obtained, Figure 3d. The rate of change also increased by more than 4 times. As mentioned earlier, Figure 4 shows the cross sections of the field intensity distribution of the Bessel beam. As a reference distribution, we took the cross section that corresponded to Figure 4a, when a Gaussian beam with a diameter of 4 mm without apodization fell on the axicon. It can be seen that there were features of the zero-order Bessel function in each distribution that are of interest. The absolute value of field intensity in the analysis plane for the cases in Figure 4b-e corresponds to the intensity profile shown in Figure 3a. Apodization of the beam incident on the axicon did not in any way affect the size or shape of the central maximum (the nucleus of the Bessel zone).  For the z AP sections corresponding to the region of abrupt change in axial field intensity (vertical dashed lines in Figure 3a-d), the cross-sections of field intensity are shown in Figure 4. The results of calculating the gradient of the function are presented in yellow lines in Figure 3. The most interesting area from the entire dependence of the axial intensity distribution of the Bessel zone was the leading and trailing edges. When a Gaussian beam hit an axicon without a mask, the gradient of the function had equal rise and fall times, Figure 3a. In the case when the outer or inner part of the Gaussian beam was diaphragmed, only one edge of the Bessel beam field intensity changed, as described in Section 2. With the external diaphragm, the leading edge in the function remained unchanged, and with the internal diaphragm, the trailing edge remained unchanged. In this case, the rate of change of the intensity distribution profile in the studied areas increased by more than 4 times. In Figure 3b,c, one can see these jumps in the gradient speed. After installing the mask with annular aperture, two sharp changes in axial intensity were obtained, Figure 3d. The rate of change also increased by more than 4 times.
(c) (d)  Figure 4. The results of calculating the gradient of the function are presented in yellow lines in Figure 3. The most interesting area from the entire dependence of the axial intensity distribution of the Bessel zone was the leading and trailing edges. When a Gaussian beam hit an axicon without a mask, the gradient of the function had equal rise and fall times, Figure 3a. In the case when the outer or inner part of the Gaussian beam was diaphragmed, only one edge of the Bessel beam field intensity changed, as described in Section 2. With the external diaphragm, the leading edge in the function remained unchanged, and with the internal diaphragm, the trailing edge remained unchanged. In this case, the rate of change of the intensity distribution profile in the studied areas increased by more than 4 times. In Figure 3b and Figure 3c, one can see these jumps in the gradient speed. After installing the mask with annular aperture, two sharp changes in axial intensity were obtained, Figure 3d. The rate of change also increased by more than 4 times. As mentioned earlier, Figure 4 shows the cross sections of the field intensity distribution of the Bessel beam. As a reference distribution, we took the cross section that corresponded to Figure 4a, when a Gaussian beam with a diameter of 4 mm without apodization fell on the axicon. It can be seen that there were features of the zero-order Bessel function in each distribution that are of interest. The absolute value of field intensity in the analysis plane for the cases in Figure 4b-e corresponds to the intensity profile shown in Figure 3a. Apodization of the beam incident on the axicon did not in any way affect the size or shape of the central maximum (the nucleus of the Bessel zone).   The change in relative position of the Bessel zone core in Figures 4 was associated with camera positioning accuracy during the registration of the cross section of the field intensity distribution. However, this fact did not affect the analysis and interpretation of the results. The same can be said about the triple symmetry of the first maximum, the origin of which has nothing to do with the amplitude diaphragm set in front of the axicon, and can be explained by the presence of aberrations in the afocal optical system located after the axicon.

Conclusions
The diffraction of a Gaussian beam on a binary mask and its propagation through an axicon to form a quasi-uniform axial intensity distribution of a zero-order Bessel beam were studied. With an external or internal diaphragm of a Gaussian beam, only one edge of the axial intensity distribution of the Bessel beam field changed. The use of the mask with annular aperture allowed to obtain two sharp edges in the distribution of the axial field intensity. In addition, the internal diaphragm of the input Gaussian beam allowed to exclude the influence of axicon apex angle quality on the formed Bessel beam. Thus, the choice of the parameters of the mask and the axicon allows one to control the length of the Bessel zone of the formed zero-order Bessel beam with a sharp drop in one or both edges of the axial field intensity. This method and the laser optical system that implements it use a commercially available element base and do not have significant restrictions on radiation resistance when working with high-power radiation. The obtained research results will find application in various technologies of photonics and laser technology, such as optical tweezers, visualization of objects, material processing, etc. As mentioned earlier, Figure 4 shows the cross sections of the field intensity distribution of the Bessel beam. As a reference distribution, we took the cross section that corresponded to Figure 4a, when a Gaussian beam with a diameter of 4 mm without apodization fell on the axicon. It can be seen that there were features of the zero-order Bessel function in each distribution that are of interest. The absolute value of field intensity in the analysis plane for the cases in Figure 4b-e corresponds to the intensity profile shown in Figure 3a. Apodization of the beam incident on the axicon did not in any way affect the size or shape of the central maximum (the nucleus of the Bessel zone).
The change in relative position of the Bessel zone core in Figure 4 was associated with camera positioning accuracy during the registration of the cross section of the field intensity distribution. However, this fact did not affect the analysis and interpretation of the results. The same can be said about the triple symmetry of the first maximum, the origin of which has nothing to do with the amplitude diaphragm set in front of the axicon, and can be explained by the presence of aberrations in the afocal optical system located after the axicon.

Conclusions
The diffraction of a Gaussian beam on a binary mask and its propagation through an axicon to form a quasi-uniform axial intensity distribution of a zero-order Bessel beam were studied. With an external or internal diaphragm of a Gaussian beam, only one edge of the axial intensity distribution of the Bessel beam field changed. The use of the mask with annular aperture allowed to obtain two sharp edges in the distribution of the axial field intensity. In addition, the internal diaphragm of the input Gaussian beam allowed to exclude the influence of axicon apex angle quality on the formed Bessel beam. Thus, the choice of the parameters of the mask and the axicon allows one to control the length of the Bessel zone of the formed zero-order Bessel beam with a sharp drop in one or both edges of the axial field intensity. This method and the laser optical system that implements it use a commercially available element base and do not have significant restrictions on radiation resistance when working with high-power radiation. The obtained research results will find application in various technologies of photonics and laser technology, such as optical tweezers, visualization of objects, material processing, etc.