Circularly Symmetric Light Waves: An Overview

TL;DR

This paper provides a comprehensive analysis of vortex solutions to the Helmholtz equation, including scalar and vector beams, emphasizing their angular momentum properties and propagation characteristics.

Contribution

It introduces a rigorous operator-based framework for classifying and analyzing vortex solutions, extending understanding of OAM beams and their polarization-spatial coupling.

Findings

Characterization of scalar and vector vortex solutions

Analysis of propagation properties of OAM beams

Connection with two-dimensional rotation group

Abstract

Orbital Angular Momentum (OAM) waves were first recognized as those specific vortex solutions of the paraxial Helmholtz equation for which the orbital contribution to the total angular momentum of the beam yields an integer multiple of $\hbar$ along the propagation direction. However, this class of solutions can be generalized to include more sophisticated vector vortex waves with coupled polarization and spatial complexity, that are eigenfunctions of the third component of the angular momentum operator. In this work, a rigorous framework is proposed for the analysis of all the possible families of vortex solutions to the homogeneous Helmholtz equation. Both the scalar and vector cases are studied in depth, making use of an operator approach which emphasizes their intimate connection with the two-dimensional rotation group. Furthermore, a special focus is given to the characterization of the propagation properties of the most popular families of paraxial OAM beams.

Equations

$$\frac{1}{2\pi }{\int }_0^{2\pi }{\Phi }_j^{\ast }\left(\phi \right){\Phi }_m\left(\phi \right)d\phi ={\delta }_{jm};\sum _m{\Phi }_m\left(\phi \right){\Phi }_m^{\ast }\left({\phi }^{\prime }\right)=\delta \left(\phi -{\phi }^{\prime }\right).$$

$$J_{1}=\left[\begin{array}[]{ccc}0&0&0\\ 0&0&-1\\ 0&1&0\end{array}\right];\quad J_{2}=\left[\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ -1&0&0\end{array}\right];\quad J_{3}=\left[\begin{array}[]{ccc}0&-1&0\\ 1&0&0\\ 0&0&0\end{array}\right],$$

$${\nabla }^2\Psi \left(x_1,x_2,x_3,t\right)-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}\Psi \left(x_1,x_2,x_3,t\right)=0,$$

$${\nabla }^2\psi \left(x_1,x_2,x_3\right)+k^2\psi \left(x_1,x_2,x_3\right)=0.$$

$$i{\hat{J}}_z\psi =-i\frac{\partial }{\partial \phi }\psi =m\psi \Rightarrow \psi \left(x_1,\phi ,x_3\right)={\Phi }_m\left(\phi \right)u\left(x_1,x_3\right),$$

$$\frac{1}{h_1{h}_{\phi }{h}_3}\left[\frac{\partial }{\partial x_1}\left(\frac{h_{\phi }{h}_3}{h_1}\frac{\partial u}{\partial x_1}\right)+\frac{\partial }{\partial x_3}\left(\frac{h_1{h}_{\phi }}{h_3}\frac{\partial u}{\partial x_3}\right)\right]+\left(k^2-\frac{m^2}{h_{\phi }^2}\right)u=0,$$

$$\left[\frac{{\partial }^2}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial }{\partial \rho }-\frac{m^2}{{\rho }^2}+\frac{{\partial }^2}{\partial z^2}+k^2\right]u\left(\rho ,z\right)=0.$$

$${\hat{P}}_{z}u=\frac{\partial }{\partial z}u=-i{k}_{z}u\Rightarrow u\left(\rho ,z\right)=B\left(\rho \right)e^{-i{k}_{z}z},$$

$$\left[\frac{{\partial }^2}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial }{\partial \rho }-\frac{m^2}{{\rho }^2}+\left(k^2-k_z^2\right)\right]B\left(\rho \right)=0.$$

$${\Psi }_m^B\left(\rho ,\phi ,z,t;k_{\rho }\right)=C_m^B{J}_{|m|}\left(k_{\rho }\rho \right){\Phi }_m\left(\phi \right)\exp \left(-i{k}_{z}z+i\omega t\right),$$

$${\tilde{\Phi }}_m\left(\phi \right)=\sqrt{2}\left\{\cos ϵ\cos \left[m\left(\phi -{\phi }_0\right)\right]+i\sin ϵ\sin \left[m\left(\phi -{\phi }_0-\sigma \right)\right]\right\},$$

$$\left|\frac{{\partial }^{2}f\left(\rho ,\phi ,z\right)}{\partial z^2}\right|\ll k\left|\frac{\partial f\left(\rho ,\phi ,z\right)}{\partial z}\right|;\left|\frac{{\partial }^{2}f\left(\rho ,\phi ,z\right)}{\partial z^2}\right|\ll \left|{\nabla }_T^{2}f\left(\rho ,\phi ,z\right)\right|,$$

$$\left[\frac{{\partial }^2}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial }{\partial \rho }+\frac{1}{{\rho }^2}\frac{{\partial }^2}{\partial {\phi }^2}-2ik\frac{\partial }{\partial z}\right]f\left(\rho ,\phi ,z\right)=0.$$

$$\left[\frac{{\partial }^2}{\partial {\varrho }^2}+\frac{1}{\varrho }\frac{\partial }{\partial \varrho }+\frac{1}{{\varrho }^2}\frac{{\partial }^2}{\partial {\phi }^2}-4i\frac{\partial }{\partial \zeta }\right]f\left(\varrho ,\phi ,\zeta \right)=0.$$

$$f\left(\varrho ,\phi ,\zeta \right)=\frac{C}{\mu \left(\zeta \right)}u\left[\frac{\varrho }{\mu \left(\zeta \right)},\zeta \right]\exp \left[i\varphi \left(\zeta \right)\frac{{\varrho }^2}{{\mu }^2\left(\zeta \right)}+im\phi \right],$$

$$i\frac{\partial u\left(r,\zeta \right)}{\partial \zeta }=\hat{H}\left(\zeta \right)u\left(r,\zeta \right)$$

$$\hat{H}\left(\zeta \right)$$

$$\hat{P}=\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r};\hat{Q}=i\left(r\frac{\partial }{\partial r}+1\right).$$

$$a=\mu \frac{d\mu }{d\zeta }+\varphi ;b={\mu }^2\left[\frac{d\varphi }{d\zeta }-\frac{{\varphi }^2}{{\mu }^2}-\frac{2\varphi }{\mu }\frac{d\mu }{d\zeta }\right],$$

$$i\frac{\partial u\left(r,\zeta \right)}{\partial \zeta }=\left\{\frac{1}{4{\mu }^2\left(\zeta \right)}\left[\hat{P}+4a\hat{Q}-\frac{m^2}{r^2}\right]+\frac{b{r}^2}{{\mu }^2\left(\zeta \right)}\right\}u\left(r,\zeta \right).$$

$$i\frac{\partial u_B\left(r,\zeta \right)}{\partial \zeta }=\frac{1}{4}\left[\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{m^2}{r^2}\right]u_B\left(r,\zeta \right),$$

$$ik\frac{\partial u_B\left(\rho ,z\right)}{\partial z}=\frac{1}{2}\left[\frac{{\partial }^2}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial }{\partial \rho }-\frac{m^2}{{\rho }^2}\right]u_B\left(\rho ,z\right).$$

$$k_z=\sqrt{k^2-k_x^2-k_y^2}\approx k-\frac{k_{\rho }^2}{2k},$$

$$u_B\left(\rho ,z\right)\propto B\left(\rho \right)\exp \left(i\frac{k_{\rho }^{2}z}{2k}\right),$$

$${\rho }^2\frac{d^{2}B\left(\rho \right)}{d{\rho }^2}+\rho \frac{dB\left(\rho \right)}{d\rho }+\left({\rho }^2{k}_{\rho }^2-m^2\right)B\left(\rho \right)=0.$$

$$f_m^B\left(\rho ,\phi ,z;k_{\rho }\right)=C_m^B{J}_{|m|}\left(k_{\rho }\rho \right)\exp \left(i\frac{k_{\rho }^{2}z}{2k}+im\phi \right).$$

$$i\frac{\partial u_{LG}\left(r,\zeta \right)}{\partial \zeta }=\frac{1}{4{\mu }^2\left(\zeta \right)}\left[\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{m^2}{r^2}-4r^2\right]u_{LG}\left(r,\zeta \right).$$

$$u_{LG}\left(r,\zeta \right)\propto r^{|m|}LG\left(2r^2\right)\exp \left[-r^2-i\Theta \left(\zeta \right)\right],$$

$$4r^{|m|+2}\left(L{G}^{{}^{\prime \prime }}-L{G}^{{}^{\prime }}\right)+r^{|m|}\left\{2\left(\left|m\right|+1\right)L{G}^{{}^{\prime }}-\left[\left|m\right|+1+{\mu }^2\frac{d\Theta }{d\zeta }\right]LG\right\}=0,$$

$$\left|m\right|+1+{\mu }^2\left(\zeta \right)\frac{d\Theta \left(\zeta \right)}{d\zeta }=-2p.$$

Full text

1 Introduction

Electromagnetic waves carry energy and both linear and angular momenta. Whereas a contribution to the total angular momentum is realized by the spin of the photon, the fundamental physical quantity associated with polarization, an orbital angular momentum (OAM) component can also be present. Such contribution is often said to be “quasi-intrinsic” [1], since it is independent of the axis of calculation for any beam endowed with a helical wavefront apertured symmetrically about the beam axis [2]. The OAM content of helical beams in the paraxial regime was first discovered by Les Allen et al. in 1992 [3], although these waves have been known for much longer in terms of optical vortices and are now often referred to as vortex waves. OAM beams are characterized by an azimuthal phase dependence of the form $\exp\left(im\varphi\right)$, where the integer $m$ represents the topological charge of the vortex (or, more simply, the number of intertwined helices the wavefront is made of) and is related to the OAM carried by the beam along the propagation axis. For $m\neq 0$, due to the presence of the on-axis phase singularity, scalar OAM waves present a doughnut-shaped intensity profile with a central null.

It would seem natural to think that, being the OAM a global property of the beam in contrast to the spin angular momentum (SAM), the orbital and spin contributions to the total angular momentum of an electromagnetic wave are two distinct quantities. However, it has been argued that, in view of the absence of any rest frame for the photon, a gauge-invariant decomposition of the total electromagnetic angular momentum in orbital and spin parts is unfeasible (see for example [4, 5]) or achievable at most within the paraxial limit [6], the polarization and spatial degrees of freedom being in general coupled. Despite these arguments, Steven J. van Enk and Gerald Nienhuis have claimed that it is indeed possible to extract gauge-independent expressions for the spin and orbital parts of the angular momentum, but neither of these contributions alone represents a true angular momentum, since the transversality of the radiation fields affects the commutation relations for the associated quantum operators [7]. Such interpretation proves to be the currently accepted one [8].

The first experimental proof of an angular momentum transfer between light and matter dates back to 1936, when Richard A. Beth verified the interaction between a circularly polarized beam and a quarter-wave plate [9]: owing to the photon spin, the torque exerted by the polarized radiation on the plate put it into rotation. Since then, a great number of similar experiments have been performed, also involving the use of OAM beams. The angular momentum of light can be transferred to suitable trapped material particles causing them to rotate [10, 11], with relevant applications in both micromanipulation and the design and operation of micromachines [12, 13, 14, 15]. Vortex waves are also used as a tool for advanced cell manipulation [16] and to transport aqueous droplets in solution [17]. Moreover, direct transfer of OAM from light to atoms in ultracold gas clouds has been reported [18, 19, 20].

The interest in OAM waves extends far beyond the light-matter angular momentum exchange and several attractive applications can be found in optics, quantum physics and astronomy. For instance, light vortices are employed to probe different physical and biological properties of matter in optical imaging [21, 22, 23], but also for encoding quantum information via the corresponding photon states in high-dimensional Hilbert spaces [24, 25, 26, 27, 28]. OAM beams are involved in second harmonic generation processes [29, 30, 31] as well as in many other nonlinear optics phenomena [32, 33, 34, 35], being intended in terms of vortex solitons [36, 37, 38]. Further possibilities lie in the use of vortex waves to enhance the resolution in quantum imaging [39] and for overcoming the Rayleigh limit of the telescopes [40], whereas several other interesting studies are reported by the astrophysical community [41, 42, 43, 44, 45]. Many of the so far described applications together with other examples and a more complete list of references can be found in [46, 47, 48, 49, 50, 51].

Currently, a large part of the research on OAM beams is also devoted to the field of telecommunications. From a strictly mathematical point of view, vortex waves represent classes of solutions to the Helmholtz equation; since each class comprises a set of orthogonal solutions, in 2004 Graham Gibson et al. proposed the idea of using such waves to convey independent information channels on a single frequency in free-space optical communication links [52]. The authors also suggested how the intrinsic sensitivity of the OAM orthogonality to angular restrictions and lateral offsets [53] could prevent eavesdropping.

In these last few years, a strong interest in the possibility to increase the communication efficiency by vortex waves is growing [54, 55, 56, 57, 58, 59, 60] and has also been extended to the radio frequency domain [61, 62, 63, 64, 65]. In this scenario, however, the degree of innovation of OAM multiplexing with respect to other existing techniques like those based on antenna diversity and, more generally, on the use of orthogonal wavefields has been questioned [66, 67, 68, 69, 70, 71, 72].

Motivated by the astonishing amount of OAM-related publications, a theoretical study is proposed in this paper on the mathematical features of free-space vortex waves. While presenting new results and general insights that can be of interest to the optics community, a large part of the literature on this fascinating topic is carefully reviewed. In Section 2, a detailed overview on the scalar OAM waves is outlined, starting from some basic results of group theory. In Section 3, the fundamental solutions of the paraxial Helmholtz equation are derived step by step through an ansatz approach and their propagation properties are systematically analyzed. In Section 4, the vector case is studied in depth from a new perspective and some well-known examples are revisited accordingly. The last section is devoted to the conclusions. Further calculations are reported in appendix as a support to the main text.

2 Vortex solutions of the scalar wave equation

As it is known in group theory, separable coordinate systems for second-order linear partial differential equations can be characterized in terms of sets of operators in the algebra relative to the underlying continuous symmetry group [73]. In this framework, the separated solutions of the considered equations are common eigenfunctions of the symmetry operators and the expansion of one set of separable solutions in terms of another leads back to a problem in the representation theory of the Lie algebra. The Helmholtz wave equation, one of the most studied examples, is proven to be separable in eleven three-dimensional coordinate systems [74]: Cartesian, circular cylindrical, elliptic cylindrical, parabolic cylindrical, spherical, prolate spheroidal, oblate spheroidal, parabolic, ellipsoidal, paraboloidal and conical. In optics, solutions of the Helmholtz equation which are eigenfunctions of the Lie algebra generator of the translations along the $z$-coordinate are usually considered, since only for them the equation can be separated into transverse and longitudinal parts; actually, this condition is satisfied in the first four orthogonal coordinate systems.

A useful approach for building a formal description of OAM waves is to analyze their definition from a very general point of view, proceeding along the lines of the methods and techniques developed in group theory. The fundamental feature which has been referred to in defining vortex waves from the very beginning is indeed represented by the aforementioned screw phase dislocation $\exp\left(im\varphi\right)$, whose presence is soon recognized as the footprint of the underlying circular symmetry. In particular, given a fixed axis along the unit vector $\mathbf{u}_{z}$, rotations about this axis form a one-parameter subgroup of SO(3) which is isomorphic to the group of rotations in the plane perpendicular to $\mathbf{u}_{z}$, that is SO(2) [75]. Through the azimuthal integer index $m$, the function $\varPhi_{m}\left(\varphi\right)=\exp\left(im\varphi\right)$ is known to span the set of single-valued irreducible representations of SO(2), for which the following orthogonality and completeness relations hold:

[TABLE]

Associated with the subgroup algebra there is a generator, $\hat{J}_{z}$, and all elements of the given subgroup can be written symbolically as $\exp\left(\varphi\hat{J}_{z}\right)$, where $\varphi$ parameterizes the rotation.

A basis for the Lie algebra of SO(3) in standard matrix representation is provided by the three operators:

[TABLE]

which obey the commutation relations $\left[J_{s},\,J_{v}\right]={\varepsilon_{sv}}^{k}J_{k}$, where $\varepsilon_{svk}$ is the Levi-Civita pseudo-tensor and the Einstein summation convention has been considered; in (2), the subscripts $1$, $2$ and $3$ denote Cartesian components $x$, $y$ and $z$, respectively. Although the origin of these commutation relations is really geometric in nature, the generators acquire even more significance in quantum physics, where they correspond to measurable quantities [76]. In this respect, $\mathbf{J}=\left(J_{1},J_{2},J_{3}\right)$ can be seen as the vector angular momentum operator in units of $\hbar$.

Whereas the matrix representation is widely used to express the SAM operator, a differential representation is often preferable to (2) for describing the OAM operator. As suggested in [77], analogies between quantum mechanics and scalar paraxial optics can be built in order to demonstrate heuristically that some cylindrical laser modes with a $\varPhi_{m}\left(\varphi\right)$ azimuthal phase dependence carry an OAM of $m\hbar$ per photon along the propagation axis.

Let us consider an arbitrary coordinate system $\left\{x_{1},x_{2},x_{3}\right\}$ and a free-space monochromatic scalar wave of the form $\Psi\left(x_{1},x_{2},x_{3},t\right)=\psi\left(x_{1},x_{2},x_{3}\right)\mathrm{e}^{i\omega t}$, where $\omega=kc$ represents the angular frequency, $k=2\pi/\lambda$ the wavenumber, $\lambda$ the wavelength and $c$ the speed of light. Even if the full vector nature of the electromagnetic radiation is being neglected at present, $\Psi\left(x_{1},x_{2},x_{3},t\right)$ can be understood as the amplitude of a linearly polarized electric field or vector potential, which satisfies the wave equation:

[TABLE]

as follows from Maxwell theory. From (3), the scalar Helmholtz equation for the spatial coordinates dependent amplitude $\psi\left(x_{1},x_{2},x_{3}\right)$ is soon derived:

[TABLE]

We are now interested in solutions $\psi$ to (4) that are eigenfunctions of the $\hat{J}_{z}$ operator, namely:

[TABLE]

where $\left\{x_{1},\varphi,x_{3}\right\}$ is a suitable coordinate system and the differential expression of the $\hat{J}_{z}$ operator has been introduced. It is important to note that the condition of periodicity of the wavefunction in the azimuthal coordinate imposes integer values for the $m$ index. In (5), $\varPhi_{m}\left(\varphi\right)$ corresponds to the sought vortex term and $u\left(x_{1},x_{3}\right)$ is a function obeying the reduced equation in the remaining two variables:

[TABLE]

being $h_{1}$, $h_{\varphi}$ and $h_{3}$ three $\varphi$-independent metric scale factors for the given coordinates [74].

It can be shown that (6) separates in five coordinate systems: circular cylindrical, spherical, prolate spheroidal, oblate spheroidal and parabolic [73]. While leaving to the appendix a formal proof of the equation separation in the latter four systems, let us focus in detail on the circular cylindrical case $\left\{\rho,\varphi,z\right\}$, for which equation (6) reads:

[TABLE]

As pointed out above, it is common to look for solutions of the Helmholtz equation that are eigenfunctions of the $\hat{P}_{z}$ operator, i.e. the generator of the translations along the $z$-axis:

[TABLE]

where $B\left(\rho\right)$ represents a function to be determined. The imposed condition implies a further reduction of equation (7):

[TABLE]

After some straightforward manipulations, (9) is easily recognized as the Bessel equation (see, for instance, [78]) and thus $B\left(\rho\right)=J_{\left|m\right|}\left(k_{\rho}\rho\right)$, where $k^{2}=k_{\rho}^{2}+k_{z}^{2}$.

Scalar modes of the form:

[TABLE]

where $C_{m}^{B}$ is a suitable dimensional constant, are known as Bessel beams (BBs) and represent a complete set of orthogonal solutions to (3). The transverse intensity and phase profiles of a representative BB are displayed in Figure 1 for different values of the topological charge $m$.

As will be shown further on, BBs also satisfy the paraxial wave equation. Due to their divergence-free profile, the Bessel modes belong to the family of “non-diffracting” beams [79], which are not square-integrable and thus carry infinite energy; for this reason, only truncated forms of such waves can be realized experimentally. Physical OAM waves that approximate the intensity distribution of these ideal beams are usually produced via computer generated holograms [80] or axicons [81].

It is important to emphasize that, despite the leading role played by the circular cylindrical coordinate system in optics, BBs are not the only possible scalar OAM solutions and waves with a vortex term $\varPhi_{m}\left(\varphi\right)$ are also found by separating the Helmholtz equation in spherical, prolate spheroidal, oblate spheroidal and parabolic coordinates (see A for further details). In [82], an estimation of the electromagnetic OAM has also been performed for Mathieu beams, which can be intended as a generalization of BBs to elliptic cylindrical coordinates; in this case, however, the $\varPhi_{m}\left(\varphi\right)$ term cannot be present due to symmetry considerations and this results in the appearance of a more complicated vortex structure with fractional OAM mean content.

The presence of a non-integer OAM per photon in connection with some sort of symmetry breaking is a fact, by the way, not new, that leads to possible extensions of the simple original definition of OAM waves in several interesting scenarios [83, 84, 85, 86, 87, 88, 89, 90]. As a related topic, it should be noted that the very concept of wave vorticity is not limited to the isotropic case, which has been taken as a working hypothesis in defining OAM beams as vortex solutions. In particular, a vortex is said to be isotropic when the phase increases linearly from [math] to $2\pi m$ around a circle enclosing the singularity [91, 92]; in this case, the iso-intensity contour lines are exact circles and the usual $\varPhi_{m}\left(\varphi\right)$ term is found. However, vortices are in general anisotropic, meaning that the intensity contours close to the singularities are elliptical [93, 94, 95, 96, 97, 98]. The anisotropic vortex can be understood as a geometrical deformation of the isotropic one, with an angular dependence of the form:

[TABLE]

where $\epsilon$ represents the vortex anisotropy, $\sigma$ is the vortex skewness and $\varphi_{0}$ the vortex rotation angle. The possibility of representing anisotropic dislocations through a superposition of isotropic ones, as reported in [98], provides an easy derivation of the expectation value and the uncertainty of the OAM.

A further aspect to be mentioned lies in the non-existence of propagating solutions characterized by a term $\varPhi_{\mu}\left(\varphi\right)$ with non-integer $\mu$ index: when programmed to have a non-integer phase dislocation of this form as a starting condition, electromagnetic beams evolve with the creation of a group of standard vortices which break the circular symmetry of the intensity profile. Such mechanism and all the related topics have been deeply explored and characterized [83, 84, 86, 87, 92].

The study of the electromagnetic vortices and their dynamics and interactions belong to the rich field of the singular optics, to which reference should be made. In spite of all the possible generalizations that can be introduced in the definition of OAM waves as vortex solutions, only the original requirement, i.e. the presence of an on-axis isotropic vortex, will be taken into account, since a complete classification would fall beyond the scope of the paper.

3 Paraxial OAM waves

In the event that the angle between the wavevector $\mathbf{k}$ and the propagation direction ($z$-axis) is small and assuming a beam amplitude of the form $\psi\left(\rho,\varphi,z\right)=f\left(\rho,\varphi,z\right)\exp\left(-ikz\right)$, where $f\left(\rho,\varphi,z\right)$ represents a slowly varying function of the $z$-coordinate, such that:

[TABLE]

being $\nabla_{T}^{2}$ the transverse Laplacian, equation (4) reduces to its paraxial version [99]:

[TABLE]

The paraxial Helmholtz equation (13) is widely used in optics, where the above described requirements are usually met and it represents a complete enough approximation to handle transverse variations and diffraction effects of the optical beam profile [100]. The reported equation is of the Schrödinger type and proves to admit solutions with separable variables in seventeen coordinate systems [73].

A huge amount of publications regarding the analytical derivation of families of paraxial beams can be found in the literature [101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129]. The most general set of solutions of (13) in circular cylindrical coordinates has been obtained and characterized recently [105, 107]; closed form expressions are provided in terms of the confluent hypergeometric functions [78]. Among all special cases of the paraxial vortex waves, some cylindrical families are of particular relevance and will be analyzed by means of an intuitive ansatz approach [130].

First, in order to simplify the notation, it is useful to introduce a set of dimensionless circular cylindrical coordinates $\left\{\varrho,\varphi,\zeta\right\}=\left\{\rho/w_{0},\varphi,z/z_{R}\right\}$, where $z_{R}=\pi w_{0}^{2}/\lambda$ and $w_{0}$ corresponds to a characteristic length parameter for the beam. Equation (13) then becomes:

[TABLE]

For cylindrical beams satisfying the requirements of the paraxial approximation it is useful to appeal to the following ansatz, which is eigenfunction of the $\hat{J}_{z}$ operator by definition:

[TABLE]

where $\mu\left(\zeta\right)$ and $\phi\left(\zeta\right)$ are two dimensionless functions which account for the diffraction effects, $m\in\mathbb{Z}$ represents the topological charge of the on-axis vortex and $C$ is a constant. By inserting expression (15) into (14), a partial differential equation for the function $u\left(r,\zeta\right)$ is found which can be put in the form of a $\zeta$-dependent Schrödinger equation:

[TABLE]

with Hamiltonian given by:

[TABLE]

where the following two operators have been introduced:

[TABLE]

As will be shown further on, a simple characterization of the cylindrical paraxial OAM beams can be derived from the definitions:

[TABLE]

which allow to rewrite equation (16) as:

[TABLE]

It is interesting to note that both $\hat{P}$ and $\hat{Q}$ are Hermitian operators with respect to the integration measure $rdr$, therefore the Hamiltonian $\hat{H}\left(\zeta\right)$ (term in braces in the above expression) is Hermitian if $a,\>b\in\mathbb{R}$ and $\mu^{2}\left(\zeta\right)$ corresponds to a real valued function of the $\zeta$-coordinate: in this particular case, complete sets of orthogonal solutions are obtained for the stationary Schrödinger equation. These and other classes of solutions are derived from the choice of the parameters $a$, $b$ and by establishing suitable functions $\mu\left(\zeta\right)$, $\phi\left(\zeta\right)$ which satisfy (19).

3.1 Bessel beams

With the choice $a=b=0$, $\mu\left(\zeta\right)=1$ and $\phi\left(\zeta\right)=0$, equation (20) reads:

[TABLE]

where the subscript $B$ has been introduced; alternatively, going back to the original coordinate system:

[TABLE]

For a paraxial beam of the form $\psi_{B}\left(\rho,\varphi,z\right)=g_{B}\left(\rho,\varphi,z\right)\exp\left(-ik_{z}z\right)$, the following relation holds:

[TABLE]

where $k_{\rho}^{2}=k_{x}^{2}+k_{y}^{2}$. Then we can write $\psi_{B}\left(\rho,\varphi,z\right)=f_{B}\left(\rho,\varphi,z\right)\exp\left(-ikz\right)$ with $f_{B}\left(\rho,\varphi,z\right)=g_{B}\left(\rho,\varphi,z\right)\exp\left(ik_{\rho}^{2}z/2k\right)$. Taking into account expression (15), we infer:

[TABLE]

which, substituted in equation (22), gives:

[TABLE]

Equation (25) corresponds exactly to (9) and, as mentioned above, its regular solution $B\left(\rho\right)$ is provided by the Bessel function of the first kind $J_{\left|m\right|}\left(k_{\rho}\rho\right)$, therefore leading to:

[TABLE]

The paraxial modes $\Psi_{m}^{B}\left(\rho,\varphi,z,t;k_{\rho}\right)=f_{m}^{B}\left(\rho,\varphi,z;k_{\rho}\right)\mathrm{e}^{-ikz+i\omega t}$ are nothing but the above discussed BBs (10). By virtue of their completeness and of the fact that they satisfy both the exact and the paraxial wave equations, BBs are well suited for implementing non-paraxial extensions of known beams families which belong originally to the set of paraxial solutions [48].

3.2 Laguerre-Gaussian beams

Let now $a=0$, $b=-1$, $\mu\left(\zeta\right)=\sqrt{1+\zeta^{2}}$ and $\phi\left(\zeta\right)=-\zeta$:

[TABLE]

In order to solve equation (27), it is convenient to start from the following assumption:

[TABLE]

where $LG\left(2r^{2}\right)$ and $\Theta\left(\zeta\right)$ are functions to be determined. After some straightforward algebraic steps, equation (27) reduces to:

[TABLE]

with the definitions $LG^{{}^{\prime}}=dLG\left(x\right)/dx$ and $LG^{{}^{\prime\prime}}=d^{2}LG\left(x\right)/dx^{2}$. Since an equation which holds for every $\zeta$ is needed, the coefficient with $\zeta$-dependence must be set equal to a constant, a conventional choice being $-\left(2p+\left|m\right|+1\right)$, with $p\in\mathbb{N}$:

[TABLE]

The equation is easily solved and gives:

[TABLE]

Hence, expression (29) becomes:

[TABLE]

whose regular solution is the generalized Laguerre polynomial $L_{p}^{\left|m\right|}\left(2r^{2}\right)$ [78] and leads to:

[TABLE]

In expression (33), the definition $w\left(z\right)=w_{0}\,\mu\left(z/z_{R}\right)$ has been introduced, $w_{0}$ corresponds to the waist of the Gaussian beam, $z_{R}$ is the Rayleigh distance (see Figure 2) and $\Theta\left(\zeta\right)$ the Gouy phase [100].

The constant term $C_{pm}^{LG}$ ensures correct dimension and normalization of the beam profile.

The cylindrical modes $\Psi_{pm}^{LG}\left(\rho,\varphi,z,t\right)=f_{pm}^{LG}\left(\rho,\varphi,z\right)\mathrm{e}^{-ikz+i\omega t}$ are named Laguerre-Gaussian (LG) beams and represent another complete set of orthogonal paraxial OAM solutions of the scalar wave equation. The transverse intensity and phase profiles of a representative LG beam are displayed in Figure 3 for different values of the topological charge $m$. Unlike BBs, the LG modes diverge during propagation (in this regard, $z_{R}$ can be seen as the diffraction scale), they are spatially confined and carry a finite amount of energy, therefore can be physically realized, at least up to some degree of accuracy.

LG beams are shape-invariant modes: during free propagation, their intensity profile maintains the same shape and it is just scaled owing to diffraction and energy conservation. The radius of the primary intensity maximum evolves with $z$ according to the law:

[TABLE]

where $\rho_{max,0}^{LG}=\rho_{max}^{LG}\left(0\right)$ vanishes for $m=0$. At large distances, for $m\neq 0$, we infer $\rho_{max}^{LG}\left(z\right)\propto z$. By replacing $\rho$ with expression (34) in the square modulus of (33), we easily get:

[TABLE]

being $I_{max,0}^{LG}=I_{max}^{LG}\left(0\right)$ the maximum intensity of the original profile. Therefore, for $z\gg z_{R}$, the asymptotic behavior $I_{max}^{LG}\left(z\right)\propto z^{-2}$ is obtained, as shown in Figure 4.

Let us now consider the decay of the intensity of a LG beam at a fixed radial coordinate $\bar{\rho}$ (for $m=0$, this also includes the case of the primary maximum at $\bar{\rho}=0$); the square modulus of (33) is given by:

[TABLE]

For large enough distances, the term which contains the Laguerre polynomial tends to a constant and so does the exponential one, thus implying an asymptotic behavior of the form $I_{LG}\left(z\right)\propto z^{-2\left|m\right|-2}.$ The same law also describes the far-field evolution of the LG beam power collected on a central surface of limited size placed at distance $z$, reported in Figure 5.

3.3 Hypergeometric beams

By choosing $a=-1/2$, $b=0$, $\mu\left(\zeta\right)=\sqrt{-\zeta}$ and $\phi\left(\zeta\right)=0$, equation (20) reduces to:

[TABLE]

If a field amplitude of the form:

[TABLE]

is assumed, expression (37) becomes:

[TABLE]

which holds true for any $\zeta$ under the following hypothesis:

[TABLE]

where $\gamma\in\mathbb{R}$. From (40), we get:

[TABLE]

and finally, from (39) and (40):

[TABLE]

This last equation admits a regular solution, given in terms of the Kummer hypergeometric confluent function ${}_{1}F_{1}\left[\left(\left|m\right|+1-i\gamma\right)/2,\left|m\right|+1;\,ir^{2}\right]$ (see, for example, [78]), which allows to obtain:

[TABLE]

where the dimensional constant $C_{m}^{H}\left(\gamma\right)$ is a function of the beam parameters and can be fixed by imposing the orthogonality condition.

The hypergeometric (HyG) beams $\Psi_{m}^{H}\left(\rho,\varphi,z,t;\gamma\right)=f_{m}^{H}\left(\rho,\varphi,z;\gamma\right)\mathrm{e}^{-ikz+i\omega t}$ constitute the third complete family of orthogonal OAM solutions of the scalar wave equation in the paraxial regime; unfortunately, like BBs, they carry infinite energy and are not square-integrable [120]. The transverse intensity and phase profiles of a representative HyG beam are displayed in Figure 6 for different values of the topological charge $m$. Equation (43) differs from the standard formula found in literature by a minus sign in the argument of the Kummer function, which is due to the choice of a propagation term of the form $\exp\left(-ikz+i\omega t\right)$ instead of $\exp\left(ikz-i\omega t\right)$, showing the two cases an opposite sign in front of the $z$-derivative in the corresponding versions of the paraxial Helmholtz equation.

3.4 Hypergeometric-Gaussian beams

Taking into consideration the possibility of using complex parameters, we can set $a=i/2$, $b=0$, $\mu\left(\zeta\right)=\sqrt{\zeta\left(\zeta+i\right)}$ and $\phi\left(\zeta\right)=-\zeta$:

[TABLE]

Let us suppose the function $u_{HG}$ to be factorized as:

[TABLE]

then (44) reduces to:

[TABLE]

Following the usual procedure, it is required that:

[TABLE]

where $\eta\in\mathbb{R}$. A quick integration gives:

[TABLE]

whereas the Kummer equation [78] arises from (46):

[TABLE]

[TABLE]

Expression (15) leads to:

[TABLE]

The modes represented by $\Psi_{m}^{HG}\left(\rho,\varphi,z,t;\eta\right)=f_{m}^{HG}\left(\rho,\varphi,z;\eta\right)\mathrm{e}^{-ikz+i\omega t}$ constitute the family of the hypergeometric-Gaussian (HyGG)* beams*, an overcomplete set of non-orthogonal paraxial OAM solutions of the scalar wave equation [116]. It can be proven that these configurations carry a finite power (and therefore can be experimentally approximated) as long as the condition $\eta\geq-\left|m\right|$ is satisfied.

The HyGG beams are not shape-invariant modes under free propagation; as far as the asymptotic behavior is concerned, the radius of the primary maximum of their intensity profile follows the same law as the LG case for $m\neq 0$, namely $\rho_{max}^{HG}\left(z\right)\propto z$, which implies the usual decay of the intensity maximum at large distances, $I_{max}^{HG}\left(z\right)=\left|f_{m}^{HG}\left[\rho_{max}^{HG}\left(z\right),\varphi,z;\eta\right]\right|^{2}\propto z^{-2}$. Also the asymptotic evolution of the intensity at fixed $\rho=\bar{\rho}$ is found to be $I_{HG}\left(z\right)\propto z^{-2\left|m\right|-2}$, as for the LG beams.

3.5 Hypergeometric-Gaussian beams of the second kind

Another interesting family of cylindrical modes is found by choosing $a=-i/2$, $b=0$, $\mu\left(\zeta\right)=\sqrt{1-i\zeta}$ and $\phi\left(\zeta\right)=0$. In this case, equation (20) reads:

[TABLE]

Let us suppose:

[TABLE]

which implies:

[TABLE]

Let now:

[TABLE]

whose solution is given by:

[TABLE]

Equation (54) has the same form as (46), thus $HG_{II}\left(r^{2}\right)=HG\left(r^{2}\right)$, expression (50). Finally:

[TABLE]

Paraxial modes of the kind $\varPsi_{m}^{HG_{II}}\left(\rho,\varphi,z,t;\eta\right)=f_{m}^{HG_{II}}\left(\rho,\varphi,z;\eta\right)\mathrm{e}^{-ikz+i\omega t}$ are known as hypergeometric-Gaussian type-II (HyGG-II) beams and, like the HyGG modes, represent an overcomplete set of non-orthogonal OAM solutions of the Helmholtz equation which are square-integrable for $\eta\geq-\left|m\right|$, condition under which they carry a finite amount of energy [115].

The HyGG-II beams are not shape-invariant modes under free propagation (see Figure 7). Following the standard convention, where a propagation term of the form $\exp\left(ikz-i\omega t\right)$ is considered, inducing a change of sign in the $z$-derivative in (13), the function $\mu\left(\zeta\right)$ has the form $\sqrt{1+i\zeta}$ and equation (57) changes accordingly.

Despite their name, resembling the previous family, the HyGG-II beams show significant differences with respect to the other known solutions. First, the radius of their primary intensity maximum increases with $z$ according to two possible asymptotic behaviors, depending on the choice of the $\eta$ and $m$ indices: when $\eta=-\left|m\right|\neq 0$, a square root law is found, $\rho_{max}^{HG_{II}}\left(z\right)\propto\sqrt{z}$, whereas the usual linear behavior $\rho_{max}^{HG_{II}}\left(z\right)\propto z$ is quickly reached by all the configurations in which $\eta>0$ for every $m$ or $\eta\geq 0$ for $m\neq 0$ (the doughnut-shaped asymptotic profile is always present, except for the standard Gaussian wave $\eta=m=0$). It is fundamental to understand that the square root divergence of the primary intensity maximum relative to the HyGG-II modes in the first subclass does not apply to the beam profile in its entirety, as can be inferred from the comparison reported in Figure 7.

Since the $\eta$ parameter must be greater than or equal to $-\left|m\right|$ in order to ensure square-integrability, the whole range $-\left|m\right|<\eta<0$ where $m$ is not zero can be explored: in this region, as a consequence of the very rippled and variable field pattern, the primary intensity maximum oscillates at different radii in the beam profile and it evolves alternating between the two possible above described laws as a function of the propagation coordinate, asymptotically stabilizing to the linear evolution in $z$ (see Figure 8).

In order to establish the decay law of the primary intensity maximum, the two asymptotic cases must be considered separately. If $\eta=-\left|m\right|\neq 0$, the argument of the Kummer function in (57), as well as that of the exponential, tends to a constant value for large $z$ and so we immediately get $I_{max}^{HG_{II}}\left(z\right)=\left|f_{m}^{HG_{II}}\left[\rho_{max}^{HG_{II}}\left(z\right),\varphi,z;-\left|m\right|\right]\right|^{2}\propto z^{-2}$. On the other hand, when $\eta>-\left|m\right|$ a useful asymptotic expansion of the Kummer function can be employed [78]:

[TABLE]

which applies for large $\left|x\right|$. Remembering that now $\rho_{max}^{HG_{II}}\left(z\right)\propto z$, the local far-field intensity reduces to:

[TABLE]

where expression (58) was used. Since the second term in the sum under square modulus dominates over the first at large $z$, this enables us to obtain $I_{max}^{HG_{II}}\left(z\right)\propto z^{-2}$, which then holds true for every $\eta\geq-\left|m\right|$.

Finally, at fixed radial coordinate and large enough $z$, it is easy to prove that $I_{HG_{II}}\left(z\right)\propto z^{-\eta-2\left|m\right|-2}$, which also describes the power decay in a central region of the beam as a function of the propagation distance, displayed in Figure 9. Of special interest is the case $\eta=-\left|m\right|$, corresponding to the best possible evolution, $z^{-\left|m\right|-2}$.

3.6 Bessel-Gaussian beams

In Section 3.1 it was said that BBs are not square-integrable; however, there exist paraxial solutions of the Helmholtz equation which consists in non-diffracting modes modulated by a Gaussian envelope. Such beams are known as Helmholtz-Gaussian (HzG) waves [114]. Although these modulated solutions lose the divergence-free behavior, they present the advantage of being square-integrable and therefore admit physical realization. In the case of a Bessel configuration, the corresponding HzG solutions are given by the Bessel-Gaussian (BG) beams $\varPsi_{m}^{BG}\left(\rho,\varphi,z,t;k_{\rho}\right)=f_{m}^{BG}\left(\rho,\varphi,z;k_{\rho}\right)\mathrm{e}^{-ikz+i\omega t}$, achievable from expressions (15) and (20) through the choice $a=-i$, $b=1$, $\mu\left(\zeta\right)=1-i\zeta$ and $\phi\left(\zeta\right)=\zeta$, which leads to:

[TABLE]

With reference to Section 3.1:

[TABLE]

Equation (60) is easily reduced to the following:

[TABLE]

corresponding to Bessel equation in the variable $\left(\rho/\mu\right)$. Then, from (15):

[TABLE]

Once again, due to the convention employed for the propagation term, equation (63) is not of the form usually reported in the literature: in order to obtain the standard expression, the whole procedure presented in Section 3 should be repeated starting from the assumption $\psi\left(\rho,\varphi,z\right)=f\left(\rho,\varphi,z\right)\mathrm{e}^{ikz}$.

The intensity distribution of the BG beams is strongly related to the value of the transverse momentum $k_{\rho}$ as well as to the beam waist $w_{0}$ of the Gaussian envelope [111]. The evolution of the radius of the maximum intensity follows the asymptotic law $\rho_{max}^{BG}\left(z\right)\propto z$ (for certain choices of the $k_{\rho}$ parameter a doughnut-shaped profile is also generated when $m=0$, instead of the usual maximum at $\rho=0$). At large distances, $I_{max}^{BG}\left(z\right)=\left|f_{m}^{BG}\left[\rho_{max}^{BG}\left(z\right),\varphi,z;k_{\rho}\right]\right|^{2}\propto z^{-2}$. By exploiting the well-known relation [78]:

[TABLE]

which holds for small $\left|x\right|$, it is possible to infer the decay law of the BG beams intensity in $\rho=\bar{\rho}$ at large $z$: $I_{BG}\left(z\right)\propto z^{-2\left|m\right|-2}$.

To conclude the section, a comparison of the longitudinal intensity profiles of some representative solutions belonging to each of the above described paraxial families is reported in Figure 10. As it can be noticed, the energy distribution of these scalar waves as a function of the propagation coordinate varies considerably from beam to beam.

4 Vortex solutions of the vector wave equation

As shown in the two previous sections, exact and paraxial solutions of the scalar wave equation which are eigenfunctions of the $\hat{J}_{z}$ operator can be easily derived by means of standard mathematical techniques; indeed, such scalar OAM solutions have been deeply explored in the literature. However, an electromagnetic field can hardly be described in terms of a pure scalar function of space and time, since only a full vector treatment allows to take into account the presence of complex polarization structures. Complete solutions of the vector wave equation that can be directly applied to boundary-value problems in electromagnetism have been derived for certain separable systems of cylindrical coordinates and for spherical coordinates: in this context, the electromagnetic field can be resolved into two partial fields, each derivable from a function satisfying the scalar wave equation [74, 131].

Within any closed domain of a homogeneous and isotropic medium with zero conductivity and free-charge density, all vectors characterizing the electromagnetic field satisfy:

[TABLE]

being $\epsilon$ and $\mu$ the inductive capacities of the medium. On using the following relation:

[TABLE]

and assuming a time-harmonic dependence of the form $\exp\left(i\omega t\right)$, equation (65) becomes:

[TABLE]

with $k=\epsilon\mu\omega^{2}$. Equation (67) represents the vector analog of (4). As emphasized by Julius A. Stratton in [131], equation (67) can always be replaced by a system of three scalar equations, but it is only when $\boldsymbol{\psi}$ is expressed in rectangular components that three independent scalar equations are obtained:

[TABLE]

The Laplacian in (68) can of course be written in different coordinate systems, but the vector character of the original equation is inevitably lost. Three independent vector solutions of (67) can instead be built as follows:

[TABLE]

being $g=g\left(x_{1},x_{2},x_{3}\right)$ a scalar function satisfying (4) and $\mathbf{C}$ a unit norm constant vector. The three vector functions above possess several interesting properties: $\mathbf{M}=\mathbf{L}\times\mathbf{C}=k^{-1}\nabla\times\mathbf{N}$, $\mathbf{L}\cdot\mathbf{M}=0$, $\nabla\times\mathbf{L}=0$, $\nabla\cdot\mathbf{L}=-k^{2}g$; moreover, it is easy to prove that both $\mathbf{M}$ and $\mathbf{N}$ are solenoidal and, owing to this and to the fact of being each proportional to the curl of the other, can be employed to represent the electric and magnetic fields.

Since a decomposition of any arbitrary vector wavefunction in terms of (69) is completely general, families of orthogonal vector solutions of the wave equation (65) are derived via (69) from sets of solutions of the scalar wave equation (3) in the various coordinate systems. Also the field polarization is now being considered and care must be taken to properly describe the action of the rotation group generators on the whole vector field.

As can be shown by expanding to first order an infinitesimal rotation about the $z$-axis (see B), the complete form of the $\hat{J}_{z}$ operator for vector fields is provided by:

[TABLE]

where $S_{z}$ corresponds to the third SO(3) generator in three-dimensional matrix representation. When the vector field is resolved into its Cartesian components, the explicit expression for $S_{z}$ is given by $J_{3}$ in (2) and equation (70) reduces to its more common form [132]:

[TABLE]

Let $\left\{g_{m}\right\}=\left\{\varPhi_{m}\left(\varphi\right)u\left(x_{1},x_{3}\right)\right\}$ be a set of OAM scalar functions derived according to Section 2. A general vector solution of (67) can be constructed from $\left\{g_{m}\right\}$ via the superposition:

[TABLE]

where $a_{m}$, $b_{m}$ and $c_{m}$ represent suitable expansion coefficients and the vector fields in the set $\left\{\mathbf{M}_{m},\mathbf{N}_{m},\mathbf{L}_{m}\right\}$ are provided by (69) with $g=g_{m}$. Therefore, the transformation of $\mathbf{\boldsymbol{\psi}}$ under an infinitesimal rotation about the $z$-axis leads back to the action of the $\hat{J}_{z}$ operator on each basis function. In all the five separable coordinate systems for which the $\varphi$-coordinate exists, the following relations hold:

[TABLE]

as explicitly shown in B. The importance of this last result lies in the fact that it is the only requirement of separability of the scalar Helmholtz equation which ensures the vector fields $\left\{\mathbf{M}_{m},\mathbf{N}_{m},\mathbf{L}_{m}\right\}$ built from the OAM set $\left\{g_{m}\right\}$ to be eigenfunctions of the $\hat{J}_{z}$ operator.

From (69), (70) and (73) we understand that the so-derived eigenfunctions cannot be simply considered as vector OAM waves, since their structure is in general more complex, also involving polarization and thus a SAM contribution. Furthermore, for all the three eigenfunctions in (73), it is the vector field components which present the characteristic vortex term $\varPhi_{m}\left(\varphi\right)$, meaning that the phase evolution around the $z$-axis is accompanied by the evolution of the state of polarization. Such eigenfunctions are then named vector vortex beams (VVBs), bearing in mind that in this case the on-axis phase singularity is no longer exclusive and also polarization singularities may be present (see, for instance, [133, 134, 135, 136, 92] and references therein). The above described approach undoubtedly provides a rigorous method for understanding how the VVBs definition can be traced back to general symmetry arguments.

In order to give some instructive examples, let us suppose that the vector potential $\mathbf{A}$ can be represented by an expansion in characteristic vector functions, as in (72). The electric field $\mathbf{E}$ in the Lorentz gauge is then expressed by:

[TABLE]

the usual time-harmonic term $\exp\left(i\omega t\right)$ being implied. Considering a circular cylindrical coordinate system, for which the natural OAM scalar field basis $\left\{g_{m}\right\}$ is provided by BBs, and choosing the $z$-directed unit vector as $\mathbf{C}$ in (69), we get:

[TABLE]

[TABLE]

where all Bessel functions have argument $\left(k_{\rho}\rho\right)$, $\left\{\mathbf{u}_{\rho},\mathbf{u}_{\varphi},\mathbf{u}_{z}\right\}$ represent the cylindrical system unit vectors and (10) has been employed together with some known properties of Bessel functions under the assumption $m\geq 0$, for convenience.

In Figure 11, the norm and polarization plots of (75) are reported for some values of $m$, revealing two apparent deviations from the conventional OAM patterns: when $m=0$, a null is present instead of the usual on-axis intensity maximum, whereas for $m=1$ the central null is replaced by a maximum.

Both deviations are attributable to the mixing of spatial and polarization contributions, which also affects all the other modes, as evidenced by the circumstance that additional factors $\mathrm{e}^{\pm i\varphi}$ appear once (75) is rewritten in terms of the Cartesian unit vectors. It should be noted that the presence of the on-axis intensity null in the $m=0$ case results from a polarization singularity known as V point in singular optics [135].

In [137], Karen Volke-Sepúlveda et al. have shown that equations (74), (75) and (76) make it possible to describe a wide variety of vector solutions in different polarization states. For instance, a radially or azimuthally polarized beam can be obtained by setting $a_{m}=0$ and $b_{m}=i\delta_{m0}E_{0}k/\left(k_{\rho}k_{z}\right)$ or $a_{m}=-\delta_{m0}E_{0}/k_{\rho}$ and $b_{m}=0$ in (74), respectively; here $E_{0}$ represents a constant proportional to the square root of the beam power and with electric field units. Such peculiar vector waves with circular symmetric structure can be seen as free-space transverse magnetic (TM) and transverse electric (TE) modes, as can be deduced from their electric field expressions:

[TABLE]

[TABLE]

where the subscripts $R$ and $A$ stand for radial and azimuthal, respectively, and the time dependence is now explicitly written. Both (77) and (78) are eigenfunctions of the $\hat{J}_{z}$ operator with eigenvalue $m=0$, by construction.

Following the approach presented in [137], right-handed and left-handed circular polarized modes are easily derived from $b_{m}=\mp a_{m}k/k_{z}$ upon proper choice of the $a_{m}$ coefficients:

[TABLE]

being $\left\{\mathbf{u}_{x},\mathbf{u}_{y},\mathbf{u}_{z}\right\}$ the Cartesian unit vectors. In particular, since equation (79) is obtained from the $\left(m\pm 1\right)$ term of (74), we find:

[TABLE]

as can be inferred from (73) or explicitly verified using (71). Whereas the two orthogonal $x$-polarized and $y$-polarized states are simple superpositions of the right-circular and left-circular ones, the most general field is given in terms of the Jones vector Cartesian components $\alpha$ and $\beta$ [5]:

[TABLE]

For $\beta=\pm i\alpha$, equation (81) leads back to (79). Any other choice of the polarization parameters $\alpha$ and $\beta$ in (81) does not give rise to an eigenvector of the $\hat{J}_{z}$ operator, as evidenced by the fact that the circular symmetry of both the beam energy and polarization is broken (see Figure 12). The angular momentum density for the various cylindrical solutions described above has been rigorously calculated in [137], showing that its distribution is indeed not radially symmetric in the case of linearly polarized BBs. On the other hand, the presence of the ratio $\left(k_{\rho}/k_{z}\right)$ in all previous formulas makes it possible to extend the analysis to the paraxial regime $k_{\rho}\ll k_{z}$, with which most of the literature is concerned. Under this approximation, the electric field $\mathbf{E}$ in (81) is found to be eigenvector of $-i\partial/\partial\varphi$, the orbital component of $i\hat{J}_{z}$, and the circular symmetry of the solution norm is restored, as it is shown in Figure 12. Despite the analytic expressions for paraxial VVBs may be obtained more practically by non-separable combinations of spatial and polarization modes [48], the above derivation could be in some circumstances preferable.

Finally, among the various paraxial and non-paraxial VVBs, vector beams presenting cylindrical symmetric states of polarization, often known as cylindrical vector (CV) beams or vortices, have been extensively analyzed due to their interesting properties such as the tighter waist upon focusing (a thorough review of the applications can be found in [138]).

5 Conclusions

In this work, an in-depth overview on the exact and paraxial vortex solutions of the homogeneous wave equation has been presented. These solutions have been interpreted as eigenstates of the $\hat{J}_{z}$ operator, leading to a mixture of polarization and spatial modes in the case of a general vector vortex wave.

Whereas a scalar approach is often convenient to describe the transverse free-space modes of linearly polarized laser radiation in optics, the same does not hold for the radio frequency domain, where a full vector treatment is usually mandatory. In neither cases could the homogeneous wave equation probably be sufficient for the most realistic modeling of the OAM waveforms generated experimentally by means of resonators and antennas, for which only the complete Maxwell equations with real or equivalent source terms can ideally provide the highest level of accuracy. Nevertheless, it is unquestionably true that the homogeneous Helmholtz equation still represents the most affordable way for describing specific free-space solutions as well as for providing general basis sets to expand any arbitrary waveform. For instance, the scalar and vector OAM solutions analyzed in this paper could be easily interpreted as free-space generalizations of resonator and waveguide eigenmodes that are usually found in both laser theory and electromagnetics handbooks. The intensity, phase, polarization and propagation properties of such free-space beams are therefore considered as ideal reference models or, alternatively, as targets to be reproduced experimentally up to a certain degree of approximation via proper synthesis methods. Since many of the exact vortex solutions possess infinite transverse extension, are divergence-free and carry infinite energy, paraxial OAM waves with Gaussian envelope or truncated forms which take into account diffraction and finite-size effects must be preferred in order to ensure a more realistic description.

\ack

The main ideas and results of this work were carried out during my PhD in physics at the University of Torino. Among the various people to whom I am indebted, I would like to especially thank Prof. Roberto Tateo, Prof. Paolo Gambino and Prof. Miguel Onorato, from the Department of Physics, Dr. Rossella Gaffoglio, Prof. Francesco Andriulli and Prof. Giuseppe Vecchi, from the Polytechnic University of Torino, Dr. Assunta De Vita and Eng. Bruno Sacco, from the Centre for Research and Technological Innovation, RAI Radiotelevisione Italiana.

Appendix A Vortex solutions in different coordinate systems

A.1 Spherical coordinates

When expressed in a spherical coordinate system, the scalar Helmholtz equation (4) reads:

[TABLE]

This section is devoted to the search for spherical OAM solutions of the form:

[TABLE]

where $R\left(r\right)$ and $\Theta\left(\theta\right)$ represent functions to be determined and $C$ is a constant. It is immediate to verify that, in order for (83) to satisfy equation (82), $\Theta$ cannot be constant. As a first step, for simplicity, it can be assumed that $R=1$, so that equation (82) reduces to:

[TABLE]

provided that $r\neq 0$. By introducing the parametric coordinate $t=\tan\left(\theta/2\right)$, from which $\cos\theta=\frac{1-t^{2}}{1+t^{2}}$ and $\sin\theta=\frac{2t}{1+t^{2}}$ follow, we get:

[TABLE]

a Fuchsian ordinary differential equation with singularity in $t=0$. Since the roots of the corresponding indicial equation are $\pm\left|m\right|$, the two linearly independent solutions of (85) can be written in the following way:

[TABLE]

From substitution of the power series in (85) we easily find $k,\,d=0$, hence $\Theta^{(1)}\left(t\right)=t^{\left|m\right|}$ and $\Theta^{(2)}\left(t\right)=t^{-\left|m\right|}$, where only the first represents a regular solution in $t=0$. Then, expression (83) reduces to:

[TABLE]

which is not defined for $\theta=\pi$. A regularized version of the previous function over the whole sphere can be provided through the use of both $\Theta^{(1)}$ and $\Theta^{(2)}$, for instance:

[TABLE]

OAM waves of the form (87) and (88) may be called “helico-polar” modes, as they possess helical wavefronts and show a very simple dependence on the polar angle $\theta$. It can be shown that these non-paraxial waves represent square-integrable solutions with respect to the spherical measure $\sin\theta d\theta d\phi$.

In order to prove that the well-known spherical multipoles correspond to generalizations of the above derived waves, the radial term $R\left(r\right)$ in (83) is now reintroduced. Equation (82) becomes:

[TABLE]

where the simplified notation $\dot{X}\left(x\right)=dX\left(x\right)/dx$ has been employed. Then, for $r\neq 0$ we find:

[TABLE]

being $L$ a constant term. To solve the first equation we set $R=\sqrt{r}\mathcal{R}\left(r\right)\mathrm{e}^{ikr}$, which leads to Bessel equation in $kr$:

[TABLE]

and therefore:

[TABLE]

represent the two linearly independent radial solutions, given in terms of Bessel functions of the first and second kind, respectively. After performing the substitution $\cos\theta=t$, the second equation in (90) is easily reduced to:

[TABLE]

which admits solutions in terms of the generalized Legendre polynomials. Finally, by imposing the regularity of $\Theta$ at the boundary points, it can be proven that $L=-\ell\left(\ell+1\right)$ for an integer $\ell\geq\left|m\right|$ and equation (83) acquires the standard form:

[TABLE]

where $b_{\ell}\left(kr\right)$ corresponds to a linear combination of the two independent spherical Bessel functions:

[TABLE]

and $Y_{\ell}^{m}$ are the spherical harmonics.

A.2 Prolate spheroidal coordinates

Let $u\left(\xi,\eta\right)=A\left(\xi\right)B\left(\eta\right)$ be a function that satisfies (6) in the prolate spheroidal coordinate system $\left\{\xi,\varphi,\eta\right\}$. On using the explicit expressions of the metric scale factors [74]:

[TABLE]

where $2\varsigma$ is the distance between the foci of the family of confocal ellipses and hyperbolas, we obtain two forms of spheroidal wave equations:

[TABLE]

with solutions given by [73]:

[TABLE]

where $S_{n}^{\left|m\right|}$ represents a spheroidal wavefunction, $n\geq\left|m\right|$ is an integer and the discrete eigenvalues $\tau=\tau_{n}^{\left|m\right|}$ are analytic functions of $\varsigma^{2}k^{2}$.

A.3 Oblate spheroidal coordinates

In the oblate spheroidal coordinate system $\left\{\xi,\varphi,\eta\right\}$, the metric scale factors read [74]:

[TABLE]

and equation (6) can be separated into:

[TABLE]

whose bounded and single-valued solutions are written as [73]:

[TABLE]

A.4 Parabolic coordinates

In parabolic coordinates $\left\{\mu,\varphi,\nu\right\}$, with metric scale factors [74]:

[TABLE]

$u\left(\mu,\nu\right)=U\left(\mu\right)V\left(\nu\right)$ and equation (6) reduces to:

[TABLE]

The separated solutions take the form [73]:

[TABLE]

where ${}_{1}F_{1}$ is the Kummer hypergeometric confluent function.

Appendix B Angular momentum operator

Let $\boldsymbol{\psi}\left(x_{1},\varphi,x_{3}\right)$ be an arbitrary vector function in one of the five coordinate systems $\left\{x_{1},x_{2}=\varphi,x_{3}\right\}$ for which equation (6) separates. Under infinitesimal rotation by an angle $\varepsilon$ about the $z$-axis, the function $\boldsymbol{\psi}$ transforms as:

[TABLE]

where $\left\{\mathbf{u}_{1},\mathbf{u}_{2}=\mathbf{u}_{\varphi},\mathbf{u}_{3}\right\}$ correspond to the standard unit vectors of the considered system and a first order Taylor expansion in $\varepsilon$ has been performed. Equation (105) can be rewritten in the following form:

[TABLE]

On using some basic results of differential geometry, it can be shown that [74]:

[TABLE]

with $i\in\left\{1,2,3\right\}$ and $h_{2}=h_{\varphi}$. For the considered class of coordinate systems, none of the three metric scale factors depends on the azimuthal variable and, by means of formula (107), expression (106) reduces to:

[TABLE]

where $S_{z}$ represents a matrix operator which acts on the fundamental column vectors:

[TABLE]

and whose explicit form depends on the coordinate system:

[TABLE]

By direct calculation, we find:

[TABLE]

where $B$ represents the basis change matrix from the curvilinear coordinate system $\left\{x_{1},\varphi,x_{3}\right\}$ to the Cartesian one, $\left\{x,y,z\right\}$. Therefore, $S_{z}$ exactly corresponds to the third SO(3) generator in three-dimensional matrix representation.

Since any arbitrary rotation of the vector function $\boldsymbol{\psi}$ can be expressed through the rotation operator, the following relation holds:

[TABLE]

If we compare equations (108) and (112), we immediately obtain:

[TABLE]

which is the sought expression for the third component of the vector angular momentum operator (in $\hbar$ units).

Let now $\boldsymbol{\psi}=\boldsymbol{\psi}_{m}$ be one of the three vector functions $\left\{\mathbf{M}_{m},\mathbf{N}_{m},\mathbf{L}_{m}\right\}$ defined in (69), with $g=g_{m}=\varPhi_{m}\left(\varphi\right)u\left(x_{1},x_{3}\right)$ representing a scalar function derived according to Section 2. Owing to the fact that the three metric scale factors do not depend on the azimuthal variable in the five considered systems, the partial derivative $\partial/\partial\varphi$ commutes with each component of the differential operators appearing in (69), as can be shown explicitly by resorting to the formulas [74]:

[TABLE]

Making use of this property and of equation (105), it is easy to prove that:

[TABLE]

Then, from comparison with (112):

[TABLE]