Nonparaxial Propagation of Bessel Correlated Vortex Beams in Free Space

The nonparaxial propagation of partially coherent beams carrying vortices in free space is investigated using the method of decomposition of the incident field into coherent diffraction-free modes. Modified Bessel correlated vortex beams with the wavefront curvature are introduced. Analytical expressions are presented to describe the intensity distribution and the degree of coherence at different distances. The evolution of the intensity distribution during beam propagation for various source parameters is analyzed. The effects of nonparaxiality in the propagation of tightly focused coherent vortex beams are analyzed.


Introduction
Coherent properties of fields must be taken into account in many problems of wave propagation in free space and in inhomogeneous media. This is due to the fact that real sources (laser, LEDs, etc.) generate partially coherent radiation, and purely coherent sources, as a rule, are not implemented in practice. Partially coherent beams are useful in remote sensing, ghost imaging, optical communication, particle trapping, etc.
The theory of propagation of partially coherent waves in free space, as well as in inhomogeneous media, has now been developed quite fully [1][2][3][4]. The conventional Gaussian-Schell-model partially coherent beam was studied mainly until the recent past. Partially coherent sources with specific propagation properties that lead to the formation of highly directional light beams were considered in [5][6][7].
In recent decades, partially coherent vortex beams have aroused great interest due to their unique propagation properties and hidden correlation features [8][9][10][11][12][13][14][15][16][17][18]. It was shown in [9] that the shape of a focused beam can be controlled by changing the initial spatial coherence length. Recently, a variety of partially coherent beams with extraordinary properties has been presented [15][16][17][18]. In [8], a sufficient condition was proposed for devising a genuine correlation function of a partially coherent beam. Based on this, various partially coherent vortex beams with nonconventional correlation functions were introduced. A new class of Laguerre-Christoffel-Darboux sources emitting shape-invariant beams is introduced in [19,20]. Vortex partially coherent beams carrying orbital angular momentum (OAM) are useful in optical manipulation and trapping, optical machining, optical communication, quantum information, and optical microscopy.
Currently, various beams are being studied, ranging from scalar vortex beams to vector vortex beams [16][17][18][21][22][23]. It was found that the polarization properties of vector partially coherent beams change during propagation in free space.
In [24], a new family of partially coherent beams was introduced, the cross-spectral density function (CSD) of which is represented as an incoherent superposition of fully coherent Laguerre-Gauss modes of arbitrary order. Recently, partially coherent vortex beams have been introduced, which are an incoherent superposition of diffraction-free Bessel

Cross-Spectral Density Function
The coherence properties of the field radiated by the source are described by the cross-spectral density (CSD) function W(r 1 , r 2 , 0) [2] which is defined as the second-order correlations of the field at two different points.
The degree of coherence is defined by the CSD function [2]: γ(r 1 , r 2 , 0) = W(r 1 , r 2 , 0) I(r 1 , 0)I(r 2 , 0) where I(r, 0) = W(r, r, 0) is the average intensity at coincident points. Partially coherent vortex beams with various CSD were considered in the past decades. The CSD function of the Gaussian-Schell-model vortex beam in the source plane is given by [16][17][18]: where r and ϕ are the radial and polar angle coordinates, respectively, w 0 is the beam width, r 0 is the coherence length. The Laguerre-Gaussian correlated Schell-model vortex beam was experimentally demonstrated in [52]: where L 0 p [z] denotes the Laguerre polynomial. In [25], the cross-spectral density of a partially coherent diffraction-free vortex beam with a Bessel-mode structure is considered: Below we consider the vortex Bessel-Gauss correlated beams [24]: where ξ = , w is the spot size at the waist of the beam, r 0 is the coherence length. Note that when r 0 → ∞ (fully coherent case) ξ → 0 and when r 0 → 0 (completely incoherent case) ξ → 1 .
The intensity can be calculated from the CSD function at coincident points r 1 = r 2 , i.e., In Figure 1, the intensity distributions of the field with l = 0 as function of the radial distance are shown for the spatially high coherent (ξ = 0.02) and nearly incoherent (ξ = 0.80) cases.
The intensity can be calculated from the CSD function at coincident points = , i.e., In Figure 1, the intensity distributions of the field with l = 0 as function of the radial distance are shown for the spatially high coherent ( = 0.02) and nearly incoherent (ξ = 0.80) cases. In Figure 2, the intensity distributions of the field with l = 1 as function of the radial distance are shown for the spatially coherent and incoherent cases In Figure 2, the intensity distributions of the field with l = 1 as function of the radial distance are shown for the spatially coherent and incoherent cases.
It is seen from the figures that the intensity distributions of the coherent beam are symmetric with the width of the order of w. For an incoherent beam, the intensity distribution is asymmetric with a long tail, the length of which is of the order of coherence length r coh . A characteristic diffraction angle of a coherent beam in paraxial regime is defined by θ d ∼ λ/w, while for an incoherent beam it is defined by the coherence length, i.e., θ d ∼ λ/r coh . Note that for a nearly incoherent beam r coh w. It is seen from the figures that the intensity distributions of the coherent beam are symmetric with the width of the order of w. For an incoherent beam, the intensity distribution is asymmetric with a long tail, the length of which is of the order of coherence length r coh . A characteristic diffraction angle of a coherent beam in paraxial regime is defined by ~ ⁄ , while for an incoherent beam it is defined by the coherence length, i.e., ~ ⁄ . Note that for a nearly incoherent beam ≪ .

Coherent Mode Representation
Using the coherent mode decomposition, the cross-spectral density function of a source can be represented as [2]: where and ( ) are the eigenvalues and eigenfunctions of the Fredholm integral equation of the second kind Note that all the eigenvalues are non-negative, and the eigenfunctions ( ) form an orthogonal set. If ( ) are chosen to be the coherent Laguerre-Gauss modes then the CSD takes a closed form (5) [24]. The eigenvalues in Equation (7) are equal to

Nonparaxial Propagation in Free Space
Below, we consider the nonparaxial evolution of vortex Bessel-Gauss correlated CSD.

Coherent Mode Representation
Using the coherent mode decomposition, the cross-spectral density function of a source can be represented as [2]: where λ n and Ψ n (r) are the eigenvalues and eigenfunctions of the Fredholm integral equation of the second kind Note that all the eigenvalues λ n are non-negative, and the eigenfunctions Ψ n (r) form an orthogonal set.
If Ψ n (r) are chosen to be the coherent Laguerre-Gauss modes then the CSD takes a closed form (5) [24]. The eigenvalues in Equation (7) are equal to

Nonparaxial Propagation in Free Space
Below, we consider the nonparaxial evolution of vortex Bessel-Gauss correlated CSD.
To do this, instead of Laguerre-Gauss modes [24] in Equation (5), we represent the CSD function as an incoherent superposition of coherent orthogonal diffraction-free vortex Bessel modes which are the solutions of Helmholtz equation.
The evolution of the LG modes is determined by the expression where ψ pl (r, ϕ) are the modal solutions and β pl are the propagation constants of the Bessel modes defined from the Helmholtz equation.
The normalized Bessel functions with radial p and azimuthal l indices can be considered as the modal solutions of Helmholtz equation within the effective depth of field [53]: where µ 1 , µ 2 , . . . are the positive zeros of the Bessel function J l (z). It follows from the orthogonality condition for Bessel functions [54].
that these modes satisfy the equation The solutions (12) form a complete set of mutually orthogonal functions in the given interval [0, R 0 ]. Hence, any field in the initial plane z = 0 can be decomposed into these modal solutions.
Solutions (12) are the non-diffractive Bessel beams with radial and azimuthal indices which are propagation-invariant in free space. Note that the non-diffractive scalar Bessel beams in free space were considered in [55], where a narrow annular slit in the screen together with a spherical lens located at a focal distance from it is proposed to create a zero-order Bessel beam. In [56], the decomposition of the field into Bessel modes was used when considering the scattering of light by small particles. In [57], the mode properties of Bessel beams of nonzero order in free space are discussed. Phase optical components are proposed in [58] for the generation of free-space Bessel modes. In [59,60], the amplitude components of the vector non-diffractive beams were obtained as solutions to the vector Helmholtz wave equation. It was shown in [53] that the normalized Bessel beams (12) with radial indices are the solutions of the scalar Helmholtz equation.
The modal coefficients a pn are determined from the integration Using the integral [61] ∞ 0 we find that The evolution of the intensity distribution is given by where The evolution of the CSD has the form For the particular case of the CSD with l = 0 we have where λ n = ξ n ; ξ = Note that the values of λ n decrease with the mode number n, and the number of modes with a noticeable contribution increases significantly with decreasing coherence radius r 0 .
The intensity profile change at propagation is defined by where Analogously, for the intensity of the source with l = 1 we have where λ n = n! (n+1)! ξ n ; ξ = The intensity change with distance is determined by the expression where

Modified Partially Coherent Vortex Bessel-Gauss Beams
Consider a partially coherent Bessel-Gauss beam with the wavefront curvature radius R f in the plane z = 0.
where R f is the wavefront curvature, k 0 = 2π/λ is the wavenumber. Unlike CSD (5), here, we introduced the new parameter-the wavefront curvature radius R f . In contrast to the conventional LG modes, here, we consider generalized Laguerre-Gauss modes with spherical wavefrontss the eigenfunctions The intensity distribution in this case is expressed by where λ n = n! (n+l)! ξ n ; ξ = The evolution of the modified CSD is given by where c pl = a p √ πR 0 J l+1 (µ pl ) .

Propagation of Coherent Vortex Bessel-Gauss Beams
Below, we apply the mode decomposition method to analyze the propagation of coherent vortex BG beams in free space.
Consider the incident beam at z = 0 with Gaussian and Bessel-Gauss (BG) spatial distributions of the intensity: where w 0 is the radius of a Gaussian beam, A 0 = , w B /µ 1 = γ −1 is the effective width of the BG beam.
The modal amplitude coefficients for these incident beams can be calculated analytically. The expressions for modal coefficients have the form: where µ p are the positive zeros of the Bessel functions, and I 0 (z) and I 1 (z) are the modified Bessel functions of the first kind. The evolution of the incident field is determined by the expression → E (r, ϕ, z) = ∑ pl a pl Ψ pl (r, ϕ, 0)e iβ pl z (34) where β pl = k 0 1 − µ pl k 0 R 0 2 1/2 are the propagation constants of the modes (12).

Bessel-Gauss Beam with l = 0
Consider a Bessel-Gauss beam with l = 0 at the source plane z = 0: where A 0 = 2 , w B /µ 1 = γ −1 is the effective width of the BG beam.
In Figure 3, the intensity distributions of Bessel-Gauss beam with l = 0 as function of the radial distance are presented at different propagation distances.
It follows that the diffraction spreading of the BG beam is significantly less than for Gaussian beam, i.e., the BG beams have a sharp radiation pattern. Note the intensity distributions are in good agreement with the results of numerical simulations obtained using the Fresnel diffraction integral [62].

where
= exp ⁄ , / = is the effective width of the BG beam. In Figure 3, the intensity distributions of Bessel-Gauss beam with l = 0 as function of the radial distance are presented at different propagation distances.  It follows that the diffraction spreading of the BG beam is significantly less than for Gaussian beam, i.e., the BG beams have a sharp radiation pattern. Note the intensity distributions are in good agreement with the results of numerical simulations obtained using the Fresnel diffraction integral [62].

Bessel-Gauss Beam with l = 1
Let us now consider a Bessel-Gauss beam with l = 1 in the plane of the source z = 0: In Figure 4, the intensity distributions (a, c, e) and powers (b, d, f) of Bessel-Gauss beam with l = 1 as function of the radial distance are presented at different propagation distances.

Bessel-Gauss Beam with l = 1
Let us now consider a Bessel-Gauss beam with l = 1 in the plane of the source z = 0: where The total power P of the incident beam is normalized, i.e., P = |E(r, 0)| 2 rdrdϕ = 1 In Figure 4, the intensity distributions (a, c, e) and powers (b, d, f) of Bessel-Gauss beam with l = 1 as function of the radial distance are presented at different propagation distances.
It follows from simulations that the beam profile retains its original shape when propagating over long distances.  It follows from simulations that the beam profile retains its original shape when propagating over long distances.

Nonparaxial Propagation and Focusing of a Gaussian Beam
The propagation of a Gaussian beam in free space has long been well studied using various methods, including analytical, asymptotic and numerical approaches to integration. Usually, the problem is solved by evaluation of diffraction integrals. However, numerical simulation of diffraction integrals is often time consuming, and asymptotic methods have been developed for calculations. A hybrid integration method combining numerical integration with asymptotic methods was proposed in [63]. In [64], the accuracy and computational savings of this hybrid technology were examined.

Nonparaxial Propagation and Focusing of a Gaussian Beam
The propagation of a Gaussian beam in free space has long been well studied using various methods, including analytical, asymptotic and numerical approaches to integration. Usually, the problem is solved by evaluation of diffraction integrals. However, numerical simulation of diffraction integrals is often time consuming, and asymptotic methods have been developed for calculations. A hybrid integration method combining numerical integration with asymptotic methods was proposed in [63]. In [64], the accuracy and computational savings of this hybrid technology were examined.
Below, we show that the effects resulting from the diffraction-free mode decomposition method are in good agreement with the known results.
Consider the incident field at z = 0: where A 0 = 2 π 1 w ; k 0 = 2π λ ; R f is the wavefront curvature. Evolution of the field with distance is given by where a p = √ 2 ) . Figure 5 shows the radial intensity distributions I(r, z) = |E(r, z)| 2 of a Gaussian beam with the wavefront curvature radii R f = 1000 µm and R f = 100,000 µm at different propagation distances.
It can be seen that the beam width decreases significantly at a geometrical focusing plane z = R f = 1000 µm. Note that tight focusing takes place if the wavefront curvature radius is not much different from the width of the incident beam. In this case, the nonparaxial effects become significant.
It is of interest to analyze the behavior of the beam profile near the focusing plane taking into account nonparaxial effects. It is known that in the paraxial approximation, the propagation of a Gaussian beam with the preservation of the profile is observed. However, shape-invariant propagation is disrupted due to nonparaxial effects. Figure 6 shows the changes in the intensity of the field with distance along the axial direction (Figure 6a,c) and the intensity distributions of the focused Gaussian beam in the focusing plane in the radial direction (Figure 6b,d).
It follows that the focal planes of the nonparaxial and paraxial beams do not coincide. The focusing plane is shifted in the opposite axial direction compared to the geometric focusing plane if nonparaxiality is taken into account. For an incident Gaussian beam with a width w = 30 µm and a radius of the wavefront curvature R f = 100 µm, we obtain a displacement of the focal plane by 4.7 µm (Figure 6a). For an incident Gaussian beam with a radius of the wavefront curvature R f = 50 µm, a displacement of the focal plane is 4.0 µm (Figure 6c). There is a significant difference in the axial intensity distributions in front of and behind the focus. Before the focus plane, there are significant oscillations in the field intensity. This asymmetry is caused by nonparaxiality. The beam intensity profile in the focus plane does not correspond to the Gaussian profile. It can be seen, that a noticeable sidelobe appears in the profile of the beam (Figure 6d). There is no sidelobe in the paraxial approximation. Note that these effects were also observed with nonparaxial focusing of light beams in a graded-index medium [36,37]. The observed effects may be important in optical trapping and manipulation of nanoparticles. . Figure 5 shows the radial intensity distributions ( , ) = | ( , )| of a Gaussian beam with the wavefront curvature radii Rf = 1000 µm and Rf = 100,000 µm at different propagation distances. It can be seen that the beam width decreases significantly at a geometrical focusing plane z = Rf = 1000 µm. Note that tight focusing takes place if the wavefront curvature radius is not much different from the width of the incident beam. In this case, the nonparaxial effects become significant.
It is of interest to analyze the behavior of the beam profile near the focusing plane taking into account nonparaxial effects. It is known that in the paraxial approximation, the propagation of a Gaussian beam with the preservation of the profile is observed. However, shape-invariant propagation is disrupted due to nonparaxial effects. Figure 6 shows the changes in the intensity of the field with distance along the axial direction (Figure 6a,c) and the intensity distributions of the focused Gaussian beam in the focusing plane in the radial direction (Figure 6b,d).  It follows that the focal planes of the nonparaxial and paraxial beams do not coincide. The focusing plane is shifted in the opposite axial direction compared to the geometric focusing plane if nonparaxiality is taken into account. For an incident Gaussian beam with a width w = 30 µm and a radius of the wavefront curvature Rf = 100 µm, we obtain a displacement of the focal plane by 4.7 µm (Figure 6a). For an incident Gaussian beam with a radius of the wavefront curvature Rf = 50 µm, a displacement of the focal plane is 4.0 µm (Figure 6c). There is a significant difference in the axial intensity distributions in front of and behind the focus. Before the focus plane, there are significant oscillations in the field intensity. This asymmetry is caused by nonparaxiality. The beam intensity profile in the

Conclusions
Thus, the nonparaxial evolution of a closed-form CSD of a Bessel-correlated beam in free space is represented by an incoherent superposition of diffraction-free Bessel modes. Instead of partially coherent beams with Laguerre-Gauss modes [24], here, we have considered a family of diffraction-free Bessel vortex beams. It is shown that the decomposition of arbitrary incident fields into Bessel beams with a truncated aperture is an effective method for analyzing nonparaxial propagation and tight focusing of light in free space.
It was shown in [24] that the cross-spectral density is shape invariant during propagation in free space, due to the fact that the LG modes are shape invariant on paraxial propagation. Our results show that shape-invariant propagation is not observed when nonparaxial effects become significant.
The observed nonparaxial effects (asymmetry in the intensity distribution in the axial direction and the appearance of sidelobes in the transverse field intensity distribution) at a beam focusing in free space were also shown earlier by numerical modeling [64]. Our results show that the beam profile of a partially coherent vortex beam can be shaped by changing its initial radius of wavefront curvature, which is useful for optical trapping.
Here, we have considered the nonparaxial propagation of scalar fields in free space. However, electromagnetic fields have a vector nature, and polarization effects play a significant role in the propagation of partially coherent fields. Therefore, it is of interest to extend the results to vector fields. This is especially important for describing tightly focused beams. The extension of coherence from the scalar to the vector domain can be formulated using generalized two-point Stokes parameters [65]. It is of fundamental and practical interest for studying the coherence-induced polarization effects in free space [66][67][68][69][70][71][72][73][74][75].
Future research may be related to the consideration of the propagation of vector vortex partially coherent and partially polarized beams. Of particular interest is the consideration of structured vortex flows with an orbital angular momentum [76][77][78][79][80][81][82] and the effects of nonparaxiality and depolarization during propagation [83,84]. Other important topics are the coherence and polarization effects in gradient-index medium [85], plasmonic structures [86,87], in turbulent media, etc.
In summary, it is shown that the coherent mode decomposition is an effective method for analyzing the nonparaxial propagation of partially coherent vortex beams in free space. A modified cross-spectral density function corresponding to the family of Bessel-correlated beams is introduced. Explicit analytical expressions for the mode decomposition weights of the CSD function are obtained. The possibility of analyzing nonparaxial propagation and focusing of Bessel-Gauss vortex beams in free space using the mode decomposition method is demonstrated. A noticeable asymmetry in the axial intensity distribution in front of and behind the focus caused by nonparaxiality is shown.
The results may be useful for trapping microparticles where a focused beam spot with a special beam profile is required, and may also be of interest for information transmission, optical imaging and free space optical communication.

Conflicts of Interest:
The author declares no conflict of interest.