Bessel Beams: A Unified Perspective

The literature on Bessel beams is both fragmented and incomplete. We present here a unified and completed perspective of such beams, irrespective to their orbital angular momentum (OAM) -- zero, integer or noninteger -- and mode -- scalar or vectorial, and LSE/LSM or TE/TM in the latter case. The unification is based on the integral superposition of constituent waves along the angular-spectrum cone of the beam, and allows to describe, compute, relate, and implement all the Bessel beams in a universal fashion. The paper first establishes the integral superposition theory. Then, it demonstrates the previously unreported existence of noninteger-OAM TE/TM Bessel beams, compares the LSE/LSM and TE/TM modes, and establishes useful mathematical relations between them. It also provides an original description of the position of the noninteger-OAM singularity in terms of the initial phase of the constituent waves. Finally, it introduces a general technique to generate Bessel beams by an adequate superposition of properly tuned sources. This global perspective and theoretical complement may open up new avenues in applications such as spectroscopy, microscopy, and optical/quantum force manipulations.

Electromagnetic Bessel beams represent a fundamental form of structured light. They are localized waves that carry orbital angular momentum (OAM) along their propagation axis. Localized waves are waves that propagate without spatial dispersion (or diffraction) and without temporal (or chromatic) dispersion. Their energy is thus uniformly confined and invariant perpendicular to and along, respectively, their direction of propagation. The OAM is a beam property whose macroscopic manifestation is an isophase surface that has the form of a vortex along the axis of the beam. It may be integer or noninteger. In the former case, the wave has the transverse phase dependence e inφ (n ∈ Z), corresponding to an OAM of n per photon 1 , while in the latter case the wave is made of a superposition of integer-OAM waves that combine so as to produce a noninteger-OAM per photon [2]. This property association of localization and OAM confers to Bessel beams specific capabilities for manipulating light that may be exploited in diverse applications, such as nanoparticle guiding [3] and trapping [4], spectroscopy [5], microscopy [6], and quantum key distribution [7]. Figure 1 shows a general classification of Bessel beams that pertain to the sequel of the paper.
Bessel beams are the simplest form of light OAM after the Laguerre-Gauss beams and the most studied localized waves. They are monochromatic beams with a transverse amplitude pattern that follows Bessel functions of the first kind, J n (αρ), multiplied by the phase function e inφ (n ∈ Z), or combinations of such waves in the noninteger OAM case. Their simplest forms are the scalar Bessel beams, also existing in acoustics and restricted to the paraxial approximation in optics. Such beams were discovered by Durmin in 1987 [8] for n = 0 (no-OAM), while their integer-OAM carrying versions were reported shortly thereafter [9]. Noninteger-OAM scalar Bessel beams were reported only 15 years later [10]. In the case of electromagnetic waves, such as light, Bessel beams are vectorial, and may be either LSE/LSM 2 or TE/TM, as mentioned in Fig. 1. Vectorial Bessel beams, the only exact forms of electromagnetic Bessel beams, were introduced in integer-OAM form in [13] and generalized to the noninteger case in [14]. Bessel beams of various complexities have been generated by illuminating an axicon with a Laguerre-Gaussian beam [15], by a single helical axicon (azimuthal phase plate) [16], by a spatial light modulator [17], [10], and by a metasurface [18]. Unfortunately, the literature on Bessel beams is very fragmented. It is composed of different theoretical perspectives that are often restricted to a specific type of Bessel beam or deprived of physical insight. Moreover, the circular cylindrical Bessel modes have received surprisingly little attention to date, with their noninteger version still remaining on open question to date, despite their fundamental character. Finally, the aforementioned implementation techniques are as specific of the theoretical approaches, cumbersome or/and inefficient. This paper resolves this double issue of fragmentation and incompleteness by presenting a unified and completed perspective of Bessel beams. The unification is based on the integral construction of the beam along its angular-spectrum cone, and allows to describe, compute, relate, and implement all the Bessel beams (see Fig. 1) in an equally simple fashion. The paper first establishes the theory of this approach. Then, it demonstrates the existence of noninteger-OAM TE/TM Bessel beams in addition to their LSE/LSM counterparts. Finally, it introduces a general technique to generate Bessel beams by an adequate superposition of properly tuned sources.
A Bessel beam may be generally described by the scalar function of space (x, y, z) and time (t) where ψ ν (φ G ) represents a continuous set of waves that propagate at the angle φ G along a circular cone towards its apex, so as to form a transverse interference pattern in the form of a Bessel function, as illustrated in Fig. 2, and where ν is equal to the OAM of the beam when it is integer or half-integer (see Supp. Mat. A). These waves have the mathematical form with the oblique wave vector where so that α 2 + β 2 = k 2 0 = (ω/c) 2 , and δ is called the axicon angle, and with the linear phase In these relations, w(ξ) is the transverse apodization of each wave with respect to its propagation direction, k(φ G ), with ξ being the radial variable of the corresponding local coordinate system, frac(·) is the fractional part function [19], and φ G,0 is related to the position of the phase singularity for ν ∈ R\Z, as will be seen later. Note that Eq. (2d) represents a sawtooth function of φ G that reduces to γ ν (φ G ) = νφ G for φ G ∈ [0, 2π[ and φ G,0 = 0. Practically, the waves ψ ν (φ G ) must be spatially limited (e.g. Gaussian cross-sectional apodization w(ξ)), but they will be  initially considered as plane waves (w(ξ) = const.) 3 because, as we shall see, an analytical solution can be derived for Eq. (1) in this case 4 . Moreover, their azimuthal separation (∆φ G ) is infinitesimally small, so that they effectively merge into the azimuthal continuum corresponding to the integral (1). Finally, these waves may be temporarily considered as a scalar, but they will later be seen to represent any of the components of fully vectorial electromagnetic waves.
The Bessel beam superposition in Eq. (1) is plotted in Fig. 3 using 20 equi-spaced plane waves for ψ ν (φ G ) in Eq. (2), with the top row showing the superposition of the maxima and minima, and the bottom row plotting the superposition of the actual waves with continuous sinusoidal gradients. This figure shows that a perfectly smooth Bessel pattern is obtained around the axis of the beam with a restricted number of constituent waves 5 , and illustrates the increasing structural complexity of the Bessel beam without OAM (ν = 0), with integer OAM (ν = 1) and with non-integer OAM (ν = 1.5).
Substituting Eq. (2) with w(ξ) = A PW (const.) into Eq. (1), simplifying the resulting integral (see Supp. Mat. C), and decomposing the field into its transverse-dependence and longitudinal/temporal-dependence parts as U ν (ρ, φ, z, t) = 3 Note that the beam is still localized in this case due to the decreasing envelope of the Bessel interference pattern. 4 Eq. (1) may then be alternatively expressed as the inverse Fourier transform Uν (ρ, φ, z) Mat. B). 5 The top row of Fig. 3 shows that this radial effect is due to the radially decreasing density of the constituent waves. U ν,⊥ (ρ, φ)e i(βz−ωt) , yields then (3) For ν = n ∈ Z, the integral in Eq. (3) has a tabulated closed-form primitive [20], leading to U n,⊥ = 2πi −n e −inφG,0 A PW J n (αρ)e inφ , which is the conventional integer-OAM Bessel solution for circular-cylindrical problems. In contrast, for ν ∈ R\Z, the integral does not admit a simple primitive, and we must devise a strategy to lift this restriction so as to find the most general beam solution. This can be accomplished by the following procedure. First, we replace the generally non-periodic (ν ∈ R\Z) complex exponential function e iνφ in Eq. (3) by its expansion in terms of the complete and orthogonal set of periodic (ν = m ∈ Z) complex exponential functions e imφ (see Supp. Mat. D) [21], i.e., Then, we substitute Eq. (4) into Eq. (3), which leads to (5) Finally, we eliminate the integral by applying the Bessel identity 2π 0 e −iαρ cos(φ−φ ) e imφ dφ = 2πi −m J m (αρ)e imφ [22], and find the general solution to Eq. (5) in the closed form with the complex weighting distribution Equation (6) contains the integer-OAM Bessel solution, since ν = n ∈ Z transforms Eq. (6b) into A BB m = 2πi −m e −inφG,0 A PW δ mn , which reduces the sum in Eq. (6a) to the single term U n,⊥ = 2πi −n e −inφG,0 A PW J n (αρ)e inφ . But it also contains noninteger-OAM (ν ∈ R\Z) solutions, where the satisfaction of the circular periodic boundary condition is realized by a superposition of integer-OAM Bessel waves with proper phase (e imφ ) and weighting coefficients (A BB m (ν, φ G,0 )). In this noninteger-OAM case, the sum in Eq. (6a) must be practically truncated to an integer m = ±M that is large enough to provide a satisfactory approximation of the Bessel beam.
Equation (6a) reveals that the parameter α of the cone in Fig. 2 corresponds to the spatial frequency of the Bessel pattern. Since this parameter is proportional to the axicon angle, δ, according to Eq. (2c), we find that increasing the aperture of the cone in the integral construction of Eq. (1) compresses the Bessel ring pattern towards the axis of the beam. Figure 4 plots the magnitude and phase of the Bessel beam given by Eq. (6) for different integer and noninteger OAMs. The cases ν = 0, 1 and 1.5 correspond to the instantaneous field plots in Fig. 3. The OAM-less beam ν = 0 has simultaneously azimuthally symmetric magnitude and phase. The integer beams ν = n ∈ N have azimuthally symmetric magnitudes but asymmetric phase (OAM). The beams ν ∈ R\Z have simultaneously asymmetric magnitude and phase. The parameter φ G,0 , which is here 0, corresponds to a dummy initial phase of the integer-OAM and the position of the discontinuity of the noninteger-OAM in the individual waves, with increasing φ G,0 clockwise rotating the asymmetric magnitude of the noninteger-OAM pattern (see Supp. Mat. E). The general 6 scalar Bessel solutions described above are restricted to acoustic waves, quantum waves, and vectorial waves under special conditions such as the paraxial condition (δ π/2) in the electromagnetic case. On the other hand, they fail to describe Bessel beams with a large axicon aperture (angle δ in Fig. 2), which are relevant to applications such as microscopy and optical force manipulations. Therefore, we extend here the previous scalar generalization to the vectorial case.
For the scalar case, we have established two alternative solutions, the integral solution of Eq. (1) and the analytical solution of Eq. (6). In the present extension to the vectorial case, we shall restrict our treatment to the integral approach 7 , because it provides more insight into the physical nature of the beam and because it will constitute the basis for the practical implementation to be discussed later. We shall still assume a plane wave construction (w(ξ) = A PW (const.) in Eq. (2b)), for simplicity.
As indicated in Fig. 1, the vectorial Bessel beams can be LSE i /LSM i with i = {x, y} or TE z /TM z , where the subscript denotes the field component that is zero. In the former case, the electric/magnetic transverse field component that is nonzero is set as the scalar Bessel beam solution U ν in Eq. (1), while in the latter case, it is the magnetic/electric longitudinal component of the field that is set as U ν , and the other components are found from Eq. and where η is the impedance of free-space. The most striking difference between the LSE/LSM and TE/TM beams resides in the simplest feature of their respective constituent waves.
The former have a linear transverse polarization, while the latter have of a constant transverse magnitude.
A detailed investigation of these solutions reveals that the axicon angle (δ) distinctly affects the LSE/LSM and TE/TM modes (see Supp. Mat. G). In both cases, increasing δ compresses the Bessel ring pattern towards the axis of the beam; however, this variation also breaks the symmetry of the transverse LSE/LSM patterns, even for ν = 0, whereas it leaves the TE/TM pattern azimuthally symmetric. It is also interesting to note that, for a small axicon angle, i.e., δ π/2, the LSE y mode can essentially reduces to its E x and H y components, similarly to the scalar form. Figures 5 and 6 depict the LSE y and TM z Bessel beams of global order ν = 1.5 corresponding to the solutions of Eqs. (7) and (8), respectively. The vectorial field distributions plotted in the panels (a) and (b) of the two figures represent samples of the constituting waves of the integral construction of the beam (Fig. 2). Their strong vectorial nature starkly contrasts with the configuration of the scalar solution, except for the LSE y case in the aforementioned axicon limit (δ π/2). Note that the electric field of the constituent waves of the LSE y mode is linearly polarized in the x-direction, while that of the TM z modes is radially polarized. Nonzero and noninteger ν vectorial modes are obtained from their fundamental counterpart by simply setting the ν parameter in the initial phase of the constituent waves, i.e.,   7)) and TM + z (Eq. (7)) modes with δ = 25 • .
The LSE/LSM and TM/TE electromagnetic vectorial Bessel beams are related by the following relations, which may be easily verified upon comparing Eqs. (7) and (8): (9c) where the second subscripts (ν and ν ± 1) correspond, as usual, to the global OAM ν in Eq. (2a). Note the OAM conservation between the LSE/LSM and TM/TE modes in each of these relations, and the interesting mediation of the axicon angle. Note also that, in the paraxial approximation (δ π/2), the TM/TE beams, with their complex cylindrical (radial/azimuthal) polarizations, be realized superpositing two transversally linearly-polarizaed LSE/LSM beams, according to Eq. (9a).
Several techniques have been proposed for generating Bessel beams experimentally. The main ones are axicon lenses illuminated by a Laguerre-Gauss beam [15], spatial light modulators [17], open circular waveguides with selectively excited modes [24], antenna arrays with proper phase feeding network [25], metasurfaces illuminated by plane waves [18], and 2D circular leaky-wave antennas [26]. Unfortunately, these techniques are restricted to simple beams, excessively complex to implement, bulky and expensive, or suffering from poor efficiency.
The unified integral formulation presented in this paper ( Fig. 2 with Eq. (1) for the scalar case, and Eqs. (7) and (8) for the vectorial case) naturally points to a generation technique that is immune of these issues and that offers in addition a universal implementation framework. Indeed, circularly distributing a set of sources with the phases, amplitudes and polarizations of the derived modal field solutions (e.g., top panels of Figs. (5) and (6)) would exactly and efficiently produce the corresponding Bessel beams, irrespective to their order or complexity.
Specifically, the integral-formulation generation technique consists in the following design steps: 1) select a sufficient number of sources (N ) to properly sample the desired OAM according to the Nyquist criterion, 2) determine an appropriate beam apodization (w(ξ)) for each of the constitutent waves to be radiated by these sources, and 3) adequately set the phase, magnitude and polarization of each of the sources, and orient them so as to launch the constituent waves along a cone with the selected axicon angle (δ). This is mathematically expressed by the formula where w(φ G ) is the apodization of the constituent waves, E PW (φ G ) is their plane-wave modal field solution (e.g., Eqs. (7) or (8)), and ∆φ G = 2π/N . In the case of a (typical) Gaussian apodization, we have where w 0 is the waist of the beam, and (x • (φ G ), y • (φ G )) represents the local conical coordinates which are related to the radial conical coordinate Note that apodization of the plane wave E PW by the function w(ξ) results into a localization of the beam in a restricted of extent L = w 0 / sin(δ) about the center of the cone at (z = 0). Moreover, the discretization of the integral induces a distortion of the Bessel pattern, which grows with the distance from the axis of the beam, as previously explained, so that N may have to be increased to provide a satisfactory beam approximation across the transverse area of interest. Figure 8 depicts the experimental implementation of the integral-formulation Bessel beam generation. Figure 8(a) represents a direct incarnation of this formulation, which consists of a circular array of laser beams with proper magnitudes, phases and polarizations, as illustrated in Fig. 8(b). Such an implementation, involving N independent lasers with respective magnitude, phase and polarization controls, is quite complex and cumbersome. Fortunately, recent advances in metasurface technology suggests the much more practical implementation shown in Fig. 8(c). Indeed, this metasurfacebased Bessel beam generator requires only one laser source, while being ideally compact and inexpensive.
The metasurface required in the implementation of the Bessel beam generator depicted in Fig. 8(c) can be easily realized using latest metasurface synthesis techniques [27]. The simplest implementation strategy would consist in cascading metasurfaces that separately tailor the amplitude, the phase, the inclination, and the polarization of the incident wave in the transverse plane of the system. Specifically, assuming a linearily polarized incident wave, such a design would then consist in three cascaded metasurfaces. Two of these metasurfaces would be common to the LSE/LSM and TE/TM cases, with one metasurface providing the required azimuthal phase distribution via azimuthal sectors made of particles inducing progressive transmission delays, and the other providing the required conical inclination via a constant radial phase gradient. In contrast, third metasurface would be different for the LSE/LSM and TE/TM cases. In the former case, given the linear transverse polarization (Eq. (7) and Fig. 5), there is no required polarization processing, and the third metasurface needs to provide the proper transverse magnitude distribution, which can be accomplished with dissipative particles, while in the latter case, given the constant transverse magnitude (Eq. (8) and Fig. 6), there is no required magnitude processing, and the third metasurface needs to provide the proper transverse polarization distribution, which can be accomplished with birefringent particles. We have presented a unified perspective of Bessel beams of arbitrary OAM (zero, integer and noninteger) and nature (scalar, LSE/LSM and TE/TM) based on an integral formulation, and deduced from this formulation a universal and efficient generation technique. The proposed formulation may be extended to other conical beams, such as Weber and Mathieu beams, both to increase the insight into their characteristics and to facilitate their generation. This global perspective opens up new horizons in structured light for a variety of applications, such as spectroscopy, microscopy, and optical/quantum force manipulations.

SUPPLEMENTARY MATERIAL A. OAM PER PHOTON FOR A BESSEL BEAM OF GENERAL ORDER
The total time-averaged angular momentum (AM) per photon for a monochromatic vectorial beam with propagation axiŝ z is [28] where g = Re {D * × B} /2 is the timeaveraged linear electromagnetic momentum density, u = Re {E * · D + B * · H} /4 is the time-averaged energy density, and V is a cylindrical volume enclosing the beam.
Wave beams, including Bessel beams, are typically satisfactorily described in terms of their paraxial approximation. In this regime, the fields D, B, H can be expressed, from the Maxwell equations and from the constitutive equations, in terms of the field E, and the longitudinal part of E, E z , can be expressed in terms of its transverse part, E T . Expressing then g and u in Eq. (A1) in terms of E T , neglecting the terms in 1/k 2 and 1/k 4 [28], and integrating the remaining terms by part yields which splits into where L z is the orbital AM (OAM) and S z is the spin AM (SAM).
In the case of a Bessel beam with transverse linear polarization (case of LSE and LSM modes, but not TE and TM modes, in the paper), i.e., E T = Ux, Eq. (A2b) reduces to where the radial integrals have been truncated to a finite radius, R, to account for the finite energy of a practical beam. Figure S1 plots, using Eq. (A3), the OAM per photon for a Bessel beam [U = U ν in Eq. (1) with ψ ν (φ G ) in Eq. (2)] for R = 25 m. As can be seen, the OAM is a continuous and wiggly function of the parameter ν, which coincides with the linear function L = ν only at integer and half-integer values of ν.

SUPPLEMENTARY MATERIAL B. VERIFICATION OF THE CONICAL SPATIAL ANGULAR SPECTRUM
We shall show here that the conical spatial angular spectrum where the ± solutions in each of the two expressions correspond to the ranges 0 ≤ φ G,0 ≤ π (+ sign) and π ≤ φ G,0 < 2π (− sign), which appears in the the spatial inverse Fourier transform corresponds to the integral representation of Eq. (1) with Eq. (2), or Eq. (3).

SUPPLEMENTARY MATERIAL C. DERIVATION OF EQ. (3)
Substituting Eq. (2) with w(ξ) = A PW (const.) into Eq. (1) yields where, assuming 0 < φ G < 2π, Inserting the Eq. (C2) into Eq. (C1) splits the integral of the former as (C4) Finally, the two integrals in this relation can be grouped into a single one via the definitions φ 1 := φ and φ G := φ , which results into (D2) Inverting the order of the integral and the sum in the righthandside term and using the orthogonality property of the integer-order comlex integral function gives then φG,0+2π φG,0 which integrates into where the dummy prime has been dropped. Substituting Eq. (D4) into Eq. (D1) finally yields the explicit expansion The physical role of the phase term φ G,0 in Eq. (2d) and Eq. (6) is not trivial. This section clarifies this issue, using two distinct and complementary strategies: 1) analyzing the problem in the spatial-frequency spectral (k) domain, where it is simpler, and 2) inspecting a beam made of truncated constituent waves in a cross-sectional plane where these waves weakly overlap and therefore only locally (i.e., only with nearest neighbors, see Fig. 2) interfere, which might be considered as a spatial alternative to 1).
The spatial-frequency spectrum solution of the general solution in Eq. (6) is (see Supp. Mat. B) (E1) Equation (E1) describes a cone with different initial phases in the k-space. If ν = n ∈ Z, the term e in2π = 1 in the top expression disappears so that this expression reduces to the second; as a result, the initial phase φ G,0 ± π does not depend any more on k φ , and represents therefore an dummy initial phase that can be set to zero, which reduces the spectrum tõ the situation is obviously more complex, with a phase (k φ ) discontinuity of 2πν appearing in Eq. (E1), specifically (E2) Thus, the parameter φ G,0 corresponds then the position of the spectral phase discontinuity as k φ,0 = φ G,0 + π. Figure S2 illustrates these results, Moving now on to the second strategy, Fig. S3 plots the effect of φ G,0 on the direct fields associated with Bessel beam with circular-cylindrically truncated constituting waves. Here, the parameter φ G,0 , has also an impact of the phase, with the neighboring waves canceling each other at the angle φ = φ G,0 (z < 0).  Fig. S3. Complex amplitude and phase of the direct fields associated with a Bessel beam with Gaussian circular-cylindrically truncated constituting in a plane of overlap only the nearest neighbours (z G < z < 0) and in the overlaping plane (z = 0) (same parameters as in Fig. S2). Figures S4 and S5 show the effect of the phase discontinuity, as the Figures S2 and S3 respectively, for a given φ G,0 = 0 and different global order ν. Figure S5 shows how the amplitude (of neighboring waves) opening is affected by the phase discontinuity 2πν.

A. Maxwell Equations for a Plane Wave
The free-space Maxwell equations for a harmonic plane wave propagating in the direction k read where η = µ 0 / 0 is the impedance of free-space and k = ω √ µ 0 0 is the wavenumber of free-space, with 0 = 8.854 · 10 −12 F/m and µ 0 = 1.257 · 10 −6 N/A 2 .

B. Construction of the LSE/LSM Modes
Let us start with the LSE y (E y = 0) modes. In this case, we have from which we find, upon solving Eq. (F1c) for E PW z , which becomes, upon setting E PW x (φ G ) = ψ ν (φ G ) (transverse field in the scalar case) in Eq. (2a) and using Eq. (2b), The corresponding components of H PW are found by inserting Eqs. (F1) and (F2b) into Eq. (F1a), and successively solving for the different components, which results into The LSE y Bessel beam corresponding to the integral form of Eq. (1) is then formed, upon grouping the previous results, by the field components The LSE x Bessel beam field components are similarly found as For the sake of completeness, we finally give the LSM x and LSM y modes, which are respectively and

C. Construction of the TM/TE Modes
Let us start with the TM z (H z = 0) modes. In this case, we have H PW z (φ G ) = 0.
Additionally, we find from the z-component of Eq. (F1a) the following relation between the transverse components of the electric field: (F3d) The corresponding components of H PW are found by inserting Eqs. (F3c) and (F3d) into Eq. (F1a), and successively solving for the different components, which results into The TM z Bessel beam corresponding to the integral form of Eq. (1) is then formed, upon grouping the previous results, by the field components The TE z Bessel beam integral field components are similarly found as Essentially, increasing δ increases the transverse spatial frequency α (see Eqs. (2b) and (6a)) of the Bessel waveform of each of the field components and hence compresses the Bessel ring pattern toward the beam axis. Figure S8 shows the effect of δ on the norm of the timeaveraged Poynting vector and time-averaged energy density u : Increasing δ breaks the azimuthal symmetry of the Poynting vector and the energy density of the LSE y mode, but not that of the TM z mode, as could be expected from the previous resutls.