Close relationship between Bessel-Gaussian and conical refraction beams.

We demonstrate that the conical refraction of the input elegant Laguerre-Gaussian beams can be effectively described through generalized Bessel-Gaussian light beams. We performed numerical simulations and show good agreement between the exact solution and our proposed Bessel-Gaussian approximation model. Physical clarity of the proposed model has allowed us to explain the transition of the classical double-ring pattern of conical refraction in the Lloyd plane into a multi-ring one and predict new phenomenon such as the Raman spot shift and dependence of the conical refraction ring radius on the value of the orbital angular momentum.


Introduction
Internal conical refraction (CR) is a phenomenon, which was predicted by Hamilton in 1832. It can be observed for the light propagating along the optical axis of a biaxial crystal. In this case, a narrow beam evolves as a hollow double-walled cone and emerges as a cylinder of light behind the exit facet of the crystal. A strong interest in the CR investigations over the last 20 years was driven by the two significant factors. Firstly, as the 21st century began, the Belsky-Khapalyuk-Berry paraxial theory of CR [1,2] has been finally formulated, which laid a solid theoretical foundation for subsequent progress in the experimental field. Secondly, developing the technical capabilities of modern laboratories made it possible to produce a substantial number of the CR practical applications for optical communication [3], ultra-efficient CR lasing [4][5][6], polarimetry [7], high-harmonic generation [8], particle trapping and optical manipulation [9,10], and for singular beams generation and annihilation [11]. See also recent review [12]. However, utilization of vortex input beams for CR now is one of the most novel and intriguing phenomena in the CR field, because it can significantly change the familiar CR patterns [13]. New specific properties of CR have already been demonstrated, such as the formation of a multi-ring image in the Lloyd's plane [13]. In this sense, studies of CR with elegant Laguerre-Gaussian beams, which are the complete set of solutions to the paraxial wave equation in free space, can be very fruitful.
However, it must be said that expressions describing the spatial evolution of the CR field have an integral structure, which not only complicates numerical calculations but leads to a lack of physical clarity. Thus, the above-mentioned paraxial model of conical refraction contains complete information about CR features, but due to its integral structure, it is extremely unclear from what mechanism certain specific properties of the studied CR phenomenon are formed. This problem can be solved by describing CR through a superposition of beams whose spatial distribution would be determined by a simple and understandable law. With this approach, many properties, such as the multi-ring intensity distribution in the focal plane, Poggendorff's rings, and Raman spots could arise "naturally" as a consequence of the propagation properties of such superposition beams.
In this work, we demonstrate the analogy between CR of the elegant Laguerre-Gaussian beams and Bessel-Gaussian beams in precise quantitative terms. We show that CR of such beams under the condition R 0 >>w 0 , where R 0 is the radius of the CR ring beyond the crystal and w 0 is a beam width, can be effectively described in terms of generalized Bessel-Gaussian beams. Also, we determine the parameters of the generalized Bessel-Gaussian beam from comparing it with the corresponding CR beam and finally analyze the correspondence between the exact solution and the proposed Bessel-Gauss model for the most important CR patterns. This approach, finally, allowed us to predict the Raman spot shift and the dependence of the CR ring radius R 0 from the value of orbital angular momentum, which leads to the necessary renormalization of the proposed BG model parameters.

Theory
Let us first consider the elegant Laguerre-Gaussian beams [14]. One can write the electric field vector E at its waist (z = 0) in the form: where (r, ϕ) are the polar coordinates in the transverse plane, w 0 is the beam width in the focal plane, e in = [d L , d R ] is the polarization vector, where we choose left and right circular polarizations as basis functions and L n (x) is the Laguerre polynomial of order n and index [15]. The dual-cone model of CR is used for a description of the CR beam evolution [16,17]. Under this model, the electric field vector behind the biaxial crystal can be represented as a sum of two CR cones E = C (+) +C (−) as shown in Figs. 1(a) and (b): In these equations (ρ, ϕ, ξ) are the cylindrical coordinates normalized to the beam waist and Rayleigh length of the incident beam, i.e. ξ = z/z R and ρ=r/w 0 , ρ 0 =R 0 /w 0 is a normalized radius of the CR ring beyond the crystal, which plays a central role in the determination of the resulting CR intensity distribution, J +m is the Bessel function of the first kind of order +m, m=0, ±1 is an integer number. It can be seen from Eqs. (2) that, for the input field with the orbital angular momentum (OAM) of , the CR beam comprises two extra components with the OAM of ±1, depending on the input polarization. This phenomenon was first studied by Peet [13] and then properly discussed in [18], where the simplest Laguerre-Gaussian optical vortices were considered. It is worth noting, that two components C (+) and C (−) have a clear physical interpretation as two cones that diverge and converge behind the exit facet of the CR crystal, correspondingly [16,17]. Axial intensity distribution of the individual cone C (+) is presented in Fig. 1(b) for n=0 and =0 (input Gaussian beam).
In this paper, we suppose that if condition ρ 0 >> 1 is met, then CR cones C (±) of arbitrary elegant Laguerre-Gaussian modes can be approximated by a superposition of Gaussian beams, whose centers are placed on a circle of radius R 0 in the focal plane and whose wave vectors lie on a cone of semi-aperture ±α, as shown in Figs. 1(b) and (c) for C (+) . We are going to show that (c) Geometry of the approximation Gaussian beam, where R 0 is a displacement vector of the beam center in the focal plane (z=0), which varies by changing the value of the angle θ, k 0 is a vacuum wave vector, whose tilt with respect to the z direction amounts to the angle α, k ⊥ is a transverse wave vector, which is the projection of k 0 onto the x-y plane. this approach gives a simple propagation law of CR cones C (±) for any values of coordinates (r, z). Furthermore, this approach, as we will see, elegantly explains many CR features, such as multi-ring Lloyd's pattern, the Poggendorff's rings, and the Raman spots, the transverse intensity distributions for which are shown in the lower part of Fig. 1(a).
The electric field vector of the approximating Gaussian beam, thus, in our paraxial treatment has the form: where, k 0 is a vacuum wavenumber, k ⊥ = k 0 ·α[cos(θ), sin(θ)] is a transverse wave vector that appears due to the tilt of the Gaussian beam direction by an angle α, R 0 = R 0 [cos(θ), sin(θ)] is a displacement vector of the beam center in the focal plane (z=0), W 0 is a Gaussian beam width in the focal plane, W(z) = W 0 1 + iz/k 0 W 2 0 is responsible for diffraction broadening and wavefront changing during the spatial propagation. The displacement vector R 0 and the transverse wave vector k ⊥ are the functions of θ. Furthermore, it must be pointed out that in the last term of the Eq. (5) we neglect phase contribution exp(ik 0 z), which is obtained from exp(ik 0 cos(α)z) under the assumption α<<1. The vector A(θ) defines polarization dependence, which is common to the CR beam: where, A 0 is a complex amplitude. We need to integrate over the θ angle from 0 to 2π in (5) to obtain a total field: where we derive the same polarization pattern as in formula (2) and define beam components with different OAM as: In Eq. (8) we use normalized coordinates (ρ, ξ) and introduce new definitions, as follows: a 0 =αkw 0 is a normalized tilt parameter of CR cones, σ(ξ) = σ 2 0 + iξ is a spatial propagation function by analogy to W(z), where σ 0 =W 0 /w 0 is defined as a ratio of the approximating Gaussian beam's width (5) and width of the input elegant Laguerre-Gaussian beam (1), ρ(ξ)=ρ 0 +a 0 ξ is a ring radius for arbitrary ξ value, I +m (x) is the modified Bessel function of the first kind of order +m. To obtain formulae for the second CR cone C (−) , one needs to replace the Eq. (8) by its complex conjugate and substitute ξ → -ξ.
The Eq. (8) is known as a generalized Bessel-Gaussian (BG) beam of order +m. BG beams were introduced in [19] to overcome problems relating to an infinite amount of energy required for Bessel beam formation. Generalized BG beams with ρ 0 0 and a 0 0 for the first time were considered in [20]. However, as it turns out, many of the CR features are typical for generalized BG beams, such as ring distribution with a radius ρ 0 in the focal plane, which is forming axial spike during the spatial propagation along the ξ-axis, that corresponds to a Raman spot. Furthermore, for a 0 0 generalized BG beams have a nonzero projection of k 0 onto the transverse coordinate plane, which leads to the interference picture in the focal plane, caused by the intersection of two BG beams with an opposite value of a 0 , that corresponds to two CR cones C (+) and C (−) . It is worth mentioning, that the idea of the relationship between CR beam and BG beam was introduced in [21] by Peet. However, he considered a modified BG beam, which is a special case of a generalized BG beam for a 0 =0, under the assumption ρ 0 ≈0. . . 2. As we show further, for the case of ρ 0 >>1, one needs to use generalized BG beams with a 0 0.
To determine the parameters of the BG beams, we will consider the CR field in the focal plane, under the condition ρ 0 >>1, which corresponds to a well-defined CR ring. In this case, we can replace the Bessel function under the integral (3) with its asymptotic expansion and obtain a formula for cone propagation, that is valid far away from the propagation axis, i.e. for ρ>>0: where we introduce f (κ) as a one dimensional Fourier transform of the function F(x) that is, for CR integrals, equal to: where θ(κ) is the Heaviside step function. From (9)-(11), one can immediately see that for ρ>>0, CR cone components with different index m are equal to each other. It should be noted that the lower limit of integration in Eq. (10) can be extended to -∞. This is to demonstrate the presence of converging and diverging waves from the ξ-axis with κ>0 and κ<0, respectively. From formula (11), it is pretty obvious that cone S (+) comprises only converge waves, because of the Heaviside function θ(κ). By this approach, the Hankel transform for the components of the BG beam, which are given by formula (8), will be as follows: where we define µ=ρ 0 +ia 0 σ 0 . Using the same procedure as in (9)- (11), we can obtain from the Eq. (13) the following expression for a one-dimensional Fourier transform: where we use the subscripts of the functions f (κ) to draw attention to the fact that we have two different Fourier transforms. In these Fourier transforms, "C" in (11) means that we are dealing with function S m (+) from Eq. (3), and "U" in (14) is related to the function U m from Eq. (12). Thus, to find three parameters of the proposed BG model: a 0 , σ 0 , and A 0 , one needs to approximate the one-dimensional Fourier transform (11) with the Gaussian function (14). This approach in the area of classical optics was applied in [22], where the replacement of the Gaussian and power function with the shifted Gaussian function was necessary to demonstrate the asymptotic equivalence of the BG beams to elegant Laguerre-Gaussian beams with the index n>>1. This is the main reason why the elegant modes fit better to the stated problem than their 'standard' counterparts. Thus, using the method proposed in Appendix B of [22], we can obtain the following formula: where it was assumed n>>| |>>1. Comparing (15) with (14) and (11) one can get all parameters of the BG model: As will be shown later, we propose a method for determining the model parameters that will be valid for arbitrary values of n and , and under assumption n>>| |>>1 will give directly the expressions (15) and (16). Of course, the exact equality between expressions (11) and (15) will be valid only for the case n>>| |>>1, but for arbitrary n and , we can understand this substitution as an effective description of the CR cones by the generalized BG beams with a tilt a 0 and a beam width σ 0 in the normalized coordinates. This procedure is similar to introducing the M 2 parameter for the effective description of the multimode laser radiation. Thus, we suggest the following definition for the BG model parameters, where we use the original one-dimensional Fourier transform of the CR cones (11): This definition is rigorous for the Gaussian distribution (15). Substituting the Fourier spectrum of the CR beam (11) to (17), after simple algebra, one can obtain: where s=2n+| | and Γ(x) is a Gamma function. For s=0: a 0 ≈ 0.886 and σ 0 2 ≈2.33. Also, with the growth of s, a 0 grows as a root of s, and σ 0 2 decreases to 2, which corresponds to the asymptotic behavior of the parameters (16). A comparison of the modulus of the normalized Fourier spectrum f (κ) from (11) and (15) with parameters (18) is shown in Fig. 2. In this case, the complex amplitude A 0 , which determines the normalization of the BG beams relative to the CR cones, is naturally to obtain from a comparison between the maximum intensity of the beams at the ρ=ρ 0 and ξ=0, using the approximation ρ 0 >> 1 and the Eqs. (9) and (10). Thus, we get: As a result, all the model parameters for arbitrary values of n and have been obtained. In the next section, we analyze the behaviour of the proposed BG model for certain space regions, such as the Lloyd plane with a characteristic multi-ring distribution, the Poggendorff's rings, and the Raman spots.

Lloyd's plane and transitional region
Furthermore, the proposed BG model elegantly explains the transition of the classical double-ring pattern of CR in the Lloyd plane into a multi-ring one reported earlier for the simplest optical vortex with >0 [13,18]. As mentioned above, the electric field behind the output facet of the crystal similarly consists of two BG beams converging and diverging, correspondingly, i.e.

E(ρ,ϕ,ξ)=U(ρ,ϕ,ξ)+U(ρ,ϕ,-ξ)*.
We will assume that the light is unpolarized since the exact CR beam and BG approximation have the same polarization feathers and we will be more interested in the space propagation properties of such beams. Since the transverse wave vector k ⊥ has the opposite sign for converging and diverging BG beam, we will observe an interference pattern in the focal plane where they intersect. Its intensity I = E·E* which, in the approximation ρ 0 >>1, will be, up to a constant, equal to: where the Gaussian function determines the width of the intensity profile, the tilt parameter a 0 is inversely proportional to the period of the interference, and the phase determines the shift of the interference picture relative to the point ρ = ρ 0 . It was previously discussed that the width σ 0 weakly depends on the parameters and n, but a 0 grows as a root of 2n+| |. Due to this fact, the interference period with the growth of Laguerre-Gaussian mode indices will decrease, and, with a fixed area of interference, we will observe an increase in the number of dark and light rings, as shown in Figs. 3(a)-(c). Since the phase depends on , the maximum of the interference pattern will be either to the right of the point ρ = ρ 0 for even , or to the left for odd . It should be noted that the formula (20) was previously obtained in [23], for the case n=0, =0. It was used for the phenomenological description of CR laser beams with M 2 >1. Also, Figs. 3(d)-(f) shows that if we move away from the focal plane by a distance ξ, CR cones and corresponding BG beams will start to spread from and to the axis, respectively. This will lead to their spatial separation and formation of the Poggendorff's rings. Moreover, the separation will be bigger, the larger the value of the mode indices n or | |. One of the CR cones subsequently focuses near the axis, forming a Raman spot, the properties of which we will consider in the next section.

Axial spike and the Raman spot shift
In the focal plane and near it, in the approximation ρ 0 >>1, the CR field C m (±) for different values of m do not differ from each other, as is obvious from formulas (9)- (11). Nevertheless, near the ξ-axis where the radiation is focused, different components will form different spatial intensity distributions, since in (3) and (13) the order of the Bessel function depends on the index m. It should be noted that the "focusing" here is employed in a more specific meaning of the light intensity growth in the near-axis region, since the intensity value directly on the axis for +m>0 will be strictly zero. Let's consider how the components of the BG beam U(ρ,ϕ,-ξ)*, which correspond to the CR cone C (−) , behave for ρ≈0 and ξ>>1. First, consider the behaviour of the argument of the modified Bessel function in (8), explicitly describing the dependence on the longitudinal coordinate ξ: Since the condition ξ≈ρ 0 /a 0 is fulfilled near the Raman spot, the real part of expression (21) is reset to zero and only the imaginary part remains, which is correspondingly equal to -iρρ 0 /ξ and does not depend on σ 0 .
Since the argument Φ is completely imaginary, it is necessary to replace the modified Bessel function I +m (-iρρ 0 /ξ) with the Bessel function J +m (ρρ 0 /ξ), which is consistent with the results of Berry's paper [1], where the asymptotic of the CR field near the Raman spots was considered.
Besides, as shown earlier, the tilt parameter a 0 increases with the growth of n and , which is why the Raman spot should focus earlier, shifting closer to the focal plane. To demonstrate this, we need to consider the argument of the Gaussian function of expression (8) and, using the approximation ρ≈0 and ξ>>1, we get: From (22), it becomes obvious that the maximum value of such a function is reached at the point ξ≈ρ 0 /a 0 . This shifts the Raman spot closer to the focal plane with the growth of the tilt parameter a 0 . This statement is illustrated in Fig. 4, on which the CR axial intensity distribution for different values of the index of the input Laguerre-Gaussian mode is calculated. Fig. 4. The axial distributions of the CR beam intensity are computed numerically from the CR integrals (3) (first column), and the BG beam Eqs. (8) (second column) for different values of the OAM =0 (first row), =1 (second row) and =2 (third row) and ρ 0 =15. One can see that with the growth of , the Raman spots move closer to the focal plane. Also, we calculated the Raman spot intensity profiles (third column) from the exact solution (blue dotted curve) and the BG model (red solid curve) for different values of the OAM. Furthermore, for profiles we neglect contribution of the Bessel function in (22), which zeroes the intensity strictly on the axis for higher-order modes.

Renormalization of the Bessel-Gaussian model parameters
The transition from formula (13) to (14), in which the Bessel function is replaced by an exponential function, in the general case of an arbitrarily large orbital moment (within the paraxial approximation) becomes unsatisfactory for the fixed value of the CR parameter ρ 0 . As a result, we predict that the radius of the CR ring is no longer exactly equal to R 0 (ρ 0 in dimensionless coordinates) and increases with the growth of , as shown in Figs. 5(a)-(c). We can give a simple explanation for this shift if one remembers that the maximum intensity of the optical vortex (1) is located on the ring of radius r max = w 0 | |(ρ max = | |). Therefore, when one moves from the parameter area where ρ 0 >> ρ max >>1 to ρ 0 ≈ 0, there should be a smooth change from the maximum of the CR ring at the point ρ=ρ 0 to ρ max , due to which the center of the ring begins to depend on as shown in Figs. 5(b) and (c). Thus, for the case of large , it is necessary to renormalize the parameter ρ 0 in our BG model. As will be shown below, we will need to perform this renormalization for the tilt parameter a 0 and for the normalization complex constant A 0 . We assume that the description in terms of generalized BG beams is valid and the Fourier-Bessel-Hankel transform for the renormalized BG model beam will have the form: where we introduce the notationμ =ρ 0 + iσ 0 2ã 0 with already renormalized parameters to be found. It is worth mentioning that, we will calculate only first nonvanishing correction to the parameters of the BG model, since the approximation ρ 0 > > 1 is still valid. Thus, the rough criterion for the applicability of our calculations is the condition < ρ 0 .
So to take into account the effects of large values of OAM , it is necessary to use the asymptotic expansion of the Bessel function for large complex arguments and retain component proportional to 1/κ [15]: Decomposing the function 1/κ near κ=a 0 up to a linear term and choosing the renormalized parameters so that (23) equals exactly to (14), after simple algebra, we get: where we introduced the notation ∆= 2 /(2a 0 2 ). From the formula (25), as a result, we can conclude that the center of the CR cone actually shifts from the point ρ=ρ 0 with the growth of , and the tilt parameter a 0 will increases with the growth of , not as fast as before renormalization. Thus, the dependence of the CR ring radius R 0 on the value of OAM can be naturally included in our proposed BG model.

Conclusion
In conclusion, we demonstrate that the CR of the input elegant Laguerre-Gaussian beams can be effectively described through generalized BG beams. This representation has an integrationfree structure and is valid in all spatial regions, including the multi-ring Lloyd's plane, the Poggendorff's transitional region, and near the Raman spots. The quantitative connection between CR of the elegant Laguerre-Gaussian beams and generalized BG beams is established in terms of the proposed determination procedure for BG beam parameters. We performed numerical simulations for different mode indexes of the input elegant Laguerre-Gaussian beam. It was demonstrated a good agreement between exact solution and our proposed Bessel-Gaussian approximation model for small indexes. For large mode indexes, our solution perfectly matches with the exact one. Furthermore, the application of the proposed BG model elegantly explains the transition of the classical double-ring pattern of CR in the Lloyd plane into a multi-ring one reported earlier for simplest optical vortexes with >0 [13]. It is noteworthy that mathematical elegance and physical clarity of the proposed BG model have allowed us to predict new phenomenon such as Raman spot shift and dependence of the CR ring radius on the value of OAM. Now we outline some implications of the above-mentioned results. The proposed BG model of CR, in the future, can be extremely beneficial for calculating stable modes of laser resonators with the CR crystal as the active medium [4,5,24]. Besides, it may boost new studies on the CR laser resonators, where axicon-shape mirrors are used [25]. Also, because elegant Laguerre-Gaussian beams are the complete set of solutions to the paraxial wave equation in free space, the proposed BG model may be very helpful for multi-mode CR applications such as for input incoherent or semiconductor laser sources [23,26].

Disclosures
The authors declare no conflicts of interest.