Gaussian mode families from systems of rays

Gaussian mode families, including Laguerre–Gaussian, Hermite–Gaussian and generalized Hermite–Laguerre–Gaussian beams, can be described via a geometric optics construction. Ray families crossing the focal plane are represented as one-parameter families of ellipses, parametrized by curves on an analog Poincaré sphere for rays. We derive the optical path length weighting the rays, and find it to be related to the Pancharatnam–Berry connection on the Poincaré sphere. Dressing the rays with Gaussian beams, the approximation returns the Laguerre-, Hermite- and generalized Hermite–Laguerre–Gaussian beams exactly. The approach strengthens the connection between structured light and Hamiltonian optics, opening the possibility to new structured Gaussian beams.


Introduction
Gaussian beams are ubiquitous in contemporary optics as the simplest model of paraxially propagating, monochromatic light beams, as are the mode sets based on them, the Hermite-Gaussian (HG) and Laguerre-Gaussian (LG) modes [1][2][3]. The fundamental Gaussian and HG modes are familiar from the laser laboratory as modes of laser cavities, and LG modes epitomize structured light beams, carrying quantized orbital angular momentum and optical vortices on their axis. As such, HG and LG modes are the most familiar and studied examples of structured light beams.
This familiarity desensitizes us to a remarkable property all of these beams share-up to a scaling, their intensity profile is the same in every plane, even to the far field. This is readily accounted for mathematically by the fact that Gaussian beam families are eigenfunctions of (isotropic fractional) Fourier transformations, or in the analogy between paraxial optics and 2+1-dimensional quantum mechanics, that two-dimensional harmonic oscillator eigenfunctions spread, under evolution by the free Schrödinger equation, in a formpreserving way [4][5][6]. The two-dimensional quantum harmonic oscillator itself is analogous, via the 'Schwinger oscillator model' [7,8], to the quantum theory of angular momentum: Gaussian beams become analogous to spin directions in an abstract space (LG vertical, HG horizontal) with mode order proportional to total spin [9]. The simplest example of mode order N=1 therefore resembles the Poincaré sphere for polarization [10], giving rise to the celebrated Poincaré sphere for modes [11][12][13][14][15][16]. Modes corresponding to an arbitrary abstract spin direction are the generalised Hermite-Laguerre-Gaussian (GG) beams [17,18].
In [19] we proposed a new and broad definition of structured Gaussian beams, using the mathematics of semiclassical ray optics, based on the self-similarity of propagating Gaussian beams. This involved identifying each Gaussian beam with a two-parameter family of light rays, potentially overlapping, each of which is weighted by a complex amplitude: the sum of ray amplitudes at each point approximates the complex amplitude of the propagating wave field. This approximation can be improved by replacing each linear ray with a Gaussian beam whose axis corresponds to the path of the ray, for example according to the 'SAFE' method [20][21][22][23]. The main emphasis of [19] was the study of the geometric connection between the representation of the structured Gaussian beams over a Poincaré sphere and the beam's caustic distributions. However, several aspects of the theory presented in [19] were rather complicated, possibly detracting from the simplicity of the basic construction, which is firmly rooted in the Hamiltonian approach to optics.
Our aims here are to simplify the construction of the wave approximations in [19], in a way that shows a direct connection to geometric phases, and to demonstrate two semiclassical ways of approximating complex amplitudes with rays, based on two methods: WKB (where each ray contributes to the amplitude of the points it passes through) and SAFE (where rays are dressed with Gaussian widths). We show in particular that the latter gives the HG, LG and GG beams exactly. Such an approach provides a deeper understanding of the complex amplitude structures of Gaussian modes as interference patterns, exemplifying their geometric structures, especially Cartesian and polar coordinates for HG and LG beams respectively, and their caustics. This complements the insights of [19] where properties such as self-healing and geometric phase of parametrized families of modes were considered. We aim to provide a self-contained description that brings out the new features emphasized here.
In the next section, we describe the approach in [19], and give a definition of the GG beams, which are less familiar than their HG and LG counterparts. In section 3 we discuss the ray construction for Gaussian beams in detail, derived from first principles from the basic optics of Gaussian beams. In section 4 we determine the correct optical path length phase factors weighting each of the rays, finding the optical path length accumulated by ellipses on the Poincaré sphere path to be analogous to the Pancharatnam-Berry phase accumulated by optical polarization. The ray approximation is refined with a wave approximation in section 6, which is used in section 7 to show that our ray-based wave approximation is exact for LG, HG and GG beams.

Summary of ray approach to Poincaré sphere for Gaussian beams
The ray-based representation of Gaussian beams directly invokes a construction analogous to the Poincaré sphere for polarization [24], whose Cartesian coordinates are the Stokes parameters. More generally, the Poincaré sphere of ellipses is the space of closed elliptical orbits of a classical 2D isotropic harmonic oscillator (2DHO). As such it is effectively the reduced Hamiltonian phase space of the full 4D phase space of the 2D oscillator, via the Hopf fibration [25,26].
In the ray approximation of a Gaussian beam in [19], each member of a two-parameter family of rays, with parameters labeled η, τ, crosses the focal plane at a point: these points are arranged to be on ellipses ('elliptic orbits') parametrized by τ, themselves represented by points on a curve on a Poincaré sphere, parametrized by η. The total ray family is then represented by this curve on the Poincaré sphere: in the focal plane the ray family corresponds to a smooth family of ellipses. The particular curves for LG, HG and GG beams are circles whose center direction is the point on the mode sphere parametrizing that beam. A representation for some examples is shown in figure 1.
The most striking feature of the ellipse families is their envelope curves, which are caustics. The LG and HG modes have caustics which are concentric circles and rectangles respectively, reflecting the fact that these modes can be written in separated polar and Cartesian coordinates [27] (equations (7.1) and (7.4) later). By contrast, the caustics of GG beams are rather more complicated, and in general GG beams cannot be written simply in terms of any coordinate system. By the analogous quantum angular momentum algebra, a GG beam corresponding to the direction θ, j (representing polar and azimuthal angles on the Poincaré/mode sphere) can be expressed as a Figure 1. Examples of Gaussian beams and their corresponding elliptic orbit families on the Poincaré/mode sphere. The pink sphere represents the Poincaré sphere of ellipses/modes, on which are drawn four circles, corresponding to Gaussian modes and their elliptic ray family representation: LG 2,2 (blue), HG 5,5 (green), GG 2,2 (orange), GG 5,5 (yellow), the latter two corresponding to the point θ=π/4, j=0 on the sphere. As in the main text, modes here are labeled by N and μ. quantum angular momentum-like sum of LG beams [4,9,17] In (2.1), N represents the mode order (analogous to a quantum spin N/2) which is common to the GG beam and LG beams in the sum, and μ is the GG order, μ=−N,−N+2, K, N. The order index μ labels all GG modes, and for LG beams (θ=0) corresponds to the angular momentum label ℓ; the radial order being . For HG beams with x and y indices m, n, μ=m−n [4,9,17]. Throughout this paper we will label all Gaussian modes in terms of N and μ, see r LG N,m¢ ( ) in (2.1). Our approach here differs from that in [19] in several ways. The key aspect to being able to use rays to approximate the field is to calculate the correct reduced optical path length (ROPL) S(η, τ), which provides a phase factor e S i , h t ( ) weighting each ray labeled by η, τ 5 . We find this here by considering the phase relationship between different ellipses on the Poincaré sphere analog, deriving a phase connection analogous to the Pancharatnam-Berry phase [28]. This is significantly simpler than the calculation in [19], which used the caustics for reference. Note that, when using the WKB method, a detailed analysis of the Maslov indices incurred as caustics are crossed is required as an extra step, while the SAFE method integrates such phases naturally. One of our main results is equation (6.6), which shows that the Gaussian-dressed ray-based wave approximation can be considered as a continuous integral over GG modes with μ=N. We prove that the approximation (6.6) is in fact exact for HG, LG and GG modes, following a direct argument based on the geometry of the Poincaré sphere and the expressions (2.1), (2.2). Some notation is changed from that in [19].

Determining ray families for Gaussian beam families
In free space, a light ray is a straight line; parametrized by propagation distance z, its path can be written where z is the optical axis unit vector, q is the transverse position of the ray in the focal plane, and p the transverse component of the ray direction, i.e.the ray momentum, such that p arctan| | is the angle the propagating ray makes with z . In the paraxial regime, as we assume here, p 1  | | . In the sense of Hamiltonian optics [29], the rays' transverse position q and direction/momentum p uniquely determine the path of the ray, and are canonically conjugate.
Gaussian beam families are distinguished in the fact that they have a well-defined beam quality factor 2  [30], which can be expressed, for a paraxial light beam of width w 0 , wavenumber k, by For a beam, 2  can be considered as an operator (with the momentum p k i , , for which a fundamental Gaussian in the waist plane, e r w 2 0 2 -, is an eigenfunction with eigenvalue 1, as are its HG and LG counterparts, with eigenvalues N+1, where N is the mode order of the beam [30]. The beam quality operator is equivalent, in the quantum mechanical analogy, with the 2D quantum harmonic oscillator Hamiltonian.
All rays representing a Gaussian beam should therefore share the same value for 2  with q 2 as r 2 , which should be equal to the desired value N+1. It is therefore convenient to define the dimensionless ray position Q and momentum P such that so that for all rays representing the Gaussian beam, Q 2 +P 2 =1. This defines the three-parameter set of all possible rays with constant quality factor: rays correspond to points on a unit three-sphere in four-dimensional, dimensionless phase space. Any real light beam is represented by a two-parameter family of rays, parametrized by η and τ, which is a smooth two-surface in the three-sphere of rays with fixed 2  .
The main construction of [19] was to use key properties of Gaussian beams (including high-order HG, LG and GG families) to motivate the choice of ray family: beams propagate on a hyperboloid, and for each constant z, the intensity pattern is self-similar to that in the focal plane up to a scaling. This can be realised by ensuring that each ray belongs to a one-parameter ray family which sweeps out a ruled hyperboloid on propagation. Provided the phase weightings of each ray around the hyperboloid are arranged consistently, the contribution from the ray family at a given point in the transverse beam profile will be the same for each value of the propagation parameter z. For a given ray with dimensionless position and momentum Q P , , such a hyperboloid is swept out, by varying τ, by all of the rays with transverse waist plane position , as shown in figure 2: the hyperboloid has an elliptical cross-section. Each ray on the hyperboloid has the same skewness z q p L k = , equivalent to the angular momentum of the linear ray about the beam axis, and τ=τ 0 gives the original ray Q P , . Remarkably, this is precisely analogous to the Jones vector description of a polarization ellipse: each hyperboloid of rays in 3D has an elliptic 'orbit' in the transverse plane, represented by a complex, normalized Jones vector Q P i + .
The definition above suggests the convention that the ray family parameter τ=0 corresponds to the dimensionless Jones vector v A B i = + , where A is along the ellipse major axis, B along a minor axis, and ; the phase of v has been chosen such that v v · is real and nonnegative; in fact v v cos J = · . Many of the following calculations are facilitated by using the latitude angle ϑ on the Poincaré sphere, rather than the more familiar polar angle θ, The azimuth on the Poincaré sphere is denoted by j to distinguish it from the azimuthal angle f in configuration space, as used already in (2.1).
The dimensionless position Q t ( ) and momentum P t ( ) corresponding to any ray on the orbit can be found from These conventions closely follow the use of the Poincaré sphere for polarized light.
Each hyperboloid of rays is represented by an ellipse given by the Jones vector and represented by a point on a Poincaré-like sphere. The other ray family parameter η determines a curve ϑ(η), j(η) on the sphere, which we call the Poincaré curve (or Poincaré path). Following [19], we assume that this Poincaré curve closes (i.e. its endpoint is the same as its initial point), so we can set 0 2   h p: the two-parameter ray family is topologically a torus in the three-sphere; the cycle in τ wraps around each ellipse (i.e. a Hopf circle in the three-sphere), the other traces out a loop on the Poincaré sphere. The fact that all rays were chosen to have the same value of 2  ensures that each elliptic ray orbit has the same size, equivalent to the normalization of the Jones vectors on the path. Elliptic orbits of rays propagating on a hyperboloid. The rays are colored by hues according to phase, which increases along the ray with optical path length (distance). Cross-sections of the hyperboloid are all given by the same ellipse, and at different propagation distances, the phases around the ellipses correspond.
A major result of [19] was that only the Poincaré curves are allowed which enclose a solid angle Ω which is an odd multiple of 2π/(N+1), i.e. Therefore we can expect curves with different areas correspond to differently labeled modes within a Gaussian beam family with fixed N. We will show here that this quantization follows directly from a simpler argument involving the analog of the Pancharatnam phase for polarization on the Poincaré sphere analog for ray orbits. The resulting family of ellipses is bounded by an envelope of caustics. As discussed in the Introduction, LG, HG and more generally GG modes are represented by elliptic orbits following circular Poincaré paths, whose centers corresponding to that mode on the mode sphere, and with radii determined to enclose one of the allowed areas (3.6), as illustrated in figure 1. These caustics can be considered as the edge of the projection of the phase space torus (parametrized by 0 , 2   h t p) into real space, which is two concentric circles for LG modes and a rectangle for HG modes, and, for GG modes, a more complicated structure interpolating these. As the parameters η, τ vary across a caustic, in addition to the ROPL, the phase acquires a discrete jump: a 'Maslov index' (see e.g. [31,32]) which we will consider in detail here for HG, LG and GG modes. However, we will not develop the investigation in [19] of the intricate relationship between the Poincaré sphere curve and the caustic curves.

Determining the ROPL for each ray
We now turn to finding the ROPL S(η, τ) corresponding to each ray with position q , h t ( ) crossing the focal plane z=0. Associating complex amplitudes with rays is standard semiclassical Hamiltonian optics, although it is somewhat beyond elementary geometric optics. In particular, rays correspond to saddle points in a diffraction integral, and via a straightforward WKB approximation, each ray acquires a weighting of the complex amplitude is the Jacobian matrix transforming between configuration space Q and the space of rays parametrized by η, τ. We will primarily focus on the phase here, given by the ROPL S(η, τ).
Between neighboring rays, the infinitesimal ROPL S d is given by the standard ray connection p q k d · [29]. From the relationships (3.4) and (3.5), we therefore have where ϑ and j are functions of η, and d d j Since the ray family is a normal congruence, the difference in ray amplitude (4.1) between any pair of rays does not change under smooth deformation of the path with fixed endpoints. From the parametrization discussed above, the major axis of the ellipse labeled 0 h = also has 0; t = we denote this ray by O. We denote the path difference to any other ray η, τ by . This determines the phase weighting of each ray in the family (relative to the (arbitrary) phase of O).
First, we consider the relative phase around the ellipse from When τ=2π, the accumulated ROPL is 2π (N+1) as it continuously returns to the ray labeled by τ=0, which agrees as it must. In fact, this is a semiclassical quantisation condition that forces the mode order of a Gaussian beam to be an integer, so that the phases weighting the rays around each elliptical ray orbit are singlevalued: for every ray orbit in a Gaussian beam where N 1 2  = + the phase cycles around each orbit N+1 times. Equation (4.3) shows that the phase weighting around the orbit is not just proportional to the phases around an elliptic orbit (i.e. to τ): there is an extra oscillatory decoration vanishing at the major and minor axes, depending on the latitude ϑ (i.e. the ellipticity of the ellipse). The Poincaré sphere in figure 3 shows ellipses with N=2 according to this phase weighting.
We now consider the dependence of the path length on η, i.e.the ROPL between neighboring ellipses on the Poincaré sphere. It might be natural to expect this to be given by a Pancharatnam-Berry geometric phase connection familiar from polarization optics [28], which requires two neighboring Jones vectors v , Indeed, the j¢-dependent part of the connection (4.2) is the geometric phase connection dG times N+1; we recognise the extra J¢-dependent part as the ϑ derivative of the decorated phase (4.3) at τ (which is missing from the standard Pancharatnam-Berry phase). Taking as reference the points with 0 t = as h increases from O to η, the geometric phase is and the associated ROPL is therefore By considering the path in the ray family as starting at O, increasing h to η gives a path length S 2 (η), and then increasing t to τ around the ellipse labeled by η gives ROPL S 1,η (τ), therefore giving the total As continuity of the ray family in (4.3) gave us our quantization N+1 around the ellipse, we now derive the area quantization (3.6). It is standard in geometric phase theory that the area Ω enclosed by a Poincaré curve is associated with the total geometric phase of the circuit, 2 s i n d [28]. When the Poincaré curve does not enclose the north pole (in a right-handed sense), the total geometric phase On the other hand, on a path enclosing the north pole (which involves a full circuit of j, 0 2   j p), Γ is equal to half the area enclosed between the curve and the equator of the Poincaré sphere, 2 1 2 p G = -W ( ) (discussed in sections 3 and 4 of the Supplementary Material of [19]). Note further in the latter kind of circuit, the ellipse executes half a turn, so A

WKB approximation, Maslov phases and solid angle quantization
In order to give a WKB approximation of the complex amplitude, as the sum over the rays through each point correctly weighted by their complex amplitude contributions, it is necessary to include contributions due to the Maslov phases: at each position a path in η, τ space crosses a caustic, the phase jumps by ±π/2 [31]. This can be understood as originating from (4.1): at a caustic, the absolute value of α goes to infinity as the inverse Jacobian matrix is singular. As η or τ varies across the caustic the sign of the Jacobian therefore changes, giving a net phase of ±π/2 due to the square root, with the sign as discussed below. These Maslov jumps therefore need to be added to the ROPL for the correct continuity between ray weights. The Maslov indices acquired on a general path in η, τ space therefore depend on the particular caustics associated with the Poincaré curve. Firstly, over a cycle in τ with η constant, we follow the path around a particular elliptic orbit in real space, which touches a caustic at four points, leading to a total Maslov contribution of −2π subtracting from the 2π(N+1) from the ROPL. Now consider the Maslov phases over a cycle in η, with τ=0 fixed. This path follows the ellipse major axis A h The number of caustics touched by this path can most easily be counted by looking at the Poincaré path, or more conveniently, at the representation of this path onto the plane with Cartesian coordinates , j J ( ). Evidently there are two types of caustic contact points for A h ( ). The first corresponds to segments of the Poincaré path that are parallel to the equator, that is, to values of η for which , since these correspond to infinitesimal bundles of ellipses that are simply rotated versions of each other. The sign of their Maslov contribution depends on the local curvature, the phase being given by sign 2 J j p -¢¢ ¢ ( ) . Therefore, for circular GG paths that do not enclose one of the poles (so that the major axis forms a closed loop) the total Maslov phase is −π, while for circular GG paths enclosing a pole (where A )) the two contributions cancel and the phase is 0. The second type of caustic contact point corresponds to crossings of the equatorial line, namely to values of η for which ϑ(η)=0, so the orbits are linear. These contact points always appear in pairs in a way that their effect to the total Maslov phase over the Poincaré path cancels. Examples of these contact points, the neightbourhoods of ellipses and the local caustic, is shown in figure 4.
This argument also can be used to show that (3.6) holds, i.e.Ω must be an odd integer multiple of 2π/(N+1 where the initial contribution of −Nπ accounts for the phase weighting differences at A and A -. Examples of ray families with correct ROPL and Maslov weightings can be seen in figure 5. These correspond to three cases shown in figure 1; HG beams are not shown as their ray families are quite degenerate, with each ellipse occurring twice but with different handedness and weighting. Clearly, the phases are discontinuous at the caustics due to the Maslov indices, and otherwise the ROPL phases smoothly join around the loops, such that the net phase change around each ellipse is 2π N. There are four types: • minima in ϑ, an exemplar of which is shown in (i); • maxima in ϑ when the loop does not include a pole, with exemplar in (ii); + maximum in ϑ for a curve enclosing a pole, as in (iii); × a curve crossing the equator, as in (iv). In (i)-(iv), several elliptic orbits are shown in the neighborhood of the crossing; the ellipse whose major axis meets the caustic is plotted thicker, with the major axis as an arrow. A section of the caustic is also shown in each case, with an arrow indicating the direction the caustic is swept out as η increases (anticlockwise for (i), (iii), (iv), clockwise for (ii)). All of the ellipses are right-handed in (i)-(iii), with orientation indicated by hollow arrows in (iv), where the handedness changes on passing through the linear orbit on the equator. These considerations can be used to give a WKB approximation for the complex amplitude at each point Q in the transverse plane, where at each point, the amplitude is given by the sum of complex amplitudes (4.1) corresponding to each ray through that point, including the Maslov phases. This crude approximation diverges on the caustics, but otherwise gives a good approximation of the mode behavior; figure 6 gives some examples of this approximation for the LG modes with N=5. These discontinuities can be improved by more refined semiclassical arguments, replacing the caustics with the appropriate diffraction catastrophe function, i.e.the Airy function (along folds) and the Pearcey function (at cusp points) [33]. In general the method is hampered by the fact that it is hard to find Q Q , h t ( ) ( ), that is, the parameters of the rays going through a particular point, and further details of these calculations are omitted. However, it can be shown that the HG case reduces to the product of WKB approximations of the one-dimensional quantum harmonic oscillator in each of x and y.

Weighting ray families with Gaussian beams
The WKB representation of a complex wave field by a sum of complex amplitude-weighted rays at each point is rather crude, particularly leading to an approximate amplitude which diverges at the caustics (and is discontinuous in phase due to the Maslov index). Another approach, which is directly wave-theoretic, is to treat each ray as a fundamental Gaussian beam whose center is at q and propagation is p. Such a beam, in the waist plane, has amplitude r q p r q r q p g w w k ; , 1 exp i , 6.1 where w 0 is the same width as we considered in the expression for the beam quality factor 2  (3.2). We write the fundamental Gaussian r g ; 0, 0 ( ) as r g 0 ( ). The beam (6.1) has an expectation value of 2  of N+1 when q and p satisfy (3.3), as we required for our rays representing high-order Gaussian beams: our field approximation will dress the rays as Gaussians of width w 0 .
There are many schemes to assign complex amplitudes to Gaussian-dressed systems of rays. Our choice of the SAFE method [20] is motivated by the fact that the amplitudes assigned to the Gaussian-dressed rays are optimized for the Gaussian widths. In this case [19], the Gaussian ray labeled by η, τ is weighted by the complex amplitude , where h(η) is some positive real function of η determining weight along the Poincaré curve (usually set to be unity), and S(τ, η) is the total path length (4.7), (4.5). The term is the complex Jacobian determinant of the components of the Jones vector with respect to τ, η, This term automatically accounts for the Maslov index smoothly: on a cycle 0 2   t p, it accumulates −2π, equivalent to the four Maslov indices around the ellipse. Similarly, for a cycle 0 2   h pthe total phase equals half the angle accumulated by the tangent vector to the curve (j(η), ϑ(η)), which also agrees with the prescription described in the previous section.
All together, this gives the field for our Gaussian-dressed ray approximation of a high-order Gaussian beam, where  is an undefined normalizing factor. Bringing together the integral of the terms in τ together with the Gaussian factor, we define the function U N , where the last equality follows from a complex residue argument ( [19], supplementary information (S26)), and H N is a Hermite polynomial. r U , ; N j J ( ) is a normalized, structured Gaussian beam, depending on N and a Jones vector v (i.e. a point on the Poincaré sphere). These beams have recently been investigated by Kotlyar et al [34] who called them 'vortex Hermite Gaussian beams', describing their relationship with various kinds of structured Gaussian beam and mentioned they are special cases of GG beams, but without further description or proof. They were also considered in the context of quantum mechanical harmonic oscillator eigenstates by Pollet et al [35].
Remarkably, this U N mode is exactly a GG beam with μ=N, r r U , ; GG , , , 6.5 as shown in the appendix. We can therefore understand GG beams with μ=N as continuous sums of Gaussian beams whose centers follow an elliptic orbit, phased appropriately according to (4.3). The representation (6.4) is a surprisingly simple closed-form expression for the GG N,N beam, hitherto only written explicitly in forms equivalent to the angular momentum sum (2.1).
), which define the caustics at which the ray amplitudes diverge. The f dependence of the WKB approximation is e imf , just as for exact LG modes.
We can therefore rewrite the wave approximation as This equation is one of the main results of this paper. It states that the Gaussian-dressed ray approximation to our Gaussian beam is equivalent to the continuous superposition of GG N,N beams, whose position on the Poincaré sphere follows the Poincaré curve, and that are weighted by a complex amplitude that includes the geometric phase accumulated by the curve. Thus the ray approximation has been completely incorporated into the fact that each Gauss-weighted elliptic orbit is exactly represented by the corresponding μ=N GG mode. In [19], U N was not connected to a GG beam, nor was the corresponding path length identified with a geometric phase. We also note that the integration of the square root Jacobian factor has been incorporated into the definition of U N and (6.6) smoothly: the sign change associated with Maslov indices has automatically been incorporated. Notably, this means that the phase change around the axis of the GG N,N mode changes N times, not N+1 as in the corresponding elliptic ray orbit (due to a factor of e it from the square root Jacobian (6.2)).
Two GG states arising from this integral over U GG N NN , = modes are shown in figure 7. Equation (6.6) is a convenient representation of this superposition, but it is not unique. If we were to define the Jones vectors on the Poincaré path to have continuous phases by the Pancharatnam connection, say Replacing v with V does not affect the argument of the Hermite polynomial since equal powers of V occur in the numerator and denominator; the multiplying factor of v v N 2 ( · ) accounts for N times the geometric phase factor e i h G( ) . After removing a constant and the τ-dependent factor from the Jacobian , one is left with v v ¶ h (the cross product understood here to give the scalar component perpendicular to the v plane). The replacement of v with V takes account of the final geometric phase factor, meaning that the approximation (6.6) can be rewritten purely and succinctly in terms of the path-dependent Jones vectors Here the geometric phase weighting has been completely absorbed into the definition of V . In practice, in numerical implementations of this integral or (6.6), care needs to be taken to ensure the integration path is free of branch cuts and the phase of the integrand varies smoothly with η.

The wave approximation for LG, HG and GG modes
We are now in a position to calculate the wave field approximation (6.6) around the circular paths representing LG, HG, and GG modes. We begin with LG modes, which are the most straightforward, and will pursue a similar strategy for HG and GG modes below. Normalized LG modes are given in the focal plane by where w 0 is the waist width in the focal plane, r and f denote radial and azimuthal plane polar coordinates respectively, and L p ℓ | | is an associated Laguerre polynomial [36]. The angular momentum label ℓ here corresponds to the quantum number μ, ℓ=−N,−N+2, K, N, and the radial number LG modes are eigenfunctions of the beam quality factor operator 2  with eigenvalue N+1, analogous to the angular momentum eigenstates of the 2D quantum harmonic oscillator.
From the prior discussion, we expect to be able to get an approximation of (7.1) using (6.6), on a circular path on the Poincaré sphere around the vertical axis, on which η=j at height ℓ/(N+1), i.e.constant N arcsin 1 . J = + ℓ ( ( )) In the expression (6.6), h(η)=1 and the Jacobian square root factor is constant: the integral only involves U N times the geometric phase factor. The associated ROPL is the last equality following from the fact that (μ−ℓ) is necessarily an even integer, so the integral over η gives a Kronecker symbol δ μℓ . In this calculation, the undefined  absorbs numerical constants in each step. We have therefore shown that for constant latitude ϑ (i.e. constant θ, Z), ψ is exactly a LG mode. We will approach the calculation of ψ for HG and GG modes in a similar way, replacing U N with the appropriate GG mode. In (7.3), the calculation was made straightforward as LG modes appear in the GG sum (2.1), so the integral around the Poincaré path reduced to a Kronecker symbol on the correct LG mode. Since GG modes are written in terms of angular momentum superpositions of angular momentum basis states, (2.1) is not a unique representation. In particular, we will proceed by finding a representation of GG modes in terms of HG modes. The normalized HG modes are defined as An interpretation of the GG mode (2.1) is by analogy with quantum angular momentum addition, since Gaussian modes, as analogous to quantum isotropic harmonic oscillator states, are realisations of Schwinger's oscillator model of quantum angular momentum [9]. The LG mode LG N,μ is analogous to the spin state with total spin j=N/2 and z-component of spin m=μ/2, where the i factor ensures the phase is correct, according to the Condon-Shortley phase convention [4,8]. The coefficients are then quantum rotation matrix elements for Euler angles χ, θ, f, which in general are given by where R n b ( ) represents a rotation about the n-axis by β. As the rotations in the present context are in the abstract space of the Poincaré sphere, Cartesian axes are referred to as X, Y, Z. The LG mode corresponds to the state in the Z-direction, HG to the X-direction, and (7.5) represents the rotation by χ about the Z-axis, then through θ about the Y axis, then by j about the Z axis. The second two rotate to the point labeled by θ and j; including the extra χ rotation would give an extra factor e i 2 mc from rotating the original LG mode.
A feature of the quantum treatment of rotations is that any basis of spin states given by fixed j, m=−j, K, j in some direction can be used to represent spin states in other directions. Extending this to the sphere for modes, we should be able to express each GG mode as a superposition of HG modes with the same mode order, viz.
where , q j are spherical polar coordinates with the X-direction as poles corresponding to the usual θ, j with the z-direction as poles, so arcsin cos cos J j J = -( ) and arctan cos cos sin ) ensure the HG modes follow the Condon-Shortley convention. In order that this new GG mode set agrees with that previously defined, we require an extra phase factor related to the third Euler angle c in the (1, 0, 0)-based spherical coordinates, as we now discuss.
The third Euler angle gives a measure of rotation about its own axis with respect to a 'flag' (e.g. the direction of the real part of a circular polarization, such as the x-direction for the right handed circular polarization state x y i +  ) on the conventional Poincaré sphere for polarization. If χ=0, as for our previous GG state (2.1) then a flag initially in the X-direction points in the q  -direction. Therefore, the angle between such a flag for Euler angles rotating from the Z-direction, and the X-direction, is the angle between the meridians with respect to Z and X at that point on the sphere. Straightforward spherical trigonometry gives this angle as where, remarkably, c is the Euler angle defined in (7.7). Now, the Jacobian must also be calculated; the term inside the square root is cos i cos e i J j J J ¶ + ¶ = j j c -, where in the argument exactly the extra Euler angle of (7.7) occurs again, and the modulus is independent of the parameter j, so can be absorbed into  . Therefore the net j-dependent weighting of U N in the integral for ψ is N exp ) . This χ-dependent factor exactly cancels the prefactor of U N identified as the HG-summed GG mode in (7.6) with μ=N. The remaining sum resembles (7.3) exactly, with the integral now over j, and the Condon-Shortley-weighted LG modes replaced by the appropriate HG modes: the ψ expression is therefore exact for HG modes. This was previously shown [19] using separation into Cartesian coordinates; the spin formalism here provides an alternative route, and is directly analogous to the LG proof. Finally, we consider the GG modes. Following the same strategy, we now wish to represent a GG mode corresponding to a point on the Poincaré sphere as a superposition of GG modes based at some other fixed point given by polar angles θ 0 , j 0 in the usual spherical coordinates. We denote this rotated coordinate system in our abstract space X Y Z , , ¢ ¢ ¢, with Z cos sin , sin sin , cos This sum resembles the other two; the GG modes automatically satisfy the Condon-Shortley convention, so no extra prefactors are necessary for the modes. The overall phase of the GG mode here is the same as in (2.1), (7.6), with the inclusion of the Euler angle c¢. As before, c¢ is the angle on the Poincaré sphere between the meridian at the point , q f ¢ ¢ from z ¢  , and the meridian at , q f ¢ ¢ from z , and again by tedious but straightforward spherical trigonometry, it can be shown that arctan sin sin sin cos cos cos sin , which is consistent with (7.7) for which θ 0 =π/2. For applications only involving the intensity, of course c¢ can be ignored. The integral in (6.6) now involves a path of constant N arcsin parametrized by j¢, in the spherical coordinate system with the chosen fixed GG state at ϑ=ϑ 0 , j=j 0 as poles. Although the expressions are more complicated, it turns out that the argument follows exactly the same structure as the HG case just considered: on these paths, the geometric phase can be shown to be N 1 ) , with c¢ being the third Euler angle defined in (7.10). The square root Jacobian has constant modulus and argument 1 2 c -¢, and the prefactor of (7.9) cancels the remaining c¢-dependence; integrating exp imj -¢ ( )times the sum in (7.9) picks out the GG N,μ times the constant  as required.

Discussion
We have reformulated and deepened the construction of [19] by which families of rays are associated with higher-order Gaussian modes. We have seen how the optical path length associated with elliptical ray orbits on the analog ray Poincaré sphere is related to the Pancharatnam-Berry phase on the Poincaré sphere, and get continuous ray families when the Maslov indices are correctly incorporated. Using this in a wave approximation with Gaussian-dressed rays (6.6) removes the discrete Maslov indices; this representation is effectively a coherent-state representation of N m = GG beams integrated on the Poincaré path on the appropriate mode sphere, again including the geometric phase. Finally we demonstrated that the Gaussian wave approximation is exact for LG, HG and GG beams.
Our proof of the latter has its own algebraic elegance (e.g. the cancelling phases from the third Euler angles , c c ¢ ); however, in the present work we have made no use of the rich conformal structure described in [19] between the Poincaré curves (realised as caustics of families of circles once projected into the equatorial disk of the Poincaré sphere) and the caustics (envelopes) of the ellipses families in real space. A more detailed analysis of the caustics associated with GG modes, exploring the 2D-projected torus caustics [9], might reveal further properties of the morphology of GG beams, such as the characteristic 'vortex rows' close to the axis when θ is small [18].
A major feature of the present work is the demonstration that the U N kernel in the (6.6) integral is in fact a GG N,N beam. This makes GG N,N modes coherent states, both as integrals of a fundamental Gaussian on an ellipse (in the usual optical sense), and as equivalent to a spin coherent state on the two-sphere [37]. In fact, the form of U N in (6.4) has been considered previously [34] as 'vortex Hermite Gaussian beams', reflecting the fact that the beam is a Gaussian times a Hermite polynomial with a complex-optical vortex-argument, but the connection with GG beams was not made (we observe that vortex rows property [18] is also related to representing GG beams by Hermite polynomials with complex argument). U N beams have also been considered as special cases of the Lissajous coherent state-related beams [38]. A different approach to GG beams, based on generating functions, has also been investigated [39]. These suggest future directions in which to extend to extend the ray-based analysis of Gaussian beams.
Our focus on LG, HG and GG beams has been more restrictive of the formalism than is necessary; more richness will come from new families of Gaussian beams based on more complicated Poincaré curves than the circles considered here. These may be understood in the sense of classical and semiclassical Hamiltonian optics, equivalent here to the quantum and classical mechanical 2D harmonic oscillator, as follows. The system has three classical constants of the motion, representing the orthogonal axes of the ray Poincaré sphere, i.e.Stokes parameter analogs s 1 , s 2 , s 3 of the rays. The one corresponding to s 3 is the ray skewness σ, i.e.the ray's angular momentum around the beam axis. A Poincaré curve may therefore be considered as a contour of some function J(s 1 , s 2 , s 3 ) of the Stokes parameters: the circles corresponding to a LG beam when J=s 3 , HG when J=s 1 , and GG in direction j, ϑ when J is a linear combination s s s cos cos sin sin ) (every linear combination of Stokes parameters corresponds to some GG mode). In a Hamiltonian sense, the Poincaré sphere is the reduced phase space on which J is an action, constant for all the elliptic orbits on a path; the conjugate variable is an angle ϖ around the contour, which may be η (if not, the reparametrization is included by the function h ). The quantization rule (3.6) means there are always N+1 distinct beams with quantized J.
The wave approximation (6.6) always gives an exact Gaussian beam; it is an approximation in the sense that the resulting beam only approximates the appropriate eigenfunction of the quantum operator Jˆcorresponding to the classical constant of the motion J(s 1 , s 2 , s 3 ) whose contour gives the given Poincaré path. For arbitrary curves, this will semiclassically improve for N 1  . It would therefore be interesting to investigate some simple examples of J beyond linear combinations of s 1 , s 2 , s 3 , and explore the connections between the Gaussian beams, the corresponding Poincaré curves, ellipse families and caustics, to investigate further relations beyond the GG beams considered here.
Since U N appearing in the semiclassical expression for ψ in (6.6) are GG modes, (6.6) and ( where here  is independent of m¢. Appropriate modifications can be made to represent the curve in terms of HG or GG basis sets.
We have found here a key role played by the geometric phase implicit in the underlying ray structure of LG, HG and GG beams. We note that in general, a beam evolving between different GG states acquires an overall geometric phase [13], which strikingly can be measured without interferometry using this ray picture [16]. It would be interesting to see how this geometric phase evolves under more general deformation of the underlying Poincaré path, which might reveal new aspects of geometric phases for more general Gaussian beam families.
The rich and elegant structure associated with the ray picture of Gaussian beams is related to the underlying symmetries of the classical and quantum 2D isotropic harmonic oscillator: its spin-like SU(2) structure is both superintegrable and originates from a hidden symmetry, allowing for the semiclassical quantization conditioneffectively the old quantum theory of Bohr and Sommerfeld-to be exact. We have seen how this quantization is purely geometric, in terms of explicit classical ray families and the corresponding geometric phase. One can speculate whether other families of structured light beam might have representations in terms of rays and symplectic structures with a similarly geometric quantization [40]. whose rhs is equivalent to (2.1). QED.