Generalized ray optics and orbital angular momentum carrying beams

In classical optics the Wolf function is the natural analogue of the quantum Wigner function and like the latter it may be negative in some regions. We discuss the implications this negativity has on the generalized ray interpretation of free-space paraxial wave evolution. Important examples include two classes of beams carrying optical orbital angular momentum—Laguerre–Gaussian (LG) and Bessel beams. We formulate their defining eigenfunction properties as phase–space symmetries of their Wolf functions, whose analytical form is shown, and discuss their interpretation in the ray picture. By moving to a more general picture of partly coherent fields, we find that new solutions displaying the same symmetries appear. In particular, we find that mixtures of Gaussian beams (thus fully describable using classical ray optics) can mimic the basic properties of LG beams without the need for negativity, and are not restricted to quantized values of angular momentum. The quantization of both the l and p parameters and negativity of the Wolf function are both inevitable and, indeed, arise naturally when a requirement on the purity of the solution is added. This work is supplemented by a set of computer animations, graphically illustrating the interpretative aspects of the described model.


Introduction
We know that most phenomena of light can only be explained when its wave nature is taken into account. Within the classical context, this has effectively eliminated attempts at a broader adoption of Newton's corpuscular description [1]. Nevertheless, geometric ray optics remains a good compromise between clarity and exactness in many real-world situations. Historically, its domain had been restricted to the limit of short wavelength, applicable to a good extent to many important fields including photography, microscopy, and telescopy, but failing at scales where interference phenomena can no longer be neglected.
The concept of generalized ray description of light waves, providing an exact mathematical model of light propagation and detection accounting even for interference phenomena, but keeping the intuitiveness of a geometric picture, dates back to the 1970ʼs. First attempts can be identified in [2], after which the topic has been broadly-and to a good extent independently-developed by Bastiaans [3][4][5] and Sudarshan [6][7][8]. Simon [9] has noted that the ray picture becomes particularly plausible and well-behaved in the paraxial approximation, which was an ad hoc initial assumption in [2].
In order to successfully extend the ray picture of light to describe interference phenomena, some assumptions must be relaxed. Most notably, one must allow to consider rays carrying a negative intensity. These form the basis for describing destructive interference phenomena without the notion of a phase.
In this work, we apply the above ideas to the case of paraxial beams carrying optical orbital angular momentum (OAM) [10]. Laguerre-Gaussian (LG) beams, in particular, have been the target of a broad interest lately in connection to their use in, for example, optical tweezers and optical wrenches [11], multimode optical communication [12], or surface analysis [13], among many others, along with a recent development in the accessibility of computer-generated holograms as a primary means of their fully customizable generation [14,15]. Many important generalizations of practical applicability have been thoroughly studied, for example, beams with optical vortices displaced from their centroid [16]. The topic in turn is a part of a much broader field, extending also into the quantum optical domain where the spatial distributions are reflected in the sense of an abstract harmonic oscillator representation, allowing for methods similar to those described here to be employed [17,18].
There certainly remains space for research in some of the theoretical aspects of OAM-carrying beams in their different representations. For example, a picture of a LG beam as formed by a sheaf of co-rotating rays had been repeatedly used in literature (see, for example, [19,20], or in electron optics, [21]) inspired by the approximate shape of the integral curves of the Poynting vector [22,23], only reaching this behaviour asymptotically for both azimuthal and radial indices simultaneously large [24]. This picture does not address the curvature of Poynting vector lines nor the divergence of its axial component near beam axis [25], the latter leading to a theoretical paradox [26]. By providing a complete generalized ray description of LG and Bessel beams, along with remarks on their interpretation, we aim to strengthen the methods currently in use as well as to resolve the aforementioned paradox and to provide new insights into the topic.

Paraxial wave equation and the Wolf function
Consider a weakly divergent monochromatic beam propagating along the z + axis such that in its decomposition into plane waves, most of the intensity is concentrated in a close neighbourhood of k k 0, 0, , Let us ignore its polarization (or select one fixed polarization direction) and let the beam be described by a scalar function x y z , , ( ) y obeying the Helmholtz equation Due to the paraxial structure of the beam, it can be assumed that most of the evolution along the z axis will be constituted by a linear phase ramp of e . kz i Extracting this factor, let us define , , e , kz i ( ) ≔ ( ) y and represent the assumption of paraxiality by [27] z u x y z k z u x y z , , , On rewriting (1) in terms of u x y z , , ( ) and applying neglections conforming to (2), we obtain the paraxial wave equation [27,28,30]: In contrast with (1), this is a first-order evolution equation in z and thus only requires the specification of u x y z , , ( ) at a single plane z z 0 = as the initial condition.
There is an exact formal agreement between the form of (3) and the (2 + 1)-dimensional free particle Schrödinger equation The agreement is complete if we identify the evolution parameter t with z and k with the Compton wave number m/ (both using c = 1). Many quantum mechanical methods can then be directly translated into paraxial wave optical domain. One of the notable concepts is the Wigner function [29], , , e d d , 4 x y x p y p x y In the following, we shall assume both functions to be normalized in the L 1 norm (although it is more natural for the latter to integrate to the full cross-sectional intensity). This is automatically satisfied if x y t , , ) are L 2 -normalized, respectively. Further, we shall drop the reference to the evolution parameter z in W x y k k z , , , x y ( |)where its value has been fixed and no ambiguity can occur. Equations (4) or (5) assign a real-valued phase-space function to any given wave function x y t , , or a paraxial beam cross section u x y z , , , ( ) respectively. In many respects these functions can be interpreted as a probability (in the case of W qm ) or intensity (for W) distribution in a restricted sense. Notably, x y x y x y x y where ũ denotes the two-dimensional Fourier transform of the wave function u, among many other properties.
Of particular interest in phase-space formulation of quantum mechanics is that expectation values of observables can be computed averaging their phase-space representations using W qm as the weight. This representation, called the (Weyl's) symbol f A of an operator A, is in this case obtained by Note in particular that W qm is, up to a prefactor of 2 , 2 ( )  pthe symbol of the density matrix of the state. A similar theorem holds also for the Wolf function with obvious interpretational differences, primarily expectation values becoming ensemble averages. The mapping of the two formalisms is achieved by replacing ÿ by unity, time by the z coordinate, and components of momentum by those of wave vector perpendicular to the direction of propagation. In this respect, we will also speak of density operator z ( ) r corresponding to a wavefront u x y z , , ( ) x y z x y u x y z u x y z x y z x y z , , , representing the equal-time second-order cross-correlation function of the monochromatic field [30], and we will carry the notation of a symbol of an operator to paraxial optical systems. The evolution equation (3) can be rewritten for the Wolf function as follows: and explicitly solved as x y This result can be interpreted in a direct correspondence to the phase-space intensity density picture: the local density is preserved along straight trajectories of the slope k k k , , x y ( )in the x y z , , ( ) space, keeping k x and k y constant. Thus, in terms of measurable physical quantities, any solution to the paraxial wave equation (3) can be decomposed into a set of ideal, non-interacting, linear rays going through every point of the cross section in every direction and carrying a positive, zero, or negative contribution to local intensity.
The possibility of negative values is an inherent property of the Wolf function and the only solutions of (3) with a completely non-negative Wolf function are Gaussian beams [31]. It is important to emphasize at this point that negative rays can never be directly observed. This is precluded by marginalization (6) in position or transverse k detection, which always yields non-negative values, or in general by the positivity of measurement operators when more intricate measurements are performed. In an intuitive understanding, positive rays are always stronger in total intensity in any phase-space cell of x k y k 1, 1 x y   D D D D than negative ones and trying to separate the latter would inevitably overcome the diffraction limit, compromising the validity of the equations and approximations used.
LG beams A LG mode is a prototypical example of an optical beam carrying OAM. Its form is stable under free-space propagation (only the scale of the intensity profile in a cross section changes) and each photon carries a constant integer number of OAM quanta, L l .  = Here and in the following, unless explicitly stated otherwise, we shall always assume that coordinates are chosen such that z coincides with the beam axis and z = 0 is the point of the narrowest cross section, the beam waist. An LG beam is then completely specified by l ,  Î an integer 'radial' parameter p 0  , and the width at the waist w 0 , and often denoted LG .
l p , The mathematical form of its wavefront in the scalar theory is [10] u r z w z denotes the beam width at z, reducing to the beam waist 3 w 0 at z = 0, the radius of curvature, and This function can be found as a joint eigenfunction of the differential angular momentum operator corresponding to eigenvalues l and p l 2 1 , | | + + respectively.
As both operators are quadratic in position and momenta observables, their respective eigenvector equations correspond to phase-space symmetries of the Wolf functions of the LG beams: in general, for any operator Â satisfying this condition, the Wolf function of the commutator A Âr r corresponds, up to a constant factor, to the Poisson bracket of the Wolf function assigned to ρ and the symbol of the operator A. For u x y z , , ( ) satisfying the eigenvector equation the density operator (7) commutes with A, and therefore the Wolf function (5) satisfies The symbols of operators (11) and (12) are Corresponding here to an intensity drop by a factor of e 1 in the Gaussian case l p 0. = = The more usual e The Wolf function of a LG beam in the beam waist can be computed using further properties of these symmetries, as demonstrated very elegantly by Simon and Agarwal in [32], or using a set of ladder operators for the l and p numbers, an approach taken by Vanvalkenburgh [33]. The result can be written in a compact form, manifestly symmetric with respect to (14) and (15), The coordinates μ and ν, which both range from zero to infinity, can be supplemented by two angular coordinates , 0,2 to parametrize the full phase-space:  f and 2 f simultaneously by α while (15) corresponds to displacing , 1 2 16) and the definitions (17) we deduce that rays too strong in either displacement or angle are superexponentially suppressed. As mentioned earlier, the case l p 0 = = (Gaussian beam) is the only Wolf function that is completely positive. As soon as u or v are non-zero, there will be a sign change when μ or ν cross each root of the corresponding Laguerre polynomial. An example of the Wolf function of LG 11 is shown in figure 1. A typical cross section in the ray decomposition provided by interpreting the Wolf function as a generalized phase-space intensity is illustrated in figure 2. The graphic representation displays several features of the Wolf function when interpreted as a phase-space intensity density. As expected, a majority of the total positive intensity lies in the half-plane corresponding to the sign of the λ coordinate matching that of l (rays co-rotating with the beam), but the separation is not exclusive. Also, one can clearly see that the function is regular and bounded over its entire domain. This is contrary to the behaviour of Poynting vector near the origin, leading to a famous problem of anomalous momentum kick discrepancy [26]. This problem has been solved by an exact treatment of a test particle's uncertainty relation but the generalized ray picture shows that its very occurence can be attributed to the limitation to Euclidean space. The problem is not observed when full phase-space is considered.  1, | | + + = + + matching the respective eigenvalues in the wavefunction picture. This example, albeit being a simple demonstration of the more general theorem mentioned above, underlines the plausibility but also practicality of treating the Wolf function (16) as a phase-space quasi-density of intensity.
The Supplementary Material to this work presents several cases of LG beams as computer animations of figure 2 in evolution along the z-axis.

Bessel beams
Bessel beams are rotationally symmetric, non-diffusive solutions to the Helmholtz equation (1) is also a solution to the paraxial wave equation

Like the LG beams, they satisfy
The eigenvalue property of (12) is replaced by which grants the stability of the solution with respect to (3). Bessel beams are non-normalizable as a beam with a wavefront of the form (19) with any non-zero amplitude would carry an infinite amount of energy. The Wolf function can still be introduced but some of its properties will be modified in connection to its own unnormalizability. In order to compute it, we rewrite u r, , 0 ( ) f as a generalized superposition of two-dimensional plane waves: J r e 1 2 e e d .  k + are reachable. For all other phase-space points the value is zero. The assignment of α and β to a given k x and k y within the supporting domain is not unique but any choice leads to the same value of the righthand side of (22).
Bessel beams do not feature a naturally defined beam width so they do not have a particularly simple form in the coordinates , , , , ) m n f f defined above, for any w 0 . They can be thought of as having an infinite value thereof, and thus an infinite Rayleigh range, in agreement with their zero diffraction. This is further supported by the fact that a Bessel beam can be obtained as a limit of LG beams as p  ¥ and w p 2 0 k = [34, equation (22.15.2)], as also noted in [24]; also the eigenvalue equation which is the symbol of the differential operator on the left-hand side of (21), conveniently denoted Ĥ 4 . The value η, interpreted as a coordinate in phase-space, can be thought of as twice the limit of the above introduced ò in the infinite beam waist limit, rescaled by an appropriate power of w 0 : The coordinate 2 ( ) l m n = is still natural to the system due to (20). We perform an analogous transform on the angular coordinates, replacing them by their average for all w 0 and thus also x y k k , , , , x y The Cartesian phase-space coordinates can be expressed using the new four-tuple as (for a graph, see figure 3) and thus exhibits a symmetry with respect to translations in both θ and χ. The former is the same rotational symmetry as in (14). The latter generates the transform group x y x x y y which can be seen as a limit of (15) as w 0  ¥ and 0 b  simultaneously, maintaining w .  (22) and (25) were chosen so as to satisfy this condition.
One might naturally expect that the Wolf function (25) could be obtained from (16) using the limit procedure outlined above. We note that this is possible but is analytically challenging as the convergence takes   3). The horizontal axis was rescaled accordingly to the limit procedure (23). As p further rises, the apparent horizontal ripples get denser but do not decrease in amplitude.
place in a weak sense only; the latter converges to the former as a distribution but neither pointwise nor in any Lnorm. The character of the convergence is illustrated in figure 4. It follows from the symmetry with respect to (26) that the ray decomposition of a Bessel beam, as described by its Wolf function interpreted as an intensity distribution, will comprise parallel rays of constant intensity arranged in planes defined by the wave vector and the propagation direction. In a cross section, we would observe this as lines of constant flow along their transverse k-vector. The whole space is thus covered by positive and negative rays, non-decreasing in intensity even at high radii but rather asymptotically cancelling each other perfectly. In areas near the beam axis positive rays are more abundant due to the detailed properties of (22), resulting in the characteristic concentric ring intensity profile. All of these aspects can be fully appreciated only at extremely high sample rates, or in motion, as provided by two examples in the Supplementary Material.

Mixtures of Gaussian beams
The preceding examples dealt with individual, perfectly coherent or 'pure' wavefronts. The conditions on purity along with those on rotational symmetry and symmetry with respect to fractional Fourier transform, or to freespace propagation, translate to eigenfunction equations of (11) and (12) or (21) and leave LG or Bessel beams as the only solutions, respectively.
A Wolf function, however, can also be calculated for a statistical mixture of different beams, opening new possibilities for constructing mixed beams exhibiting rotational symmetry and propagational stability. In a sense, the full characterization is trivial: rotational symmetry of a Wolf function alone implies that the field is represented by a single wave function u x y , ( )which is an eigenfunction of (11) or a statistical mixture thereof; similarly, the symmetry with respect to (15) constrains the individual wave functions composing the field to eigenfunctions of (12) (not necessarily corresponding to the same eigenvalues). If these two conditions are put together, the field is restricted to be a statistical mixture of LG beams with various l and p but of the same beam waist w 0 . However, we might still seek for Wolf functions which are special in certain ways.
An interesting class of paraxial field solutions are those which can be written as a statistical mixture of Gaussian beams only. As the Wolf function is linear in the second-order correlation function (represented by the density operator), statistically mixing different components results in an affine combination of their respective Wolf functions, in turn the Wolf functions of any mixture of Gaussian beams is positive in all points. The axial Gaussian beam shows both the required symmetries, as shown above by being a special case of LG , 00 and we can use it to construct more complicated fields while maintaining both symmetries by forming equal-weight mixtures of its displacements to phase-space points along their phase-space orbits.
The principal case is obtained by displacing the Gaussian Wolf function W x y k k , , , 1 e x y x y w k k w 0 2 where (see (18)  by either a direct computation using (31) or by averaging (28) and (29) with the same weight distribution. The mixtures of Gaussian beams described in this section are capable of exerting a torque in a manner similar to LG beams without the presence of helical wave fronts or optical vortices. However, proposed information encoding capabilities of LG beams [10,12] are hindered by the mutual non-orthogonality of the solutions for different l p , . The Supplementary Material visualizes the propagation of a generic case of a mixture of Gaussian beams corresponding to 2 m = and 0 n = (l = 1, p = 0) for a comparison with the corresponding LG mode.

Conclusions
LG and Bessel beams are two important classes of scalar optical beams carrying OAM. They are characterized by rotational symmetry along with preservation of spatial structure under propagation (LG beams) or a full symmetry under free-space propagation (Bessel beams). We presented their respective Wolf functions and how these two symmetries naturally manifest themselves in phase-space. We discussed the role of negativity in the induced generalized ray decomposition and compared our approach to introducing rays as vector lines of the Poynting field.
Furthermore we introduced a simple class of statistical mixtures of Gaussian beams which exhibit similar properties as LG beams: they are rotationally symmetric and preserve the shape of their intensity profile as they propagate along the beam axis. They can be assigned an effective value of non-negative l and p parameters, which are not restricted to integers. More solutions can be obtained by further mixing of these fields.
We know that negativity of the Wigner function is an inevitable effect if non-Gaussian wave functions in quantum mechanics are studied [31], and correspondingly one needs to introduce negative ray pencils in order to describe LG and Bessel beams in terms of generalized ray optics. However our last studied case shows that the negativity of the Wolf function is not a prerequisite for a non-zero l (or p). It only becomes inevitable, along with quantization of the beam parameters, as a consequence of the restriction on 'pure' wavefronts.