Off-axis vortex beam propagation through classical optical system in terms of Kummer confluent hypergeometric function

The analytical solution for the propagation of the laser beam with optical vortex through the system of lenses is presented. The optical vortex is introduced into the laser beam (described as Gaussian beam) by spiral phase plate. The solution is general as it holds for the optical vortex of any integer topological charge, the off-axis position of the spiral phase plate and any number of lenses. Some intriguing conclusions are discussed. The higher order vortices are unstable and split under small phase or amplitude disturbance. Nevertheless, we have shown that off-axis higher order vortices are stable during the propagation through the set of lenses described in paraxial approximation, which is untypical behavior. The vortex trajectory registered at image plane due to spiral phase plate shift behaves like a rigid body. We have introduced a new factor which in our beam plays the same role as Gouy phase in pure Gaussian beam.


I. INTRODUCTION
The question of propagation of light fields containing optical vortices [1][2][3] is getting attention in many fields of modern optical science nowadays. The fundamental case of such a field is vortex beam -the well-defined, single beam (as for example LG beam [4,5]), carrying the optical vortex of any order. Many different problems concerning propagation of the vortex beam have been considered in the literature so far. Some of them consider optical fields containing the lattice of vortices generated for example by three or more plane or spherical waves interference [6][7][8][9]; or so called composed vortices, i.e. vortices which are generated by two or more overlapping beams [10][11][12][13][14][15][16].
The study on single vortex beam propagation has started with the most basic problem, that is propagation of the fundamental vortex beam in a free space [17][18][19]. Next, the problem of vortex beam propagation through the simple optical element revealing circular symmetry was reported [20,21]. Here, the most important part for our considerations are papers devoted to Gaussian beam propagation through the spiral phase plate (SPP) [22][23][24][25][26][27]. SPP is now one of the most common ways of introducing optical vortex into the laser beam. In more advanced approaches the propagation of vortex beam with broken symmetry or through the system with broken symmetry (like for example diffraction by half-plane [28,29] or a phase step [30,31]) was studied. Most of these asymmetrical cases were studied combining numerical and/or strongly approximated analytical method, especially in case of higher order vortices [32][33][34][35][36][37][38][39][40]. Another highly asymmetrical problem is a vortex beam propagation through the turbid media (e.g. atmosphere) [41][42][43].
In the paper [33] the analysis of the Gaussian beam propagation through the off-axis SPP in Fraunhofer approximation is studied. Authors have focused their attention on the vortex point displacement measured by inspecting the asymmetry in intensity distribution at the far field. In the present paper we analyze asymmetrical optical system with Fresnel diffraction theory, which is more general than Fraunhofer one. The analysis of the offaxis high-order vortices is a difficult task. The integrals become highly complicated and some typical tricks often used in the calculations cannot be applied. The good example is stationary phase method [44], which cannot be used since the phase changes very fast in the vicinity of the vortex point. That is the reason why there are only few publications regarding the exact solutions of asymmetric higher order vortex propagation. In paper [45] the elliptic vortex beam propagation is studied. The paper [46] describes the generation of the higher order vortex beam by discretizing spiral phase plate. In paper [47] the generation of vortex beam through fractional spiral phase plate is studied. In papers [48,49] the propagation through off-axis hologram generating the optical vortices is analyzed, also including the effects of misalignment. In papers [50,51] we have provided a solution for asymmetrical vortex propagation in optical vortex scanning microscope (OVSM) presented schematically in Fig. 1. In this paper we propose more general solution in terms of Kummer confluent hypergeometric function which can be used for a system of arbitrary number of lenses.
In this paper the optical system shown schematically in Fig. 2 is studied. The system represents the object arm of the OVSM shown in Fig. 1 (in experiment, the reference arm is necessary for reconstructing both amplitude and phase of the object beam. Obviously, in analytical and numerical calculations we do not need it). The system is divided into blocks. Each block consists of one or more elements and is represented by its transmittance function (in case of single element) or the product of the transmittance functions (when there are two or more elements in the given block). In our first approach [50] the OVSM was reduced to a single block consisting of three elements, i.e. incident Gaussian beam, SPP and focusing lens considered as a single thin element. The image was calculated at the sample plane (noted as sample in Fig. 1). It should be noticed that the SPP can be moved perpendicularly to the optical axis, which breaks the system symmetry. In result the vortex point moves inside the focused beam, but the range of this movement is highly reduced due to focusing lens. The inset in Fig. 1 shows the exemplary vortex trajectories as registered in our experimental system. In this way the sample can be scanned with the vortex point (i.e. point where the phase is singular). This technique is named the Internal Scanning Method (ISM) [36,38,[50][51][52][53]. In the paper [51] the system built of three blocks was analyzed. The first block consisted of incident Gaussian beam, SPP and focusing lens, the second was just the sample plane, and the third contained a single imaging lens. Here, we extend the analysis to the fully expanded system shown in Fig. 2 Our analysis was performed within the frame of scalar diffraction theory in Fresnel approximation [54]. The Fresnel diffraction integral for the first block shown in Fig. 2 where T v is transmittance function of the SPP, u G is an incident Gaussian beam, T f is a transmittance function of the focusing lens.
m is a topological charge of the optical vortex; it is an integer -positive or negative, w 0 is a beam waist (of the Gaussian beam), f 1 is a focal length of the focusing lens.
Instead of moving the SPP off the optical axis by the distance x c we moved the rest of the system (including the screen) by the same x c , which simplified further calculations. As was shown in [50] the integral Eq. (1) had a solution where The parameters α and β are The main part of this solution is function Kappa K.
We have obtained an interesting result showing that a system having more blocks (with lenses or just planes) can still be described by Kappa function, but with different coefficients y , C (q) . The number of blocks in the system is indicated by superscript (q) as it has been already done in Eq. (3) and Eq. (4). In paper [51] we calculated explicitly y , C (q) up to q = 3, and postulated that it can be done for any (q). In this paper we derived a recurrence formula for the coefficient for any number of blocks (any (q)), which is actually a formal proof of our previous claim.
The function Kappa was denoted previously as G [50,51]. Since that time we have improved and generalized our solution. Now the function G is rewritten in more useful form. In order not to confuse both versions, we denote the present version by capital Kappa

K.
As it was shown in papers [50,51], close to the vortex axis the n = 0 term is sufficient to evaluate the vortex beam. Figure 3 illustrates this fact, but now it is plotted for the four-block system shown in Fig. 2(b). This is a useful result for us, as we analyze the OVSM images at the central part of the beam, where the n = 0 approximation well represents phase and amplitude distribution of our beam. Our goal was to represent the entire object arm of the OVSM as presented in Fig. 2(a). Now the SPP can be separated from the focusing lens which means that focal length for the first block (i.e. SPP block) is f 1 = ∞. Next there is a focusing lens with focal length f 2 = 15 mm followed by the sample plane for which the focal length is also infinite f 3 = ∞.
The third block is the imaging lens with focal length f 4 = 9 mm and the last one is the ocular lens with focal length f 5 = 31 mm, after which we have observation plane (screen).
Unfortunately, for the reasons given later in the manuscript the full system could not be analyzed with Kappa function. The SPP cannot be separated from the focusing lens. Therefore, we have to switch to the system shown in Fig. 2(b). In this paper we will analyze this system showing the efficiency of our formulas. As can be noticed from Fig. 2(b), in our model a sample plane is treated as a separate block. In this way we are able to enhance our calculations in order to analyze the influence of a simple phase object on the vortex beam. Thus, the system shown in Fig. 2(b) prepares us for this next step.
Certainly, all the blocks may have different focal lengths and positions than the ones shown in Fig. 2 and the Kappa function will still work. However, the important thing is that the first element must be the block containing the SPP, focusing lens and Gaussian beam.
The present paper is organized as follows. In section II the extension of our formulas as well as their closed form for the system consisting of any number of lenses is discussed. In section III we discussed the role of coefficient A. In section IV we show the efficiency of our formulas by discussing the effect of vortex trajectory rotation. Section V concludes the paper.

II. ANALYTICAL PART
In this part our previous results will be extended. In order to build more relevant model of the OVSM we have derived the explicit formula for coefficients Ξ q , y , C (q) up to q = 4, just to analyze the OVSM system shown in Fig. 2 where the meaning of z 4 and f 4 can be read from Fig. 2), d 1 ,, d 4 are lens thicknesses. In case of a simple plane we put d = 0.
Following this path a general iterative formulas for coefficient of any (q) index can be derived. The formulas have a form The derivation of the above is presented in Appendix A.
The sums in Kappa function are convergent provided that a series of conditions are hold Very similar formulas can be derived for negative vortex charge. In that case we can use the same Kappa function but with some multiplying expression and B y coefficient multiplied The derivation of this formula is presented in Appendix B.
We have also derived a closed formula, but with a special function, which is given below.
The derivation is discussed in the Appendix C.
where 1 F 1 is the Kummer confluent hypergeometric function [55]. The calculations using closed formula are much faster, but as we will show both formulas Eq. (3) and Eq. (11) are helpful in understanding the vortex beam propagation through our system.
The sum of two B coefficients in Eq. (8b)-(8c) is a complex expression which can be easily decomposed into real and imaginary part. From Eq. (8b)-(8c) we got for real part And for imaginary part The form of coefficients ξ For example for q = 4 we get Formulas in Eq. (12) and Eq. (13) show that y-coordinate is a member of imaginary part and x-coordinate is a member of real part of the B x + iB y expression. This leads us to simple formulas for vortex point trajectory. The vortex point (zero amplitude point) is at the place where both imaginary and real part are equal to zero [1][2][3].
From the above formula we may conclude that the vortex trajectory as a function of x c , for the given z is a straight line. We may also find a plane where the vortex trajectory is perpendicular to the SPP shift. Using Eq. (16b) the condition is In the paper [50] we have formulated the hypothesis that the higher order vortices (m > 1) do not split even when x c = 0. The higher order vortices are classified as structurally unstable [56], i.e. they are supposed to split into single order vortices even under small phase or amplitude perturbation [57][58][59]. Nevertheless, the formulas in Eq. (11) proof the stability hypothesis in an explicit form. The factors in front of the Kummer confluent hypergeometric function are just a vortex term.
The place where real and imaginary part of B (q) x +iB (q) y is zero indicates the position of the vortex point of the m-th order. Since the whole term is at power m, the m-order vortex does not split for any x c . The result is very interesting, which is illustrated in Fig. 4. the m-th order vortex into a set of single vortices. But our case is untypical and the higher order vortex does not split; c) when performing more complex operation as described below, we will produce even less symmetrical phase distribution, but the higher order vortex is still stable.
The situation would not be strange if the higher order optical vortex were stable for a given non-zero value x c . However, changing the x c continuously breaks the circular symmetry of the phase distribution at the SPP plane with no harm for the higher order OV stability.
We could also start our analysis at the sample plane, where both phase and amplitude distribution symmetry is broken (Fig. 5). Nevertheless, the propagation of this input beam through any number of lenses will not split the higher order vortices.
We can go even further. Adding the term inside the bracket in Eq. (18) and multiplying B x or B y coefficients by any number, but in such a way that the coefficients by x (q) remain real (or become imaginary) and coefficient by y (q) remain imaginary (or become real) may change the position of the vortex point and its phase distribution (Fig. 5(c)) but still the higher order vortices will be stable while propagating through our system. However, changing other parameters, like for example binomial factor by coefficients B in Eq. (5) splits the higher order vortices immediately (Fig. 6).
To summarize this part: We have shown that the higher order vortices when propagating through the set of classical lenses described in paraxial approximation will not split regardless of asymmetry introduced by the off axis position of the SPP. Moreover, we may perform some more symmetry breaking operation (as shown in Fig. 4(c)) in Eq. (8c)-(8d), provided that the coefficients by x remain real (or imaginary) and by y imaginary (or real). So, the conclusion is that the very basic optical system does not split the higher order vortices even if the input phase and amplitude distribution is highly non-symmetrical. In many cases such an unusual stability results from deeper physical rules. There is a question if this is also a case here. So far we have no answer to this question. What we can learn now is that classical optical system is somewhat special, at least when being described in paraxial approximation and illuminated by Gaussian beam with the vortex beam introduced by SPP. In the next section the special role of the coefficient A will be studied.

III. COEFFICIENT A
In our further discussion we will refer to two specific examples of the OVSM models. One with the separated SPP plate and focusing lens ( Fig. 2(a)) and the second one with SPP plate and focusing lens working as single thin element ( Fig. 2(b)).
The coefficient A (q) has relatively simple form. It depends neither on the x and y coordinates nor on SPP shift x c . However, it plays a crucial role in vortex beam phase evolution.
Unfortunately it cannot be totally taken outside the first sum in Kappa function. On the contrary, in Eq. (4) it can be entirely assimilated inside the first sum, but in the present form some mathematical aspects can be noticed in a more clear way.
The first important point is that the coefficient A is responsible for breaking the Kappa function convergence. In the OVSM this happens when the SPP is separated from the focusing lens, and coefficient A has a singular point ( Fig. 7(a,b)). At first we will analyze the system with two blocks, the SPP and the focusing lens (sample plane is an observation plane in this case), so we need the second order coefficient A (2) . The coefficient A (q) contributes to the Kappa function as 1/A (q) , so the term in front of the first sum takes form (for q = 2).
There exist a range of such positive z 2 for which the conditions in Eq. (9) fail. In particular there is a z 2 value that the γ 2 equals zero and the whole term in Eq. (19) becomes zero.
Thus, the whole Kappa function is equal to zero, which shows that it does not reflect the true behavior of the focused vortex beam. Moreover, the Kappa function encounters a π-jump in this point. This z 2 value can be computed from the formula Fig. 7(a,b) illustrates the problem. As we can see for the z 2 = 15.58 mm (calculated from Eq. (20)) the 1/A (2) term equals zero.
When the SPP and focusing lens are joined, we only need a coefficient A (1) which has a simple form, free of our problem (at least for the OVSM optical system). This is illustrated in Fig. 7(c,d). The higher order terms A (q) (q > 2) behave in a similar way. When the SPP and focusing lens work as a single element, they meet the conditions in Eq. (9), for a reasonable OVSM configuration. When q > 2 the set of conditions in Eq. (9) is more complicated. The detailed study of this problem would explode the volume of this paper, so we will not follow this path. Important thing is that when SPP and focusing lens work separately, condition in Eq. (9) fails for any q > 1 in case of reasonable OVSM configurations. From practical point of view it is enough to check if the A (q) factor behaves like in Fig. 7(c,d), which is a case in the system shown in Fig. 2(b). It is worth noticing here that when using the ISM all z (q) parameters are fixed, so the test for conditions in Eq. (9) is not difficult. When we move any element of our system along the optical axis the conditions in Eq. (9) have to be checked for the whole range of z (q) coordinates.
For the reason explained above we cannot use our formulas to the system shown in Fig. 2(a). Instead, we have to limit our study to the system shown in Fig. 2(b), when both SPP and focusing lens work together. Certainly for the forbidden area the numerical modeling of our optical system is still possible and effective.
We can also find A (q) coefficient inside the first sum in Kappa function. Now, we multiply this terms by 1/A, located in front of the first sum, so we have 2n+2 ; for q even (21b) The first term (for n = 0) has the largest influence on the phase and amplitude of the Kappa function The next terms rapidly drop in their values. The plot in Fig. 8 is done for q = 4, but it represents the typical curve for expressions in Eq. (21a)-(21b) for any q, provided that we avoid the forbidden area defined by conditions (9). As we can see the part for n = 0 strongly dominates over the part for n = 1. The next terms (for n = 2, 3, ...) are invisible in the figure scale. This domination is particularly strong at the beam center, when x (q) and y (q) coordinates are small (Fig. 3). When the x (q) and y (q) become larger, the coefficient grows rapidly with increasing n and things become more complicated. focusing lens themselves. This is illustrated in Fig. 9(a). If A (1) = 1 and x c = 0, there is no phase rotation when changing the position of the focusing lens along the z-axis. We can conclude that in the OVSM system the coefficient A (q) plays, in some respect, the similar role as the Gouy phase in the Gaussian beam. When the x c = 0 things become more complicated. The off-axis position of the SPP plate introduces a phase value dependence on x c being inside the coefficient B x (Fig. 9(b)).
The other coefficients (but B x , B y ) play minor role. The coefficient B x and B y are responsible for the vortex trajectory evolution as a function of SPP shift x c . This dependence is linear, but as has been already shown, the direction of vortex trajectory becomes perpendicular to SPP trajectory when the condition in Eq. (17) holds. In the previous papers [50,51] this fact was proved for the small x c . Having the formulas in Eq. (11) we can conclude that they hold for any x c . In paper [53] precise experiments were reported which confirm the theoretical results.
The coefficient Ξ q collects all constant factors. It must be observed when we analyze the rotation of the beam phase while moving the first block in the optical system along the z-axis. When the first element moves away or toward the laser source, the phase of the incident Gaussian beam changes. This incident phase is a part of coefficient Ξ q .
The coefficient C (q) multiplies the first sum in function Kappa. Due to the (x 2 q + y 2 q ) factor, the C coefficient is responsible for the equiphase lines curvature (Fig. 10). When C (q) = 1 the equiphase lines are straight (for x c = 0). Since the C (q) coefficient contributes to the Kappa function as exp(x 2 c α + IΛ), where Λ is an imaginary part of the C coefficient, and typically α Λ its contribution to the amplitude is small, especially for small x c .

IV. THE ICE-SKATER EFFECT
To see the effectiveness of our formulas we study the correlation between the focused beam radius and vortex trajectory rotation (see the inset in Fig. 1). The hypothesis was that the vortex trajectory behaves like a rigid body. Consequently we assumed that in the closed system of a focused vortex beam the angular momentum L is conserved, thus Regardless of the assumed model (point mass, flat disk, cylinder, etc.), the moment of inertia I is always proportional to the squared distance r from the rotation axis, i.e. I ≈ r 2 . Therefore, the rotational speed ω of the vortex trajectory shall satisfy ω ≈ 1/r 2 .
In our case, ω is the first derivative of the trajectory inclination angle with respect to z (ω = dθ/dz) and r is the radius of the vortex bright ring, i.e. the distance between points of zero and maximal intensity (e.g. r = 2.4 µm for the vortex analyzed in Fig. 3). Both the radius of the converging vortex beam and the rotational speed of the trajectory depend on the axial distance z. Obviously, as z approaches the focal point, the vortex radius decreases whereas the rotations speed up. As the beam focusses, the radius reaches its minimum and the speed is maximal. This is a direct analogy to the ice skater pulling their arms in for a faster spin. Such a rigid body mechanics approach to the vortex beam has been already studied by Bekshaev et. al. [11], for free propagating Gaussian beam with optical vortex.
The authors concluded that it is related to the Gouy phase dynamics. We have stated that the coefficient A plays the role of the Gouy phase in our case. Indeed, from Eq. (12)- (13) we may calculate the angle of the trajectory inclination as For q = 1 we get Calculating the angle of the equiphase line for the A coefficient we get exactly the same formula.
Here we present the results calculated at the exit of the setup (camera plane) shown in other hand, computing the vortex trajectory rotation speed ω is much easier. In paper [38] the formula for ω was given as: In fact the formula in Eq. (24) was derived for the vortex trajectory rotation at the sample plane for the approximated linear case (n = 0). But the new formulas in Eq. (11) allow to extent their applicability for general case. The classical imaging preserves the vortex trajectory orientation as was shown experimentally in [51,52]. So we expect that the angle rotation at the image plane is the same as in the sample plane.
The obtained relation between the rotational speed and the defocusing is presented in Fig. 11(b), red. Applying the above rigid-body reasoning, ω was compared with the inverse square of the vortex radius ( Fig. 11(b), black). There is a clear agreement between the two curves supporting the hypothesis of rigid-like behavior of vortex trajectory.

V. CONCLUSIONS
There is a growing interest in exact theory describing the propagation of the optical vortex beams in optical system, in particular in a system with broken symmetry. The reasons are both understanding the physics of electromagnetic waves, and practical applications. Some examples concerning both science or applications can be found in [11,14,24,33,52,60,61], but the list of papers is much longer. In this paper we have enhanced the results presented in our former works [50,51]. We have derived the coefficient for Kappa function for any numbers of elements. We have also enhanced the formulas for the vortex trajectory. Moreover, we have found a closed form of our solution in Eq. (11) using the special function. This new form proved that in classical optical system, the higher order vortices do not split even if the circular symmetry is broken. This is true under the conditions of paraxial approximation.
This result is truly surprising. It is hard to check it experimentally. We cannot use an optical system in paraxial approximation. Moreover, any real system introduces errors which break the beam symmetry in a way different than allowed by our theory. This immediately splits the higher order vortices forming a constellation of the first order ones. Nevertheless, the experiment reported in [50] suggested the mass center of such a constellation moves as ideal higher order vortex.
The study on A (q) coefficient have shown that it is responsible for breaking the conditions in Eq. (9). In result we cannot analyze the OVSM system with separated SPP plate and focusing lens. The A (q) coefficient has an important role in the rotation of vortex beam phase, which is similar to the role of the Gouy phase in Gaussian beam. The second important factor is the SPP shift x c . When we are close to vortex core we can limit the sum range (in Kappa function) to n = 0, which is also the result of the A (q) coefficient (Fig. 8).
We have shown that the vortex trajectory behaves like a rigid body. The range of the trajectory shrinks when the beam converge, so its rotation speed increases according to the rules describing the rigid body behavior. The vortex trajectory rotation is strictly related to the phase rotation of the A (q) coefficient. In [11] the same relation was shown for free propagation of axial Gaussian vortex beam and the Gouy phase.
Finally, we have to prove the implication We can write a right side of the implication as After plugging the definition of u j into the integral to obtain In this section we will prove the formula in Eq. (10).