Double and Square Bessel–Gaussian Beams

We obtain a transform that relates the standard Bessel–Gaussian (BG) beams with BG beams described by the Bessel function of a half-integer order and quadratic radial dependence in the argument. We also study square vortex BG beams, described by the square of the Bessel function, and the products of two vortex BG beams (double-BG beams), described by a product of two different integer-order Bessel functions. To describe the propagation of these beams in free space, we derive expressions as series of products of three Bessel functions. In addition, a vortex-free power-function BG beam of the mth order is obtained, which upon propagation in free space becomes a finite superposition of similar vortex-free power-function BG beams of the orders from 0 to m. Extending the set of finite-energy vortex beams with an orbital angular momentum is useful in searching for stable light beams for probing the turbulent atmosphere and for wireless optical communications. Such beams can be used in micromachines for controlling the movements of particles simultaneously along several light rings.


Introduction
Among different light beams, there are not many beam families whose propagation in free space is described analytically. Thus, to predict beams' behavior upon propagation, focusing, or other conversions, it is often necessary to use some numerical approaches which can be inexact and cannot explain the physical effects that can appear or disappear when some parameters are changed. Thus, finding a new light beam that can be described analytically is a significant achievement. Among the paraxial vortex light beams that have an analytical description, the most widely known are the Laguerre-Gaussian [1] or Bessel-Gaussian beams [2]. Other examples of such beams are hypergeometric Gaussian beams [3] and modes [4] or generalized Laguerre-Gaussian beams [5].
Starting from the well-known work by F. Gori and co-authors [2], miscellaneous modifications of the Bessel-Gaussian (BG) beams are being investigated with a growing interest. The BG beams can be generated by an axicon, a light modulator, or other elements [6][7][8].
Bessel-Gaussian beams are convenient for probing the turbulent atmosphere since, on one hand, they are of a finite energy like the Gaussian beam and, on the other hand, they manifest quasi-diffraction-free properties due to the Bessel function [9][10][11][12]. As was shown in [13], BG beams are more resistant to the distortions induced by a turbulent atmosphere than the Gaussian beam. The authors of [14] demonstrated that in a Kolmogorov-type turbulent atmosphere, BG beams conserve along a longer distance than the Gaussian beam. The free-space propagation of BG beams in the THz wavelength range was studied in [15]. Bessel-Gaussian beams are used for the manipulation of microparticles [16,17] and for generating entangled pairs of photons in quantum informatics [18,19]. Therefore, extending the set of different types and modifications of BG beams is a relevant problem. In [20], modified BG beams called "frozen waves" were studied. In [21], the propagation of asymmetric BG beams [22,23] in a turbulent atmosphere was considered.
In this work, we obtain modifications of BG beams, such as square BG beams and products of the BG beams. For such beams, we derive explicit analytical expressions that describe their evolution upon propagation in free space. Similar to standard BG beams, the beams we study are not structurally stable, but they are also composed of a set of concentric light rings. The number of rings and the light energy between them change upon the beam's propagation in free space. In the far field, all the beams studied except one are similar to the standard BG beam and contain a single light ring with almost no side lobes. Only one beam, the vortex-free power-function BG beam, contains several rings in the far field, and their number equals that of the beam order.

Bessel-Gaussian Beams and Modulated Bessel-Gaussian Beams
Bessel-Gaussian (BG) beams, first studied in [2], are not structurally stable beams, i.e., their transverse intensity distribution changes when the beams propagate in free space. However, the change is insignificant, and all the intensities keep their shape of concentric circles. The complex amplitude of the BG beam at an arbitrary distance along the optical axis z is where q(z) = 1 + iz/z 0 , z 0 = kw 2 0 /2 is the Rayleigh range, (r, ϕ) are the polar coordinates in the transverse plane, m is the topological charge, J m (x) is the mth-order Bessel function of the first kind, w 0 is the Gaussian beam waist radius, k is the wavenumber of light, c is a dimensionless (possibly complex) parameter affecting the transverse component of the wave vector (c/w 0 = k r ), and k 2 = k 2 r + k 2 z , k z is the longitudinal component of the wave vector. Here and in the following text, we will use a vertical line to separate variables and parameters.
The complex amplitude of a BG beam (1) in the initial plane can be expanded into a series of Laguerre-Gaussian (LG) beams: where the LG beams in the initial plane are described by the following expression: LG n,m (r, ϕ) = exp − r 2 Based on the associated Laguerre polynomial L m n (x). As can be seen from Equation (2), when the parameter c is small (c 1), only the first term in the series is significant, and the BG beam in this case is close to the LG beam.
In addition to standard BG beams, BG beams with quadratic radial dependence (qBG) are also known. They were introduced in [24]. Later, in [25], the transformation of these beams by an optical ABCD system was investigated, and now qBG beams are considered part of the classification of circular optical beams [26]. In the initial plane, the complex amplitude of a qBG beam is defined as As can be seen from Equation (4), the qBG beams depend on a Bessel function of an integer order when m is even and a Bessel function of a half-integer order when m is odd. It can also be seen that in contrast with BG beams, the argument of the Bessel function depends quadratically on the radial coordinate. The complex amplitude of a qBG beam in the Fresnel zone is as follows: The parameters q ± in Equation (4) depend on the propagation distance z and the parameter c: Equations (5) and (6) reveal that when the value c is large enough (c 1), the imaginary part of the factor q + q − in the argument of the Bessel and the exponential functions is small compared to the real part, and the BG becomes approximately propagation-invariant, i.e., its intensity distribution shape is almost conserved, changing only in scale.
Further, we obtain an integral transform that relates BG beams (1) with qBG beams (4). It can be shown that this transform is given by Since the complex source idea then reveals a common background of both beams, described in Equations (1) and (5).

Square Bessel-Gaussian Beams
In our works [27,28], we considered double Laguerre-Gaussian (dLG) beams and square LG beams. Both of these types of beams can be expressed via finite sums of LG beams. In [29], (Chapter X "Orthogonal polynomials", Section 10.12 "Laguerre polynomials"), the following relation between the Laguerre polynomials and the Bessel function is given: Thus, it can be seen that when the radial index tends to infinity, the number of light rings should also tend to infinity, and the radii of these rings are described by the Bessel function. This leads to the qualitative difference between the Laguerre-Gaussian and the Bessel-Gaussian beams-the former have a finite number of rings, and the latter have an infinite number of rings. A natural continuation of the works [27,28] is the investigation of whether similar solutions to the Helmholtz equation can be obtained that describe square BG beams or the products of the BG beams.
Here, we consider square BG beams and show that they can be represented as an infinite sum of BG beams (1). The complex amplitude of a square BG beam (BBG) in the initial plane reads as Here, we prefer to scale the parameters of the initial BG beam ( c → c √ 2 , w 0 → w 0 √ 2 ) but keep the Gaussian factor unchanged.
The complex amplitude of the BBG beam can be obtained by using a generating function for the squares of the Bessel functions [30]: After rather complicated algebra, we reduce the Fresnel transform of the initial complex amplitude (9) to the following series of Bessel functions: Equation (11) indicates that the square BG beams do not conserve their shape upon propagation in free space but are a superposition of a countable number of the products of BG beams of the orders whose sum is equal to 2m, i.e., the initial topological charge. It is interesting to note that in the far field, i.e., when z → ∞ , the series in Equation (11) reduces to the power function r 2m . Thus, in the far field (and in the focus of a spherical lens), the square BG beam has the shape of a single light ring without side lobes.

Product of Two Bessel-Gaussian Beams
Since the series in Equation (11) contains the product of two similar Bessel functions depending on r, it seems quite possible that for a product of two different BG beams, instead of the squared one, we can evaluate the Fresnel transform. Let us introduce the product of two BG beams in the initial plane z = 0: Then, in same way as expansion (11) was derived, the initial field (12) leads to the following solution of the paraxial equation: As can be seen from Equation (13), if n = m and a = b = c, the product of two BG beams reduces to the square BG beam. It can also be seen that for the case n = b = 0, the series in (13) collapses to the only term with ν = 0, and the product of two BG beams reduces to the standard BG beam: Some other particular cases of the dBG beams, as in Equation (13), are also interesting to mention. If n = −m, then the dBG beam is a vortex-free beam: Micromachines 2023, 14, 1029 5 of 11 Its limiting case, when b vanishes, is a product of the vortex-free BG beam by the power function: The beam from Equation (16) can be called a vortex-free power-function BG beam. Tending the parameter b to zero in Equation (15), we obtain nonzero items of the series for −m ≤ ν ≤ 0 only. Then, changing the summation index ν → −ν , we derive the Fresnel transform of the beam from Equation (16): Although the square and products of the BG beams are given by infinite superpositions (Equations (11) and (13)), expression (17) indicates that the power-function vortex-free BG beam in the Fresnel diffraction zone is a finite superposition of similar power-function BG beams of the orders ν from 0 to m.
In the far field (z z 0 ), the argument of the Bessel functions in Equation (17) becomes small. Thus, since J ν (ξ) ≈ (ξ/2) ν /ν! at ξ ≈ 0, the sum in Equation (17) transforms into an ordinary Laguerre polynomial: As result, in the far field, the pBG beam appears as m light rings.

Simulation
In this section, we describe the computational results of the beams from Equations (5), (13) and (17). All distributions are obtained in two ways: by using the numerical Fresnel transform, implemented as a convolution with adopting the discrete fast Fourier transform, and using the theoretical expression. Using the discrete fast Fourier transform is actually equivalent to computing the convolution integrals via Rieman sums, i.e., splitting the integration area into rectangles. All intensity distributions obtained in these two ways are visually indistinguishable, while the phase distributions are different only in low-intensity areas. This confirms the correctness of the Formulae (5), (13) and (17) for the complex amplitudes upon space propagation.
Shown in Figure 1 are the intensity and phase distributions of the modulated Bessel-Gaussian beam (5) in several transverse planes. To obtain a beam with several rings and to prove that it is approximately propagation-invariant, we choose a large value of the scaling factor c = 30. Figure 1 confirms that when the scaling factor c is large enough, the transverse intensity shape almost does not change upon propagation. Figure 2 depicts the intensity and phase distributions of a BG beam and of the corresponding square BG beam (9) (11) in several transverse planes.   Figure 1 confirms that when the scaling factor c is large enough, the transverse intensity shape almost does not change upon propagation. Figure 2 depicts the intensity and phase distributions of a BG beam and of the corresponding square BG beam (9) (11) in several transverse planes.  (a,b,e,f,i,j,m,n) and a square BG beam (c,d,g,h,k,l,o,p), Equations (9) and (11), in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist   Figure 1 confirms that when the scaling factor c is large enough, the transverse intensity shape almost does not change upon propagation. Figure 2 depicts the intensity and phase distributions of a BG beam and of the corresponding square BG beam (9) (11) in several transverse planes.   (9) and (11), in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist radius w 0 = 1 mm, beam orders m = n = 2, scaling factors a = b = 12, propagation distances z = 0 (a-d), z = z 0 /2 (e-h), z = z 0 (i-l), and z = 2z 0 (m-p).
According to Figure 2, the square BG beam is narrower in the initial plane, and its side lobes are suppressed due to the squared Bessel function. A narrower distribution in the initial plane leads to higher space frequencies and thus to higher divergence upon propagation in free space. This is confirmed by Figure 2g,k,o. A notable feature of the squared BG beams is the much smaller dark area in the center of the diffraction pattern. Figure 3 illustrates the intensity and phase distributions of two BG beams with different parameters, as well as of the beams described in (12) and (13), constructed as their product, in several transverse planes. For computation based on Equation (13), the series was bounded by 100 terms. radius w0 = 1 mm, beam orders m = n = 2, scaling factors a = b = 12, propagation distances z = 0 (a-d), z = z0/2 (e-h), z = z0 (i-l), and z = 2z0 (m-p).
According to Figure 2, the square BG beam is narrower in the initial plane, and its side lobes are suppressed due to the squared Bessel function. A narrower distribution in the initial plane leads to higher space frequencies and thus to higher divergence upon propagation in free space. This is confirmed by Figure 2g,k,o. A notable feature of the squared BG beams is the much smaller dark area in the center of the diffraction pattern. Figure 3 illustrates the intensity and phase distributions of two BG beams with different parameters, as well as of the beams described in (12) and (13), constructed as their product, in several transverse planes. For computation based on Equation (13), the series was bounded by 100 terms. Intensity (a,c,e,g,i,k,m,o,q,s,u,w) and phase (b,d,f,h,j,l,n,p,r,t,v,x) distributions of two BG beams with different parameters (a-d,g-j,m-p,s-v), as well as of the beams described in Equations (12) and (13), constructed as their product (e,f,k,l,q,r,w,x), in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist radius w0 = 1 mm, orders of BG beams m = 2 (a,b,g,h,m,n,s,t) and n = 3 (c,d,i,j,o,p,u,v), scaling factors of BG beams a = 8 (a,b,g,h,m,n,s,t) and b = 5 (c,d,i,j,o,p,u,v), propagation distances z = 0 (a-f), z = z0/2 (g-l), z = z0 (m-r), and z = 2z0 (s-x). The grids on the intensity distributions are shown to match the radii of the light rings in different beams.
As can be seen in Figure 3, the intensity distributions of both BG beams in the initial plane have the shape of a single bright ring (Figure 3a,c), but the first beam has a ring with a smaller radius and a pale second ring (Figure 3a). Between these rings, there is a dark, zero-intensity ring, and the multiplication of the complex amplitudes of both beams (two multipliers in Equation (12)) leads to two bright, light rings (Figure 3e), since the thick ring in Figure 3c is "cut" into two rings by the dark ring from Figure 3a. This is a key difference here from the square BG beams, i.e., instead of one dark spot in the center of a circular light distribution, the diffraction pattern contains a dark ring between two light rings.
In the initial phase distributions (Figure 3b,d,f), there are rings with phase jumps of π. The Bessel functions are equal to zero on these rings. However, upon propagation, the Intensity (a,c,e,g,i,k,m,o,q,s,u,w) and phase (b,d,f,h,j,l,n,p,r,t,v,x) distributions of two BG beams with different parameters (a-d,g-j,m-p,s-v), as well as of the beams described in Equations (12) and (13), constructed as their product (e,f,k,l,q,r,w,x), in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist radius w 0 = 1 mm, orders of BG beams m = 2 (a,b,g,h,m,n,s,t) and n = 3 (c,d,i,j,o,p,u,v), scaling factors of BG beams a = 8 (a,b,g,h,m,n,s,t) and b = 5 (c,d,i,j,o,p,u,v), propagation distances z = 0 (a-f), z = z 0 /2 (g-l), z = z 0 (m-r), and z = 2z 0 (s-x). The grids on the intensity distributions are shown to match the radii of the light rings in different beams.
As can be seen in Figure 3, the intensity distributions of both BG beams in the initial plane have the shape of a single bright ring (Figure 3a,c), but the first beam has a ring with a smaller radius and a pale second ring (Figure 3a). Between these rings, there is a dark, zero-intensity ring, and the multiplication of the complex amplitudes of both beams (two multipliers in Equation (12)) leads to two bright, light rings (Figure 3e), since the thick ring in Figure 3c is "cut" into two rings by the dark ring from Figure 3a. This is a key difference here from the square BG beams, i.e., instead of one dark spot in the center of a circular light distribution, the diffraction pattern contains a dark ring between two light rings.
In the initial phase distributions (Figure 3b,d,f), there are rings with phase jumps of π. The Bessel functions are equal to zero on these rings. However, upon propagation, the arguments of the Bessel functions become complex, and the function values are nonzero. Therefore, there are no such phase jumps in Figure 3h,j,l,n,p,r,t,v,x.
Despite the lower topological charge of the first BG beam (m = 2 vs. n = 3), its scaling factor is, vice versa, greater than that of the second beam (a = 8 vs. b = 5). Therefore, upon propagation, it diverges faster, and at a distance of z = z 0 /2, its light ring has a greater diameter than the ring of the second BG beam; it is now this ring that is cut by the minimal-intensity ring of the second beam. Therefore, the beam in Figure 3k,l also has two light rings, as in the initial plane in Figure 3e,f.
Upon propagation into the Fresnel zone and the far field (rows 3 and 4 in Figure 3), the light rings of both BG beams almost do not overlap (rings in Figure 3m,o in the Fresnel zone and rings in Figure 3s,u in the far field). Therefore, after multiplication, these rings are suppressed, and other rings appear: those that are not seen in the initial plane. Thus, the diameter of the outer light ring of the product beam (Figure 3q,w) exceeds the ring diameters of both BG beams. Figure 4 depicts the intensity and phase distributions of the vortex-free power-function BG beam (17) in several transverse planes. arguments of the Bessel functions become complex, and the function values are nonzero. Therefore, there are no such phase jumps in Figure 3h,j,l,n,p,r,t,v,x. Despite the lower topological charge of the first BG beam (m = 2 vs. n = 3), its scaling factor is, vice versa, greater than that of the second beam (a = 8 vs. b = 5). Therefore, upon propagation, it diverges faster, and at a distance of z = z0/2, its light ring has a greater diameter than the ring of the second BG beam; it is now this ring that is cut by the minimal-intensity ring of the second beam. Therefore, the beam in Figure 3k,l also has two light rings, as in the initial plane in Figure 3e,f.
Upon propagation into the Fresnel zone and the far field (rows 3 and 4 in Figure 3), the light rings of both BG beams almost do not overlap (rings in Figure 3m,o in the Fresnel zone and rings in Figure 3s,u in the far field). Therefore, after multiplication, these rings are suppressed, and other rings appear: those that are not seen in the initial plane. Thus, the diameter of the outer light ring of the product beam (Figure 3q,w) exceeds the ring diameters of both BG beams. Figure 4 depicts the intensity and phase distributions of the vortex-free power-function BG beam (17) in several transverse planes.  (17) in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist radius w0 = 1 mm, beam order m = 3, scaling factor a = 15, propagation distances are z = 0 (a), z = z0/2 (b,c), z = z0 (d,e), z = 2z0 (f,g), z = 5z0 (h,i), and z = 10z0 (j,k). Dashed squares (b,d,f,h,j) denote the areas corresponding to shown phase distributions (c,e,g,i,k).
As can be seen in Figure 4, the intensity distribution in the initial plane consists of multiple light rings (there are seven rings in Figure 4a). Upon propagation, only two rings remain; then, in the far field, the number of rings increases to three, which is consistent with the theory that predicts that there should be m rings in the far field. The phase distribution is not shown for the initial plane since it is zero, while in other planes, it is seen to be a rotationally symmetric, and the beam does not contain optical vortices.
Thus, in this section, the numerical simulation confirms our main theoretical ideas. First of all, the simulation shows that all the derived mathematical expressions are correct. Second, the hypothesis of the approximate shape-invariance of the modulated BG beams is confirmed when the scaling factor of the Bessel function is large enough. Third, the simulation confirms that the vortex-free power-function BG beam contains in the far field several light rings whose number equals to the beam order. In addition to confirming the theoretical predictions, the simulation demonstrates a feature of the product of  (17) in several transverse planes for the following computational parameters: wavelength λ = 532 nm, Gaussian beam waist radius w 0 = 1 mm, beam order m = 3, scaling factor a = 15, propagation distances are z = 0 (a), z = z 0 /2 (b,c), z = z 0 (d,e), z = 2z 0 (f,g), z = 5z 0 (h,i), and z = 10z 0 (j,k). Dashed squares (b,d,f,h,j) denote the areas corresponding to shown phase distributions (c,e,g,i,k).
As can be seen in Figure 4, the intensity distribution in the initial plane consists of multiple light rings (there are seven rings in Figure 4a). Upon propagation, only two rings remain; then, in the far field, the number of rings increases to three, which is consistent with the theory that predicts that there should be m rings in the far field. The phase distribution is not shown for the initial plane since it is zero, while in other planes, it is seen to be a rotationally symmetric, and the beam does not contain optical vortices.
Thus, in this section, the numerical simulation confirms our main theoretical ideas. First of all, the simulation shows that all the derived mathematical expressions are correct. Second, the hypothesis of the approximate shape-invariance of the modulated BG beams is confirmed when the scaling factor of the Bessel function is large enough. Third, the simulation confirms that the vortex-free power-function BG beam contains in the far field several light rings whose number equals to the beam order. In addition to confirming the theoretical predictions, the simulation demonstrates a feature of the product of two BG beams that was not predicted by the above theory. In contrast to the conventional BG beams that have one bright ring with many side lobes, the products of two BG beams contain two bright rings in the initial plane, in the Fresnel diffraction area, and in the far field since the bright ring of one BG beam in the product is being cut by an intensity at zero or at the minima of the other BG beam.

Conclusions
In this work, we obtained the following results. An integral transform (7) was derived that relates the standard BG beams, Equation (1), and the modified qBG beams, which are different from the conventional qBG beams, Equations (4) and (5). As can be seen from Equation (7), these beams can be treated as the continuous superposition of the standard BG beams, with the weight function equal to a fractional-order Bessel function. Square BG (BBG) beams (9) were proposed and studied, and their complex amplitude depends on the square of the Bessel function. We obtained the complex amplitude of the BBG beams in the Fresnel diffraction zone in the form of a series of products of three different Bessel functions (11). As a generalization of square BG beams, we also investigated double BG (dBG) beams (12), and their complex amplitudes are proportional to the product of two Bessel functions of different orders and of different scales. The complex amplitudes of such beams in the Fresnel diffraction zone were also represented in the form of a series of products of three different Bessel functions (13). We also considered modified BG beams, whose complex amplitudes are equal to a product of the Bessel function by the power function of the radial variable (16). This set of pBG beams is a subset of vortex-free BG beams. When such a beam of the m-th order propagates in free space, it becomes a superposition of a finite number of similar vortex-free power-function BG beams of the orders from 0 to m. Among the listed modifications of the BG beams, the double BG beams are the most general. They constitute a four-parametric beam family in which two parameters are integers (orders) and two parameters are complex-valued. If the complex parameters become real-valued, they define the beam scales. The square BG beams and the powerfunction BG beams are special cases of the double BG beams. When one order and one scale parameter are both equal to zero, the double BG beams reduce to the conventional BG beams from [2]. Thus, the family of the considered beams is significantly wider than the family of the standard BG beams and has two times more degrees-of-freedom. New varieties of the BG beams considered in this work will be useful for probing the atmosphere, in wireless communications, microparticle manipulation, and in quantum informatics for generating entangled pairs of photons. In micromechanics, these laser beams can be used for controlling the movements of microparticles along circular trajectories [31][32][33]. For instance, the one-dark-spot distribution of the square BG beams (from Figure 2) can be adopted for trapping a non-spherical metallic microscopic object and for rotating it around its center of mass, whereas the two-ring distributions of the double BG beams (from Figure 3) can be used for guiding metallic particles along circular paths since such particles, in contrast to dielectric ones, tend toward dark regions instead of bright regions. In addition, metallic particles tend toward an intensity minimum; thus, the squared BG beams can also allow the rotation of metallic particles (the intensity between the two bright rings in Figure 2g,k,o is nearly two times lower than the intensity on these rings). In optical data transmission, two-ring distributions can be used for redundant information encoding since it is more likely that medium-induced distortions will destroy one light ring than two.