1. Introduction

Fractional calculus is the theory of differentiation and integration to arbitrary fractional or even complex order [1, 2]. Although the theory dates back to 1695 when Leibniz presented his analysis of the derivative of order one half, for three centuries it developed mainly as a pure theoretical field of mathematics. In the last few decades there has been a revival of interest in fractional calculus because its theory has had applications in several areas, such as fractals, viscoelasticity, transmission lines, and electrodynamics, to name a few [1, 2, 3, 4]. Optics is a field in which the use of conventional calculus plays a major role, and it is of interest to see the potential application of fractional calculus in this field. In this direction, the fractionalization of the Fourier transform and its optical applications have been already studied for the last 15 years [5, 6], and more recently, the fractionalization of the curl operator was applied also in electromagnetism [7] and metamaterials [8].

The present work is part of our efforts in exploring the applications of fractional calculus and fractional operators to optics and beam theory. In particular, in a previous paper [9], we introduced fractional-order beam solutions to the paraxial wave equation (PWE) that constitute continuous transition modes between the Hermite–Gaussian beams of integral-order [10]. These new solutions were obtained formally by applying the fractional derivative operator on the fundamental Gaussian beam.

In this paper, we extend our previous work to the case of cylindrical geometry, and present new fractional-order solutions that smoothly connect the Laguerre–Gaussian (LG) beams of integral-order [10]. The new circular fractional solutions are characterized in general by two fractional indices and are obtained by fractionalizing the creation operators used to create integer-order LG beams from the fundamental Gaussian beam. Remarkably, the mathematical form of the fractional solutions in circular coordinates turned out to be much simpler than the expression for the fractional beams in Cartesian coordinates [9]. This fact allowed us to determine the width, the moments, and the M ² factor of the circular fractional beam in closed form. We discuss also the adjoint beam solutions and establish the conditions for beam biorthogo-nality. The orbital angular momentum carried by the fractional beam is a continuous function of the angular mode index and it is not restricted to take only discrete values. In this paper we focus our attention to the fractionalization of the elegant LG beams, and leave the fractionalization of the standard LG beams for a future article. We finally remark that this paper is fully self-contained as it introduces all relevant definitions as well as main motivations.

2. Fractional operators and fractionalization of elegant Laguerre–Gaussian beams

We begin the analysis by writing the complex field amplitude of the normalized Gaussian beam with time dependence exp(-iωt) propagating in free-space along the positive z axis of a coordinate system r = (x,y, z) = (rcosθ, rsinθ, z):

where w ₀ is the beam width at the waist plane z = 0, and z_R = kw ² ₀/2 is the Rayleigh distance, and k is the wave number. The field in Eq. (1) carries unit power (i.e.∫∫^∞ _-∞|U _0,0|²dxdy = 1), and constitutes a fundamental solution of the PWE

where Δ₂ is the transverse Laplacian operator.

It is known [11,12] that normalized elegant LG beams with radial n= (0,1,2,…) and angular l = (0,1,2,…) mode numbers

can be derived from the fundamental Gaussian beam U _0,0 by the repeated application of the differential creation operators

applying the following operator prescription

where L_n^l is the associated Laguerre polynomial [13], and (!) denotes the factorial operation. The basis for such a construction is the following theorem: let U be a solution of the linear operator L (i.e. LU = 0), if other linear operator A commutes with L, then the function AU is also a solution of L. For the sake of simplicity, throughout the paper we restrict the analysis to positive values of the azimuthal index. Expressions for negative index can be straightforwardly derived from the positive ones by symmetry considerations.

To explore the application of fractional operators in beam theory, we look for new fractional-order solutions U_η,λ of the PWE resulting from the operation

where η and λ are two positive arbitrary numbers denoting the fractional radial and angular indices of the beam, and c_η,λ is a normalization constant to be determined at a later stage. The operators (A ^±)^α (for arbitrary α ≥ 0) correspond to the fractionalization of the operators A ^± and satisfy the following conditions: (a) for α = 1, we get the original operators A ^±; (b) for α = 0, we obtain the identity operator, and (c) for two numbers α and β we have (A ^±)^α (A ^±)^β = (A ^±)^β (A ^±)^α = (A ^±)^α+β . The fractional operators (A ^±)^α commutes with the operator of the PWE, therefore we expect that the action of (A ^±)^α on any solution of the PWE will give a new solution of the same equation.

Evaluation of the fractional operators in Eq. (6) can be performed analytically in the frequency space. The two-dimensional Fourier transform (FT) of a function f(r,θ) and its inverse FT are written in cylindrical coordinates as

where (k_x,k_y) = (k_t cos ϕ,k_tsinϕ) denotes the transverse position in the frequency space.

To get the FT of Eq. (6), we replace ∂/∂x by ik_x and ∂/∂_y by ik_y. Therefore differential operators A ^± are replaced according to

The FT of Eq. (6) is then given by

where

is the FT of the fundamental Gaussian beam U ₀₀ (r, θ) .

Equation (9) describes the angular spectrum of plane waves for the desired fractional beam. From a physical point of view, the circular fractional beam can be derived by modulating the plane-wave spectrum of the fundamental Gaussian beam Ũ₀₀ with a modulation function proportional to k_t ^2η+λexp (iλϕ).

The field U_η,λ (r) is now determined by inverse Fourier transforming Eq. (9). Inserting Ũ₀₀ into Eq. (7b) we obtain

where d_η,λ = c_η,λ(-1)^η+λ (iw ₀/√2)^2η+λ+1/2πi√π is a constant.

Equation (11) is the integral representation for the fractional beam U_η,λ (r) and effectively corresponds to the evaluation of the fractional operators given in Eq. (6). In the following sections we will study in detail the solutions of Eq. (11) for both integral and fractional angular index λ.

3. Fractional beams with fractional radial index η and integer angular index l

In this section we consider the solutions Eq. (11) when the angular index λ becomes integer. As we will see, this condition leads to fractional beams whose transverse intensity distribution is circularly symmetrical, or equivalently, whose transverse field is separable into radial and angular parts.

Let λ= l be an integer number, for this case the angular integral in Eq. (11) can be evaluated in closed form in terms of the lth-order Bessel functions J_l , we have [14]

Substitution into Eq. (11) yields

Physically speaking, Eq. (13) corresponds to the continuous expansion of the field U_η,l in terms of l-th order Bessel beams of the form J_l(k_tr)exp(ilθ).

The evaluation of the integral in Eq. (13) [using Eq. (6.631.1) of Ref. [14]] leads to the desired expression for the circular fractional beam, namely

 

Fig. 1. Radial behavior of U_ηl (r,θ,z = 0) for the first four angular orders l = {0,1,2,3} and a continuous radial index variation 0 ≤ η ≤ 12.

Download Full Size | PDF

where Φ(a, b;z) is the Kummer confluent hypergeometric function [often denoted also as ₁ F ₁(a, b;z)], U _0,0 (r) is the Gaussian beam [Eq. (1)], and the normalization constant

has been determined to ensure unit beam power, i.e. ∫∫^∞ _-∞|U _ηl|² dxdy = 1 at any z plane. Here Γ(x) denotes the Gamma function.

Equation (14) is an exact solution of the PWE and is the main result of this work. It constitutes the radial fractionalization of the fundamental Gaussian beam in circular coordinates and indeed corresponds to the analytical evaluation of the fractional operator relation

for arbitrary positive values of the radial index η, and integers values of the azimuthal index l. The presence of the Gaussian beam in Eq. (14) ensures the physical requirement that the field amplitude vanishes for r arbitrarily large, and that the beam is square integrable.

In Fig. 1 we show the radial behavior of U_η,l (r,θ,z = 0) at the waist plane for the first four angular orders l = {0,1,2,3} and a continuous radial index variation 0 ≤ η≤ 12. At this plane μ = 1, thus the radial part of U_η,l (r) becomes a purely real function. The function Φ creates maxima, minima, and beam nulls in the amplitude distribution as η increases. In particular, U_ηl (r, θ,z = 0) has n radial zeros when η falls in the interval (n - 1) < η ≤n.

While the radial part of the fractional beam U_η,l (r) at the waist plane z = 0 is purely real, outside this plane it becomes complex leading to a continuous variation of the transverse pattern. This effect is illustrated in Fig. 2, where we show the transverse amplitude and phase of the fractional beam U _{2.5, 2}(r) at z = {0,0.5z_R,z_R} . Because of its azimuthal dependence of the form exp(ilθ), the fractional beam U_η,l (r) is azimuthally symmetric, and carries an intrinsic orbital angular momentum of lh̄ per photon which is independent of the fractional radial index η.

From the theory of the confluent hypergeometric functions [see Eq. 8.972.1 of Ref. [14]], when η becomes a positive integer n = 0,1,2,…, then Φ(-n,l+ 1;x) = n!l!L^l_n(x) /(n + l)! and thus Eq. (14) reduces to the elegant LG beams given by Eq. (3) and shown with solid lines in Fig. 1. So, effectively, the integer-order solutions have been “smoothly connected” by varying the order of fractional differentiation of the Gaussian beam.

 

Fig. 2 Amplitude and phase of the beam U_η,l (r) atz= {0,0.5z_R,z_R} for η = 2.5 and l = 2. Square image size of 6w ₀ × 6w ₀ .

Download Full Size | PDF

3.1. Moments, width, and M² factor of the circular fractional beams

The significant parameters to characterize the circular fractional beam U_η,l (r) are determined in terms of the moments of the intensity distribution and its Fourier transform. The first-order moment provides the position of the centroid of the beam on the transverse plane, and vanishes by virtue of the circular symmetry of U_η,l (r).

The second-order moments σ₀ and σ_∞ associated with the intensity distributions at the waist and at the far field, respectively, provide the irradiance spot size w = 2σ₀ and the quality factor M ² = 2πσ₀σ_∞ of the beam. For circularly symmetric beams, the moments σ₀ and σ_∞ are given by

where U_η,l and Ũ_η,l are given by Eqs. (14) and (9) evaluated at z = 0, respectively. Taking into account that Eqs. (14) and (9) are already normalized [i.e. both denominators in Eqs. (17) are unity] we have for σ₀ and σ_∞

These integrals can be evaluated analytically (see Appendix A), and the results are given by the following surprisingly simple formulas

Thus the irradiance spot size and the M ₂ factor of the fractional beam U_η,l (r) turns out to be

As expected, the M ² factor increases with orders (η,l). For integer values of η, Eq. (22) reduces to the M ² factor of the elegant LG beams [15, 16, 17]. It is important to note that the parameters defined in this section must be taken as global parameters of the circular fractional beams. They do not describe local variations of the intensity distribution.

3.2. The fractional radial creation operator

In this section we discuss the differential equation satisfied by the radial part of the fractional beams U_η,l, and derive a useful fractional creation operator for this radial part which allows to rewrite the two-variable operator prescription [Eq. (16)] using operators depending exclusively on derivatives with respect to r.

We commence by noting from Eq. (14) that U_η,l (r) can take the separated form U_ηl (r) = f_ηl (R)exp(ilθ), with f_η,l(R) being the radial part given by

where R ≡ r/ √μw ₀ is the normalized radius, and unnecessary constants have been omitted. Taking advantage of the fact that Φ fulfills the confluent hypergeometric equation [13] it can be demonstrated (see Appendix B) that f_ηl (R) is solution of the ordinary differential equation

where Q denotes the linear operator of the equation.

From Eq. (16) and the commutability of the operators A ^± , it is clear that if U_η,l is a solution of the PWE with azimuthal index l, then A ^± U_η,l is a solution with index l ± 1, that is

The special case of the repeated application of A ⁺ on a function g(r) which depends exclusively on r leads to the operator relation

where m is an integer.

For arbitrarily positive β, the individual application of the fractional operators (A ⁺)^β or (A ^-)^β on U_η,l leads to azimuthally asymmetric solutions of the PWE (see Sect. 3 for details). By applying the operator (A ⁺)_β (A ^-)^β on U_η,l it is possible to cancel out the opposite angular effects, and consequently to modify the radial index η of U_η,l by keeping constant its azimuthal index l. From the definition of A ^± in Eq. (4), the operator A ⁺ A ^- is recognized to be the same as the transverse Laplacian operator Δ₂ acting on a field with azimuthal dependence exp(ilθ) (i.e. the angular derivative ∂/∂θ is then replaced by il). We have explicitly

We then conclude that the fractional operator Δβ ₂ is indeed a radial creation operator for the field U_η,l.

The effect of acting the fractional operator Δ^β ₂ on the radial function f_η,l (r) is determined also by noting that the commutator of Δ^β ₂ and the operator Q in Eq. (24) is given by [Δ^β ₂, Q] = 4β Δ^β ₂. It follows that the new function Δ^β ₂ f_η,l satisfies the same differential equation (24) with the parameter η changed to η+β. We then conclude that the effect of Δ^β ₂ on f_η,l is simply to change the order, while remaining the original functional form, i.e. Δ^β ₂ f_η,l = f_η+β,l.

We complete this section mentioning that results discussed above allow us to rewrite the two-variable operator prescription in Eq. (16) in terms of operators depending on radial derivatives exclusively. Omitting unnecessary constants we write

Now, by virtue of Eq. (27) with l = 0, the operation (A ⁺)^η (A ^-)^η U ₀₀ is fully equivalent to [∂_r ² + (1/r)∂_r]^η U ₀₀, which yields a function depending on r only. Since l is an integer number, then from Eq. (26) we finally have

for arbitrarily positives values of η. Equation (29) allows to generate a fractional beam U_η,l starting from the fundamental Gaussian beam using only radial operators. It is worth mentioning that this operator formula (with η being an integer number) has been applied recently in Ref. [18] to determine the higher-order complex source for the elegant LG non paraxial waves.

3.3. On the adjoint radial equation and adjoint fractional beams

The adjoint equation to Eq. (24) is found to be

where ρ = R ^* = r/w ₀√μ ^* and (^*) denotes complex conjugate. The solutions to the adjoint equation are the adjoint fractional beams given by

Comparing with respect to functions f_η,l(R), note that there is no Gaussian factor associated with the adjoint functions. It is now clear that the radial equation (24) is not a self-adjoint equation, then its solutions f_η,l do not form an orthonormal set. Applying the theory of the confluent hypergeometric functions [13] it is possible to show that the biorthogonality integral ∫^∞ ₀ f_η,l(r)f̂_γ,l(r) rdr between f_η,l and its adjoint functions f̂_γl diverges for arbitrary values of η and γ unless both η and γ become integer numbers (that indeed is the case of the known biorthogonal relation for the elegant LG beams [10]). We then conclude that it is not possible to formulate an orthogonality relation for functions f_η,l (r) with arbitrary η in the semi-infinite domain 0 ≤ r < ∞. Nevertheless, it may be established in a finite domain by converting Eq. (24) into a self-adjoint form and applying the Sturm-Liouville theory.

4. Discussion of the case when the angular index λ is not integer

Let us now discuss the general case when the azimuthal index λ of the fractional beam U_η,λ (r) is not an integer number. As we will see, this case leads the transverse field to be non-circularly symmetrical.

As discussed in Sect. 2, the fractional beam results from the operation U_η,λ∝ (A ⁺)^η+λ (A ^-)^η U ₀₀, and its integral representation is given by Eq. (11), namely

 

Fig. 3. Transverse amplitude and phase distributions of the fractional beam U_η,l (r) at z = 0 for several values of (η, λ). Last row shows the z-component of the OAM carried by the beam for constant λ. Square image size of 5w ₀ × 5w ₀.

Download Full Size | PDF

The point here is that it appears that the angular integral ∫^2π ₀ exp(iλϕ)exp[ik_trcos(ϕ - θ)]dϕ cannot be evaluated in closed form for non integer values of λ. Although there is always the possibility of using numerical methods to evaluate directly Eq. (32), an alternative expression for U_ηλ (r) may be derived by expanding exp (iλϕ) in its Fourier series as follows:

After inserting this expansion into Eq. (11), we can integrate each term of the summation following the same procedure described in Sect. 3 for integer values of l. The result is

 

Fig. 4. Behavior of (a) the position of the centroid x_c, and (b) the OAM J_z carried by the beam in function of the angular index λ.

Download Full Size | PDF

where f _η′,|l| (R) is given by Eq. (23), and an overall amplitude factor has been omitted for simplicity.

Equation (34) says that a fractional beam with arbitrarily fractional indices (η,λ) can be constructed with a suitable superposition of fractional beams with integer indices l. Figure 3 shows the transverse amplitude and phase distributions of the fractional beam U_η,λ at z = 0 for several values of (η, λ) in the ranges η∈ [5,6] and λ ∈ [1,2]. The patterns were obtained by adding 101 terms of the series in Eq. (34) from l = -50 to l = 50. As the radial and the angular indices increase, the irradiance and phase patterns vary continuously exhibiting an azimuthally asymmetric shape which becomes circularly symmetrical only when λ is integer.

For fractional values of λ the centroid of the beam U_η,λ is slightly displaced from the origin in the horizontal direction. Figure 4(a) shows the position of the beam centroid x_c in normalized units of w ₀ as a function of the angular index λ. The curve was determined by evaluating numerically the definition of the first-order moment x_c = ∫∫^∞ _-∞ x|U|²dxdy/∫∫^∞ _-∞|U|² dxdy for the range λ ∈ [0,6] and η = 0.. As λ increases x_c exhibits a decreasing oscillatory behavior around the origin and vanishes for integer λ.

It is important to see what happens with the orbital angular momentum (OAM) carried by the beam when both indices vary. Within the paraxial regime, the z component of the OAM per photon in unit length about the origin of a transverse slice of a beam U (r) is given by [19]

where r _t = x x̂ +y ŷ is the transverse radius vector. From the last paragraph of Sect. 3.2, we know that U_η,λ is proportional to

Now, it is clear that the function U _η,0 depends on r only and thus it does not carry OAM. We then conclude that the existence of OAM in U_η,λ comes from the application of (A ⁺)^λ on U _η,0, and thus it follows that J_z depends on λ only, while it is independent on η. To determine J_z, we evaluated numerically Eq. (35) using a two-dimensional Gauss-Legendre quadrature for a large number of combinations of radial and angular indices (η,λ). The numerical analysis corroborated that the OAM carried by the beam is independent of the radial index η, and a continuous function of the angular index λ. For the beam shapes shown in Fig. 3, the values of J_z are included in the bottom line. Figure 4(b) depicts also the behavior of J_z as a function of the angular index within the range λ ∈ [0,6]. It is clear that by adjusting the value of λ it is possible to tune OAM carried by the fractional beam.

5. Conclusions

We have introduced new fractional-order solutions of the PWE in cylindrical coordinates that constitute continuous transition modes between the integer-order elegant LG beams. The new solutions [Eqs. (14) and (34)] were obtained from the fundamental Gaussian beam by defining appropriate fractional creation operators which commute with the operator of the PWE. In the course of obtaining the new solutions, we also determined in closed form its Fourier transform and its expansion in terms of Bessel beams. The most important results of this work are summarized as follows:

Some of the interesting issues to pursue is whether the concept of fractional beams can be extended to multipole complex-source point solutions of the Helmholtz equation [18]. Of interest are also the effect of making the fractional parameters (η,λ) complex, and the behavior of the vortex structure exhibited by the fractional beam under a continuous variation of the fractional indices. These are currently under study by the author.