A ray-optical Poincar\'e sphere for structured Gaussian beams

A general family of structured Gaussian beams naturally emerges from a consideration of families of rays. These ray families, with the property that their transverse profile is invariant upon propagation (except for cycling of the rays and a global rescaling), have two parameters, the first giving a position on an ellipse naturally represented by a point on the Poincar\'e sphere (familiar from polarization optics), and the other determining the position of a curve traced out on this Poincar\'e sphere. This construction naturally accounts for the familiar families of Gaussian beams, including Hermite-Gauss, Laguerre-Gauss and Generalized Hermite-Laguerre-Gauss beams, but is far more general. The conformal mapping between a projection of the Poincar\'e sphere and the physical space of the transverse plane of a Gaussian beam naturally involves caustics. In addition to providing new insight into the physics of propagating Gaussian beams, the ray-based approach allows effective approximation of the propagating amplitude without explicit diffraction calculations.


I. INTRODUCTION
Structured Gaussian beams are amongst the most familiar examples of paraxially propagating light beams. These include the Hermite-Gaussian (HG) beams [1,2], with intensity patterns resembling Cartesian grids, and Laguerre-Gauss (LG) beams [3], whose intensities are concentric rings and whose phase carries orbital angular momentum (OAM). A remarkable feature of Gaussian beams is that their intensity profile does not change on propagation, apart from an overall scaling; even in the far field, HG and LG modes appear the same. More recently, other self-similar beams have been studied in detail, including Airy beams [4,5], which are self-similar on propagation up to a parabolic lateral displacement, and Bessel beams [6,7] and Mathieu beams [8] whose intensity profile does not change at all on propagation. Self similarity is more than a mathematical peculiarity, and is an important aspect of many applications of structured light. For example, given the constant width of their main intensity lobe, approximations to Airy and Bessel beams have been the basis of several imaging techniques, whether for illumination to increase axial resolution [9,10] or for 3D shaping of the point-spread function to increase depth of focus [11][12][13][14][15]. Given their characteristic intensity and phase profiles, structured Gaussian and other self-similar beams have been used extensively for particle manipulation [16][17][18][19][20], and may be eigenfunctions of natural optical operators, such as Bessel beams and LG beams of azimuthal order , which carry an OAM of per photon [3].
Here we describe an approach to structured Gaussian beams in terms of ray optics. Geometric optics is usu- * alonso@optics.rochester.edu ally applied in situations where a light field has welldefined extended wavefronts with slow amplitude variations. However, it turns out Gaussian beams and their generalizations are remarkably amenable to such an analysis. As fundamental Gaussian beams, as well as HG and LG beams, are modes of laser cavities with curved mirrors [1,2], their dynamics is well-approximated by a two-dimensional isotropic harmonic oscillator representing the transverse plane, with the mirror curvature acting as the harmonic potential. Classical orbits in the two-dimensional isotropic oscillator are of course ellipses, and can be represented by points on the Poincaré sphere, more familiar in representing the polarization of a harmonic electric field [21]. In our analysis, a Gaussian beam is represented by a two-parameter family of rays; the rays are divided into subfamilies describing ellipses which propagate in a self-similar way, and which are described by points on a Poincaré-like sphere for rays. The choice of the other parameter of the ray family then corresponds to determining a closed path of ray ellipses on this sphere, which is different for different types of beams. Consistency with ray optics forces quantization conditions on these parameters, both around the ellipse and on the Poincaré sphere path. These conditions give, for certain natural choices of path, the quantum numbers associated with HG and LG modes [22], although many new kinds of structured Gaussian beam are possible. An immediate generalization is to the generalized Hermite-Laguerre-Gaussian beams (HLG) [23][24][25][26][27][28][29], which interpolate between the HG and LG families on a generalized Poincaré sphere via an anamorphic fractional Fourier transformation, realized physically by transforming HG or LG beams through a beam shaping device consisting of suitably chosen pairs of cylindrical lenses [30][31][32].
We therefore are discussing objects very familiar in modern paraxial optics: mode families, optical opera-tors, geometric optics and Poincaré spheres, although combined in what we believe is a new way. The approach describes the general behavior of the beams to the level of providing interpretations of the Gouy phase (whose significance has long been disputed) [33][34][35][36] and geometric (Pancharatnam-Berry) phase [37][38][39], in a way that reveals the hidden geometry behind the transverse spatial structure of these familiar light beams. Furthermore, well-established methods of approximating the wave fields from the ray family are highly efficient for this approach, and even give the analytic forms for HG and LG beams. In a way, it forms a more complete and intuitive approach to our operator-based description of Gaussian beams in [22]. Our emphasis throughout is in recasting known properties of Gaussian beam families in terms of rays; the methods can be readily adapted as a design tool for new kinds of structured light.
The structure of this paper proceeds as follows. In the next section we discuss elliptic families of rays and show how they are associated with a Poincaré sphere. This is followed in Section III by a discussion on their quantization, and in Section IV by their geometric representation. Families of these ellipses and their quantization are discussed in Section V, which are combined to give a general method of constructing approximate wave solutions (Section VI), which are then applied to the HG, LG and GHL beams (Section VII). Properties such as Gouy and geometric phases (Section VIII) and the generalization to other beams families such as Bessel and Airy beams (Section IX) follow, before a concluding discussion. Many additional proofs and derivations are presented in the Appendices.

II. ELLIPTIC ORBITS AND THE POINCARÉ SPHERE
In geometric optics, the complex amplitude function representing a propagating, coherent monochromatic scalar light field is associated with a two-parameter family of light rays. For fields with slowly-varying intensities such as plane or spherical waves, the rays are normal to the wavefronts and intensity is proportional to the ray density. However, for fields with more spatial structure, the ray-wave connection is more subtle, as several rays may pass through a given point, albeit with different directions. In general, families of rays are bounded by envelopes known as caustics [40,41]. Interference fringelike structures can be caused by overlapping sets of rays, propagating in different directions. Near caustics or other features of structured light, the rays can differ significantly from the wavefront normals, but there is still a tight link between the wave and ray descriptions, and it is possible to accurately reconstruct the wave field by associating a complex contribution to the rays [42]. We will see how the geometry of structured Gaussian beams can be readily understood using rays.
We assume the beam propagates in a uniform, lin-ear and isotropic medium so the rays are straight lines. Each ray is determined by the transverse coordinate Q = (Q x , Q y ) where it crosses the z = 0 plane and its direction P = (P x , P y ). The equation for the point where the ray crosses the plane of constant z is therefore Q+zP. For a beam to be self-similar on propagation, the arrangement of the rays should be the same (apart from overall scaling) as z increases. As in Figure 1(a), both the shape and the orientation of the elliptic cross-section of a ruled hyperboloid are unchanged on propagation; endowing each ray in this one-dimensional subfamily with the same amplitude indeed guarantees its self-similarity on propagation. The two-parameter family of rays making up structured Gaussian beams are therefore a oneparameter superfamily of elliptic families of rays. We will call each such elliptic family an orbit of rays. In the quantification of these elliptic orbits, it is very convenient to use the parametrization of oriented ellipses afforded by the Poincaré sphere, borrowing language from polarization optics. (a) The straight rays sweep out a hyperboloid whose cross sections at any z are ellipses of the same eccentricity and orientation. The green curve is a normal to the rays. The length of the orange ray segment must be an integer multiple of the wavelength. (b) The eccentricity and orientation of the ellipse correspond to a point on the Poincaré sphere.
The Poincaré sphere for polarization parametrizes the two-dimensional complex Jones vectors v satisfying v * · v = 1, and v and v exp(−iτ ) are associated with the same polarization state for any real τ . v is defined in terms of latitude θ (Nb. not colatitude) and azimuth φ on the unit Poincaré sphere, (1) where − 1 2 π ≤ θ ≤ 1 2 π and 0 ≤ φ < 2π. As τ varies, 0 ≤ τ < 2π, [v exp(−iτ )] traces out the ellipse. In polarization optics, v represents the transverse, harmonic electric field; with τ evolving as time, the real part gives the ellipse; for each τ , the imaginary part is the velocity of the electric field vector [22].
For the hyperboloidal orbit of rays, θ and φ are fixed parameters determining the eccentricity and orientation of the elliptic cross-section. A given ray, labelled by τ , has z = 0 position where the constant Q 0 sets the transverse scale. The ellipse's major and minor semi-axes are Q 0 cos 1 2 θ and Q 0 | sin 1 2 θ|, and its foci are f ± = ±Q 0 cos 1/2 θ (cos 1 2 φ, sin 1 2 φ). The ray's directionthat is, its transverse velocity-is given by the imaginary part where P 0 is a constant determining the beam's angular divergence. Thus, at any z, the transverse ray coordinates given by Q + zP trace the same ellipse, up to a global hyperbolic scaling: where ζ = arctan(zP 0 /Q 0 ). On the ellipse, the position of each ray changes with z (hence "orbit"), but the orientation and eccentricity are unchanged, as shown in Fig. 1(a). We stress that the parametrization of elliptic orbits of rays by a Poincaré sphere is different physically from polarization. The similarity originates from the fact that mathematically, the Poincaré sphere parametrizes the classical orbits of the isotropic two-dimensional harmonic oscillator (like a transversely oscillating monochromatic electric field). Less obviously, ray families propagating back-and-forth in laser cavities also behave like classical harmonic oscillators, as the curvature of the spherical mirrors effectively acts as an attractive harmonic potential for the rays. Structured Gaussian beams are made up of families of orbits described by paths on the Poincaré sphere. First we discuss how the ray family is made compatible with the wave picture by a semiclassical 'quantization condition'.

III. QUANTIZATION CONDITION FOR THE ORBITS
Making the ray families consistent with wave optics requires two closure conditions dictated by the field's wavelength λ. The allowed solutions with certain properties (such as quantized OAM) are discrete, and often can be expressed as eigenfunctions of certain operators. These conditions are mathematically analogous to those in quantum mechanics, so we refer to them as quantization conditions. The first condition applies to the orbits. Since the rays in an orbit are skewed, a curve normal to them does not close onto itself after tracing the orbit (such as the thick green curve in Fig. 1(a)). There is a path difference along a ray between the initial and final points (represented by the orange line segment in Fig. 1(a)). Since optical path length (OPL) times wavenumber corresponds to the phase associated to a ray, this path difference must be an integer multiple of the wavelength for the ray and wave pictures to be consistent.
This condition is expressed mathematically as follows. Let L 1 (τ ) represent the OPL along each ray in the orbit, from some reference surface normal to the rays up to the z = 0 plane. The rays' inclination is determined by P(τ ), so L 1 depends on τ as dL 1 = P · dQ. From (2) and (3), the OPL difference between any pair of rays labelled by τ 1 and τ 2 is then After tracing the entire orbit, the total OPL mismatch is L 1 (2π) − L 1 (0) = πQ 0 P 0 , so the quantization condition yields where N is a nonnegative integer. Significantly, this condition does not involve θ and φ. Since Q 0 and P 0 describe the waist size and directional spread of the beam respectively, Q 0 P 0 π/λ is the beam quality factor M 2 [43][44][45], usually defined as the ratio of the product of the spatial and directional widths of a beam to the same product for a fundamental Gaussian beam. Therefore, the beam quality factor of fields made up of orbits satisfying (6) is quantized according to M 2 = N + 1. As we will discuss later, this index is also proportional to the beam's Gouy phase shift. In the quantum mechanical analogy, the integral (5) plays the role of the semiclassical Bohr-Sommerfeld integral, whose quantization (6) corresponds to energy quantization. With light beams constructed from orbits satisfying (6) with the same N guarantees the profile will have a well-defined beam quality factor and Gouy phase, as the beam is constructed to be an eigenfunction of the corresponding Hamiltonian operator.

IV. POINCARÉ AND PHYSICAL DISKS
The shape of an elliptical orbit depends on |θ|, with sign θ determining the sense of twist of the rays around the ellipse under propagation in z (that is, the sign of the OAM of the orbit Q × P, which is positive (counterclockwise) in Figure 1(a)). The beam's intensity profile is independent of this sign, so it is convenient to project both hemispheres of the Poincaré sphere, θ ≷ 0, onto the unit Poincaré equatorial disk (PED) 1 , with coordinates s = (s 1 , s 2 ) = cos θ(cos φ, sin φ), and |s| 2 ≤ 1. In the real space describing the transverse plane of the beam, we define the normalized ray position q = Q/Q 0 . The orbits are constrained to the interior of the unit disk |q| 2 ≤ 1, which we call the physical disk, since it is a scaled version of a cross-section of the beam (for any z). As Fig. 2 shows, a point s in the PED maps to an ellipse in the physical disk with foci f ± = ± √ cos θ (cos 1 2 φ, sin 1 2 φ), whose size is such that any rectangle in which it is inscribed is itself inscribed in the unit circle.

FIG. 2.
A point s in the PED maps to an ellipse with foci f± in the physical disk. The ellipse's major and minor axes have lengths 2 cos 1 2 θ and 2| sin 1 2 θ| respectively, equal to the sides of the gray rectangle. Note that any rectangle in which the ellipse is inscribed is itself inscribed in the unit circle.
The mapping between the PED and physical disk can be appreciated mathematically by considering each as the unit disk in the complex plane, so any real vector z = (z x , z y ) corresponds to the complex number Z(z) = z x + iz y = (1, i) · z. The complex numbers corresponding to the ellipse foci f ± = ± √ cos θ(cos 1 2 φ, sin 1 2 φ) are then the two square roots of the PED coordinate s = cos θ(cos φ, sin φ), as shown in Fig. 2, This map is conformal (angle preserving) except at the origin, as shown in Fig. 2: a Cartesian grid over the PED maps onto a curvilinear orthogonal grid over the physical disk.

V. FAMILIES OF ORBITS, CAUSTICS AND SOLID ANGLE QUANTIZATION CONDITION
The complete two-parameter ray family is constructed as a continuous one-parameter set of orbits. For the global ray structure to be preserved on propagation, all orbits must be coaxial, share a waist plane, and have common Q 0 and P 0 (and hence N ), so that they all scale as (Q 2 0 + z 2 P 2 0 ) 1/2 . Such a set of orbits corresponds to a path on the Poincaré sphere, which we call a Poincaré path. The cases of interest here are those whosepaths are closed loops. For simplicity, consider first a Poincaré path confined to a hemisphere, so that its projection onto the PED is a closed loop that does not touch the disk's edge. There are two Poincaré paths (one in the upper hemisphere, one in the lower) projecting to each such PED loop, as shown in Fig. 3(a). Each point on the projected Poincaré path corresponds to an ellipse in the physical disk, so the complete closed path gives rise to a family of ellipses, as shown in Fig. 3 show how the shape of their transverse ray structures is preserved under propagation (up to a hyperbolic scaling), and that this structure is the same when the loop is in the (c) upper or (d) lower hemisphere; the hemisphere only determines the handedness (the sign of the OAM). Fig. 3(b) represents the beam as a superposition of elliptical ray orbits. This structure is determined by the path's projection onto the PED, which we will also refer to as the Poincaré path. The envelopes of the family are caustics, here an outer one enclosing all the rays, and an inner one inside of which there are no rays. The brightest intensity features of a beam are associated with these caustics, as the density of rays is highest near them. There are surprisingly simple geometrical relations between the projected Poincaré path and the caustics in the physical disk, which we now describe (the derivation can be found in Appendix A).
The geometric relation is easiest to appreciate for a Poincaré path in the PED with endpoints at the edge of the disk, such as the one shown in Fig. 4(a). The corre- FIG. 4. Medial axes of Poincaré paths in the PED map to caustics in the physical disk. Given a Poincaré path (thick black curve), one can find two medial axes as the loci of the centers of circles that touch this curve and the unit circle. The mapping Z(q) = ± Z(t), where t are points along the medial axes, corresponds to curves of points q that are the caustics of the resulting fields. Note that the caustics in (b) correspond to those in Fig. 3 sponding path on the full Poincaré sphere is symmetric in the upper and lower hemispheres (projecting to the same curve in the PED). The geometric prescription for finding the caustics is then as follows: 1) Find the set of circles that are tangent to both the Poincaré path and the unit circle. There are two such sets on each side of the Poincaré path (shown in pale blue and red in Fig. 4(a)). The centers of each set of circles define a curve equidistant from the Poincaré path and the unit circle. Each such curve, being equally close to two other curves, is a medial axis [47] (or topological skeleton), in terminology borrowed from image analysis. 2) Find the caustics by applying the square root map z(±Z 1/2 (t)) to each medial axis. Given that the square root maps each point in the PED onto two points on the physical disk, each medial axis is mapped onto two caustic segments that are identical except for a 180 • rotation about the origin. Therefore, each point along the Poincaré path gives rise to two medial axis points and therefore to four caustic points. In cases like that shown in Fig. 4(a), where the Poincaré path begins and ends at the unit circle, the two medial axes meet at the same endpoints, intersecting each other at right angles. Since the square-root mapping is confor-mal, the caustics in the physical disk also intersect at the disk's edge at right angles.
This construction is also valid for the previous case where a loop is fully within one hemisphere of the Poincaré sphere, as shown in Fig. 4(b) for the same Poincaré path as Fig. 3. Each of the medial axes is now a closed loop, as are the caustics (their square-root images). The outer medial axis (orange line), formed by the centers of the (red) circles, is constrained to the annular space between the Poincaré path and the unit circle. The inner medial axis (blue line) is formed by the centers of the (pale blue) circles that touch the inside of both the unit circle and the Poincaré path. If, as in this example, there are sufficiently small radii of curvature at some points of the Poincaré path, the inner medial axis (and its corresponding caustic) can cross itself and have cusps.
The geometric connection between the path in abstract Poincaré space and the beam's caustics in physical space is one of the main results of this work. It implies that the caustics of a structured Gaussian beam are composed of two parts that are not mutually independent: one can either prescribe a Poincaré path and determine the caustics via the medial axes, or instead prescribe one caustic (with the constraint that it must be symmetric under rotations by 180 degrees), then find the corresponding medial axis in the PED, and thus the Poincaré path, and then the second medial axis and caustic.
The Poincaré path is parametrized as s(η) = (cos θ(η) cos φ(η), cos θ(η) sin φ(η)), so overall the ray family is parametrized by 0 ≤ τ, η < 2π, topologically corresponding to a torus. Families of ellipses in the physical disk such as in Fig. 3(b) are projections of this torus, with its outline given by the caustics, consisting of either an outer and an inner loop (as Fig. 4(a)), or a quadrangle with corners at the boundary of the physical disk (as Fig. 4(b)). The OPL at z = 0 for all the rays, in terms of τ and η, may be found in Section S2.
As the Poincaré path is a closed loop, wave-optical self-consistency requires that any physical quantity (determined by OPL) must return to its starting point on a circuit of η. This gives a quantization condition around the path, just as our previous condition quantized the orbits. This condition, whose derivation is in Section S3, is geometrically remarkably simple: the solid angle Ω on the Poincaré sphere enclosed by the Poincaré path must be an odd multiple of 2π/(N + 1), namely where N/2 denotes the integer part of N/2. We may appreciate the significance of this by referring to the quantum mechanical picture. Structured light beams are usually considered as eigenfunctions of some optical operator, such as the OAM operator L = −i∂/∂φ giving the LG modes [3], or the astigmatism operator M = 1 2 (−∂ 2 /∂x 2 + ∂ 2 /∂y 2 + x 2 − y 2 ) giving the HG modes [22]. In the completely classical, Hamiltonian picture, these quantities are functions of position Q and momentum P, which define families of contours on the Poincaré sphere (the sphere of orbits of the isotropic two-dimensional oscillator). Thus the angular momentum L is simply the height coordinate of the Poincaré sphere cos θ, and for M it is the horizontal coordinate sin θ cos φ [22]; the contours are then circles concentric to the vertical or horizontal axes of the sphere. The condition (8) picks out a discrete set of these contours as the Poincaré paths, which correspond to the sets of ray families which are consistent with wave optics. We will discuss the LG and HG modes in much more detail, after we have discussed how to construct approximations to the wave fields from the appropriately quantized ray families.

VI. RAY-BASED WAVE FIELD RECONSTRUCTION
There are many methods for estimating wave fields based purely on a ray description, which are valid even in the presence of caustics. We here use an approach [48][49][50][51] in which a complex Gaussian field contribution is assigned to each ray, and the estimate takes the form of a double integral over τ and η. It is shown in Section S5 that the integral in τ can be evaluated analytically, leading to a field estimate at the waist plane of the form where x = (x, y) is the transverse position at the waist plane, A(η) is a non-negative amplitude function weighting the different orbits, v is the Jones vector in (1) parametrized in terms of η, T (η) is given in Section S3, and U N are Hermite-Gaussian elementary fields evaluated at complex values, defined as where H N is the N th order Hermite polynomial and v·v = cos θ. Up to a complex factor, U N is the wave contribution corresponding to an elliptical ray orbit specified by the Jones vector v. Fig. 5 shows, for several choices of v, the real part and intensity of each of these orbit contributions, together with the corresponding ray-optical orbit. N is evidently the number of phase oscillations around the ellipse. In fact, these elementary field contributions are themselves a subset of the HLG beams, which are associated with points over a Poincaré sphere [22][23][24][25]27]. However, these contributions are expressed not as a superposition of HG or LG beams but as a single term involving a Hermite polynomial evaluated at a complex argument proportional to the Jones vector. The expression in (10) provides a general prescription for constructing self-similar beams that are rigorous solutions to the paraxial wave equation, and that have caustics at prescribed locations.

VII. EXAMPLES: LG, HG AND HLG BEAMS
We now illustrate these ideas for the two most common families of beams of this type, LG and HG beams, as well as for the more general HLG beams. LG beams are separable in polar coordinates, and their ray structure was studied by Berry and McDonald [52]. The LG Poincaré path in the PED is a circle of radius r centered at the origin (so θ = arccos r). The solid angle enclosed by this circle over the Poincaré sphere is Ω = 2π(1−sin θ) = 2π[1−(1−r 2 ) 1/2 ], which is quantized according to (8), such that r can only take the values for n = 0, 1, ..., N/2 . The medial axes, equidistant from the unit circle and the Poincaré path, must also be circles centered at the origin, but with radii (1±r)/2. Following the square root map onto the physical disk, the two caustics are circular as well, with radii Q 0 (1 ± r)/2. More details about the ray description of these beams are given in Section S6, where it is also shown that, remarkably, the wave field estimate in (9) actually yields the exact form for LG beams. Fig. 6 shows the Poincaré and physical discs for these beams, including the Poincaré path, medial axes, caustics, and elliptical orbits, as well as the ray structure of the beam and the intensity cross section.
In terms of the operator picture, all physical quantities (Poincaré path, medial axes, caustics) must be rotation invariant, and the path quantization gives the usual angular momenta quanta −N ≤ ≤ N , quantized in integers (in steps of 2).
We now consider HG beams, which are separable in Cartesian coordinates. The Poincaré path is a straight line terminating at the edge of the PED, as shown in Fig. 7(a). Since the PED is a projection of the sphere onto its equatorial plane, the Poincaré path on the sphere is a circle centred at the s 1 axis with radius r, equal to half the length of the straight line, is also quantized according to (11) [22]. By simple geometry, both medial axes are confocal parabolas with foci at the origin, which intersect each other and the Poincaré path at the edge of the disk. These parabolic medial axes are shown in Fig. 7(a) as blue and orange curves. The caustics (square roots of the parabolas) are straight lines, as shown in Fig. 7 Finally, we consider HLG beams, which are realized by conversion of HG or LG beams through simple combinations of cylindrical lenses [30,32] or equivalent SLM implementations, which amount to rotations of the Poincaré sphere about an axis in the equatorial plane, but which cannot be expressed simply in any separable coordinate system. As for LG and HG beams, for HLG beams the Poincaré path on the sphere is a (planar) circle whose radius r is quantized according to (11). However, the centre of this circle can be at any angle β with respect to the vertical s 3 axis. We assume for simplicity that the centre lies in the s 1 s 3 plane, so that β = 0 gives LG beams while β = π/2 reduces to HG beams separable in x and y. Projected onto the equatorial disk, the Poincaré path is an ellipse centered at ([1−r 2 ] 1/2 sin β, 0) and with minor and major semi-axes given by r cos β and r respectively, as shown in the left column of Fig. 8. The medial axes and hence the caustics (shown in the figure's second column) can be found in parametric form, and are not conic sections. Similarly, the wave fields are no longer separable in a coordinate system, but they can still be In the operator picture, these beams are eigenfunctions of L cos β + M sin β, which have the same integer eigenvalues by the hidden symmetry of the isotropic 2dimensional harmonic oscillator (corresponding, in the ray picture, to rotating the spherical cap in the s 1 s 3 plane). Since the operator is linear in the coordinates of the Poincaré sphere's space, the Poincaré path is a circle with uniform weight. This simplicity of the HLG family explains why the ray-based field estimate in (9) actually yields the known exact eigenstates of the operators.

VIII. GOUY AND PANCHARATNAM-BERRY PHASES
In addition to revealing the hidden geometry behind the caustic structure of structured Gaussian beams, the description presented here provides a simple, ray-based explanation for their Gouy and Pancharatnam-Berry phase shifts. These two phase shifts turn out to correspond to shifts in each of the two ray parameters, τ and η, as follows.
Consider first the case of the Gouy phase shift. As shown in (4), propagation in z preserves the ray structure up to a shift τ → τ + ζ, where ζ = arctan(zP 0 /Q 0 ). Thus, any ray initially at a given location when z = 0 is replaced, after propagation, by another one from the same orbit, whose value of τ is larger by an amount ζ. Since a variation in τ of 2π corresponds to a path length of (N + 1)λ, this shift in τ by ζ amounts to a change in path length of (ζ/2π)(N + 1)λ, and hence to a phase of (N + 1) arctan(zP 0 /Q 0 ), namely the standard Gouy phase for a beam of this type. This effect can be appreciated from Fig. 1: all rays have roughly the same length. However, the ray that touches a given point it the orbit (say, a vertex of the ellipse) at the initial plane is not the same as the one that touches the same point at the final plane. The total phase difference is then not only due to the length of the rays but also to the OPL difference between the two rays in question. The geometric phase for beams of the HLG family under astigmatic transformations has been studied in algebraic terms by exploiting the analogy with 2D quantum harmonic oscillators [25][26][27][28] and verified experimentally for low-order beams [37][38][39]. Consider subjecting a HG, LG, or more general HLG beam to a series of optical transformations that rotate the Poincaré sphere around an axis within the s 1 s 2 plane (through a suitable combination of cylindrical lenses) or around the s 3 axis (through a beam rotator such as a pair of Dove prisms or periscopes). By choosing the sequence of transformations appropriately, the circular Poincaré path for the beam can be brought back to its initial position after its center traced a trajectory over the Poincaré sphere. However, it is easy to see that each point within the Poincaré path does not necessarily fall back onto its initial position; rather, the final state of the circle is generally rotated around its axis with respect to the initial one by some angle Θ, depending on the trajectory followed. If this trajectory is composed only of segments of great circles, as that shown in Fig. 9, then the angle Θ equals the solid angle subtended by the path. In other words, this transformation reduces to a shift η → η − Θ. Such rotation results in a phase shift for the beam that can be considered as a geometric phase, because it is not related to a change in the OPL of each ray, but to a cycling of the roles that different rays (and indeed orbits) play within the pattern. As stated in Section S3, the phase due to a complete rotation of the Poincaré loop is k∆L 2 = (N − 2n)π, so the corresponding geometric phase is (N − 2n)Θ/2 = Θ/2, where is the OAM label of the LG beam within the set.
In summary, the phase accumulated under propagation can be separated into a "dynamic" phase, due to the path length traced by each ray, and a Gouy phase, due to the cycling of rays within each orbit. If additionally the beam is subjected to a series of transformations that rotate the Poincaré sphere but that bring the beam back to its original shape, there is a third, geometric component of the phase, due to the shifting of orbits within the beam structure. However, note that while the dynamic and Gouy phases apply to any beam, the general geometric phase can only be achieved for HLG beams, given the rotational symmetry of their Poincaré path. For beams whose Poincaré paths have M -fold symmetry around an axis, a more restricted version of the same phenomenon is possible. FIG. 9. Illustration of geometric phase as a cycling of orbits. Part (a) shows the transformation of a LG beam with positive OAM transformed onto one with negative OAM by following a meridional path in the s1s3 plane over the Poincaré sphere. Five stages of this path are shown explicitly, including a HG beam at the equator. Part (b) shows the transformation of the beam back onto its initial configuration through a different meridional path. Note that the ray configurations are rotated by π/4 with respect to those on the left. While the final beam has the same shape as the initial one, the orbits (identified by color) are rotated, resulting on a geometric phase.

IX. OTHER SEPARABLE SELF-SIMILAR BEAMS AS LIMITING CASES, AND SELF-HEALING
Despite them not being explicitly Gaussian, other types of propagation-invariant beams, such as Bessel [6,7], Mathieu [8], Airy [4,5], and parabolic [54] beams, correspond to limits of the structured Gaussian beams described here. These other beams are idealized solutions that involve infinite power, corresponding to the limit N → ∞ in particular regions of the physical disk. That is, the ray families are open rather than closed loops.
Bessel and Mathieu beams correspond to a small neighbourhood of the origin of the physical disk, and the outer radius of the PED. For Bessel beams, the Poincaré path is a circle centered at the origin and whose radius is nearly equal (or equal) to unity, so that one medial axis is a small circle (or point) centered at the origin, and so is the inner caustic. Mathieu beams use the same construction, except that the large circular Poincaré path is shifted slightly from the origin but still fits within the PED. This shift de-centers the small inner medial axis, and causes the resulting inner caustic to be elliptic. (A shift larger than the difference between unity and the path's radius would make the inner caustic hyperbolic.) If instead the centred Poincaré path is slightly deformed into an ellipse, the inner caustic becomes an astroid, as in beams produced by misaligned axicons [53].
In the case of Airy and parabolic beams, on the other hand, one must focus at a small region at the edge of the PED and physical disk. Airy beams are the limit of the intersection of two medial axes (and caustics) when a locally straight Poincaré path touches the edge of the PED. This geometry is shown in Figure 10(a). The angle of this intersection determines the ratio of the spacing of the intensity lobes along the two caustic sheets. Parabolic beams result when the Poincaré path is a very small circular segment starting and ending at the edge of the PED, leading to two sets of parabolic caustics, as shown in Fig. 10(b). Propagating self-similar beams are often referred to as "self healing"; if an obstacle blocks a limited part of the beam in one plane, the blocked intensity features reappear as z increases. The effect of the block can be described to first order in terms of the ray-optical shadow projected by the obstacle, that is, the suppression of a subset of the rays composing the field. Self healing (which can occur more generally [55,56]) is then easily explained in terms of the cycling of rays within each ray orbit under propagation: the blocked rays are replaced by other rays leading to the same local ray structure. However, it is clear that rather than "healing", the beam's "wound" is simply transferred to a different part of its transverse profile. For beams such as Airy or Bessel beams, the idealized ray family is open, so the shadow is ultimately lost in an infinite reservoir of rays away from the region where the main intensity features are located. This is not the case for structured Gaussian beams whose ray family is compact.
Due to the rotational symmetry of their Poincaré path, HG, LG, and more general HLG beams can undergo local "healing" not only through shifts in τ under propagation but also through shifts in η due to rotations of the Poincaré sphere caused by the optical transformations discussed in Section VIII. Such shifts would have a similar effect of displacing the blocked regions within the beam's profile. By further abusing the already imperfect "healing" metaphor, this effect could be called "assisted healing".

X. CONCLUDING REMARKS
We proposed a ray-based description of structured Gaussian beams that reveals hidden geometrical restrictions in their spatial structure, particularly their caustics. Further, the Gouy and geometric phases that can be accumulated under propagation were also given simple explanations in terms of rays and their quantization. The description given here is based on the partition of the two-parameter ray family, one parameter giving rays around orbits with an elliptical cross-section, and the other defining a curve on the Poincaré sphere representing the elliptic ray family. This develops previous work also employing Poincaré spheres to characterize the modal structure of HLG beams [23-28, 37, 38]. However, unlike these previous studies, where each beam is associated with a point on the Poincaré sphere, in our more general construction the beam is associated with a curve on the Poincaré sphere. The shape of this extended curve not only determines the shape of the beam but also explains (and restricts) the geometric phase resulting from beam transformations.
The approach given here also differs from other raybased studies of structured Gaussian beams. For example, Gaussian beams have been described as bundles of complex rays [57][58][59], as opposed to the real rays used here. Similarly, ray-like descriptions of LG and Bessel beams have been given in terms of Wigner functions [60], but such a description uses all rays in phase space rather than a two-parameter family, so the concept of caustic is absent and the representation in the Poincaré sphere is not compatible with that treatment. Finally, descriptions also exist in terms of curved flux lines rather than rays [52].
In the complementary operator picture of our approach [22], there is a spin vector-like operator on the Poincaré sphere for which these HLG beams, described by circles whose centers are given by the vector direction of the operator. The operator approach, built around the su(2) Poisson algebra of the two-dimensional harmonic oscillator, reveals the algebraic connection between structured Gaussian beams, the classical and quantum harmonic oscillator and the Poincaré sphere, contrasting with the semiclassical approach here.
Although our focus here has been the particular examples of HG, LG and HLG beams, the ray-based approach can be applied to any beam with a Gaussian envelope with a well-defined Gouy phase -in fact, the approach allows such structured Gaussian beams to be designed from almost arbitrary paths on the Poincaré sphere satisfying the quantization condition (8). One obvious structured Gaussian family we have not explored here is the Ince-Gaussian beams [61,62], which also interpolate between HG and LG beams but which are separable in elliptic coordinates. From the other separable beams considered here, one might expect the caustics of Ince-Gaussian ray families to be confocal ellipses and hyperbolas. Indeed this is the case, and we defer a full discussion to a later article. Although no other separable Gaussian beam families exist [63], the freedom of choice of curves on the Poincaré sphere allows a huge variety of Gaussian beams with new and unfamiliar properties to be designed.

FUNDING INFORMATION
The Leverhulme Trust (MRD), and the National Science Foundation (PHY 1507278) (MAA).
It was shown in Section 4 of the main text that a point s in the PED corresponds to an ellipse (a curve) in the physical disk through a square-root mapping onto the ellipse's foci and a normalization condition for the ellipse's size. Let us now study the reciprocal map: if we prescribe a point q in the physical disk, what is the curve in the PED composed of all points s corresponding to ellipses that contain q? It turns out that this curve is a circle, as we now show. Figure 1 shows a point q = (q cos α, q sin α) within the physical disk, as well as several elliptical orbits that contain it. A representative orbit is highlighted in blue. We now use the fact that the sum of the distances (labeled by red lines in Fig. 1) from q (or any other point along the ellipse) to each of the foci equals the long axis of the ellipse (in this case 2 cos φ/2), that is, After some manipulation and using the fact that |f ± | 2 = cos θ, this equation can be reduced to By using s = (cos θ cos φ, cos θ sin φ), we can write this expression as where t = (q 2 cos 2α, q 2 sin 2α). That is, the curve in the PED described by the points s corresponding to all the ellipses that include a given point q in the physical disk is a circle centered at the point t, which is related to q through a quadratic conformal map of the form Further, since the radius of the circle is one minus the distance of its center to the origin, the circle always touches the edge of the PED (a unit disk), as shown in Fig. 1.
Supplementary Figure 1. The points in the PED corresponding to all ellipses crossing a prescribed point q in the physical disk trace out a circle that touches the edge of the PED and that is centered at the point t, which is given by the quadratic conformal map of q following Eq. (A4).
The above construction allows answering the question of which infinitesimal Poincaré path segment gives rise to a desired caustic segment. Consider an infinitesimal caustic segment joining two points q 1 and q 1 + dq 1 over the physical disk, shown in Fig. 2(a). These two points correspond to two circles that touch the edge of the PED and that are centered at points t 1 and t 1 + dt 1 , respectively. Note that, given that the quadratic map from q 1 to t 1 in Eq. (A4) is conformal except at the origin, the angle β between q 1 and dq 1 is the same as that between t 1 and dt 1 . The circles centered at t 1 and t 1 +dt 1 , respectively, intersect at two points, one of which is always at the edge of the disk. The second intersection, labeled as s in Fig. 2(a), is at the mirror image of the first intersection with respect to the line containing t 1 and dt 1 . Therefore, this second intersection is located at By construction, s is then the point in the PED associated with an ellipse in the physical disk that contains the desired caustic segment, i.e., it is a point along the Poincaré path. However, it is not sufficient to prescribe s; the local direction of the Poincaré path at s must also be specified, so that the ellipses corresponding to neighboring point are nearly tangent to the caustic segment too. This local direction is simply the direction of the circle segments at the intersection, which from Fig. 2(a) is seen to correspond to the angle 2(α − β) + π/2 with respect to the s 1 (horizontal) axis.
Consider now the case where an infinitesimal segment of the Poincaré path at a point s is prescribed. Figure 2(b) shows that there are two circles that touch tangentially both this segment and the unit circle, whose Supplementary Figure 2. (a) Two endpoints of an infinitesimal caustic segment, q1 and q1 +dq1, map onto points t1 and t1 + dt1 which are the centers of two circles whose intersection location s and direction define an infinitesimal segment of the Poincaré path. (b) Given a point s and local direction of the Poincaré path, there are two circles that are tangent to it and that touch the unit circle. The centers of these circles, t1 and t2, map onto four caustic points in the Physical disk, ±q1 and ±q2, according to Eq. (A4).
centers are labelled as t 1 and t 2 . Each of these points is related to a different caustic, and maps (under the square root map) to two symmetrically distributed points along it, ±q 1 and ±q 2 . As we move along the Poincaré path, the points t 1 and t 2 (which are the circles constrained to touch the unit circle and the Poincaré path) trace what is known as a medial axis. The caustics result from applying the square root mapping to these two medial axes.

Appendix B: OPL for the complete ray family
We now derive the parametrized OPL for the complete ray family. Let the Poincaré path be parametrized as s(η) = [cos θ(η) cos φ(η), cos θ(η) sin φ(η)], and let L 2 (η) be the OPL of a subset of rays along some curve in the physical disk that touches all the orbits. The increment in OPL for the rays along this curve is again given by dL 2 = P · dQ. By using Q = Q 0 q and P = P 0 ∂q/∂τ , this increment can be written as It is convenient to integrate this equation along one of the caustics, so that ∂q/∂η is parallel to ∂q/∂τ and their dot product equals the product of the magnitudes. Also, from Eqs. (1) and (3) we find the relation |P/P 0 | 2 + |Q/Q 0 | 2 = 1, from which we get |∂q/∂τ | = 1 − |q| 2 . With this, Eq. (B1) can be simplified to for i = 1, 2. This expression can be mapped onto the PED since points q i along a caustic map onto points t i along a medial axis. Since this mapping is quadratic and conformal, it is easy to see that |∂q i /∂η| = |∂t i /∂η|/2 |t i |, yielding so that L 2 can be written as an integral along a medial axis as The function that describes the OPL for any ray in the family is then equal to the OPL up to the orbit containing the ray plus the OPL between the ray in question and the ray in the orbit touching the caustic: where T (η) is the value of τ for the orbit labeled by η corresponding to the ray that touches the reference caustic. This value can be found as where R(φ/2) is a rotation matrix by an angle φ/2, q i (η) is the parametrized caustic, arctan(u) = arctan(u x , u y ) is the two-parameter arctangent function, and θ(η) and φ(η) follow from the specification of the Poincaré path s(η) = [cos θ(η) cos φ(η), cos θ(η) sin φ(η)].

Appendix C: Quantization of the enclosed solid angle
We now derive the second quantization condition for the ray family. For simplicity, consider a case where the Poincaré path traces a full loop within the PED, as in Fig. 4(b) of the main text. Let ∆L 2 be the total increment of optical path that results from tracing a complete medial axis (e.g., the thick orange loop inside the PED in Fig. 4(b)): It is shown in the next Section of this document that the quantity in parentheses equals the solid angle Ω enclosed between the Poincaré path and the equator over the surface of the Poincaré sphere. Therefore, by also recalling the first quantization condition in Eq. (6), Eq. (C1) can be written as where Ω is the solid angle enclosed by the Poincaré path. Given the square-root mapping, spanning the complete medial axis amounts to spanning only half of the corresponding caustic (say, the bottom-right half between the two marked dots of the orange caustic loop in the physical disk in Fig. 4(b)). That is, ∆L 2 is the optical path difference between two points along the caustic that are symmetrically located around the origin, and that therefore belong to the same orbit (indicated by the thicker green ellipse in the figure). Therefore, the phase difference due to ∆L 2 must be consistent with that calculated from the optical path difference of the corresponding endpoint rays along the orbit they belong to (i.e., along one half of the thick green ellipse joining the two orange points in Fig. 4(b)), given from Eqs. (5) and (6) by L 1 (τ 0 +π)−L 1 (τ 0 ) = (N +1)λ/2. Note, however, that when calculating the phase difference from the orbit, one must add by hand an extra phase term resulting from the fact that the ray family touches caustics between the two points. Phases resulting from caustics are known as Maslov phases, and have a magnitude of π/2. In this case, the total phase due to caustics is π, since the elliptical segment between the two rays in question touches two caustics ("half" at each contact point indicated by an orange dot with the caustic in question, and one inbetween with the other caustic, indicated by a pale blue dot). The phases along the caustic and the orbit are then consistent if where m is an integer. In other words, after using Eq. (6) and Eq. (C1) of this document, the enclosed solid angle Ω must take values: where n = m − 1. That is, the total solid angle enclosed by the Poincaré path must be an odd multiple of 2π/(N + 1).
Let us now show that this quantization for the enclosed solid angle also holds for cases where the path in the PED starts and finishes at the disk's edge, as in Fig. 4(a) Let the starting and finishing points correspond to values η 1 and η 2 of the parameter η. For simplicity, assume that this open path in the PED corresponds to the projection onto the equatorial plane of a closed loop over the surface of the Poincaré sphere that occupies both hemispheres and has mirror symmetry with respect to the equatorial plane. In this case, one must verify that the phase due to the optical path difference between the endpoints of the path segment corresponding to the s 3 ≥ 0 hemisphere is consistent with that for the path segment corresponding to the s 3 ≤ 0 hemisphere. Since the orbital angular momentum is reversed between these hemispheres, so is the accumulation of optical path difference. This means that L(η 2 ) − L(η 1 ) for the lower hemisphere has equal magnitude but opposite sign to the corresponding path difference for the upper hemisphere, so that the total optical path accumulated around the loop is ∆L 2 = 2L(η 2 ) − 2L(η 1 ). Since the closed loop (corresponding to one of the caustics) touches the other caustic twice (once at each equatorial point), an extra Maslov phase of π must be included. The quantization condition in this case can then be written as where m is an integer. However, L 2 (η 2 ) − L 2 (η 1 ) corresponds to Q 0 P 0 /4 times the solid angle within a hemisphere to one side of the loop (depending on which caustic was used), and therefore ∆L 2 = 2L 2 (η 2 ) − 2L 2 (η 1 ) = Q 0 P 0 /4(4π −Ω), where Ω is the solid angle enclosed/excluded by the symmetric loop. By now defining n = N − m + 1, we find that the quantization condition is indeed the one in Eq. (8). It is easy to see that the same condition would hold for a loop that occupies both hemispheres but is not symmetric. Finally, notice that whether the loop is restricted to a hemisphere or not, we have the relation k ∆L 2 = (N − 2n)π.

Appendix D: Solid angle interpretation of the integral
The equation for the circles centered at t that touch both the Poincaré path and the edge of the PED (the unit circle) is given in Eq. (A3). This equation can be parametrized as where u(ξ) is a unit vector at an angle ξ measured, for convenience, from the direction of t. As t runs along one of the medial axes, the circles advance touching both the unit circle and the Poincaré path, as shown in Fig. 4. The projections of these circles over the surface of the Poincaré sphere are three-dimensional curves that touch both the Poincaré curve and the equator. These curves are shown in yellow in Fig. 3, with one of them highlighted in black. Their projections onto the equatorial plane are the pale blue circles (one also highlighted in black), whose centers t trace the medial axis (blue curve). The solid angle enclosed between the Poincaré curve (thick black) and the equator can then be calculated as the solid angle subtended by either the leading or the trailing segments of the projections of the circles joining these curves as t runs through the medial axis. While using either the leading or trailing segments alone would give the same result for the integral over the whole medial axis, it would give different forms for the integrands. Using the average of the two options turns out to give a simpler expression, and more importantly, the one we set out to find.
We now calculate the average of the two solid angle elements corresponding to the changes in both the leading and trailing segments of the circle as t varies by an infinitesimal amount δt, at an angle ν with respect to t. The first step is to calculate the radial displacement of the points in the circle with respect to its center, resulting from this infinitesimal change. This radial displacement is given by As could be expected, this radial displacement vanishes for ξ equal zero and 2ν, which correspond to the inter-sections of the circles at the unit circle and the Poincaré curve, respectively, and which are shown as green and purple dots in the highlighted curve in Fig. 3. These two values are the boundaries between the leading and trailing segments of the circle. If we were to integrate this displacement in ξ times the radius (1 − |t|) over an interval of size 2π,, we would obtain the area of the two complementary crescents enclosed by the initial and final circles. However, what we want to find instead is the solid angle element projected over the surface by this area, which is given by the integral of the magnitude of this radial thickness times the length element (1 − |t|)dξ, divided by an obliquity factor h(ξ): where ξ 0 is an arbitrary limit of integration. Note that we use the absolute value of the radial displacement so that we add rather than subtract the solid angles spanned by the leading and trailing segments of the circles, and that a factor of 1/2 is included for performing an average. Given the spherical geometry, the obliquity factor h equals the height of the projection of the point in question onto the surface of the sphere: Substituting into the integral for the solid angle element, we get where the result in the integral in the last step is easiest to find by setting ξ 0 = 2ν. It is remarkable that this result is independent of ν, the local direction of the medial axis with respect to the radial direction. The main result of this Section is then the fact that the integral of the previous result over the complete medial axis equals the solid angle enclosed between the Poincaré curve and the equator: where η is a parameterization of t and the integral must extend over the complete medial axis.
Gaussian field contributions, each centered at a ray, according to where x = (x, y) are the transverse coordinates at the waist plane, Γ is a constant with units of inverse length that regulates the widths of the Gaussian contributions, and A(η) (chosen here to be independent of τ ) is a nonnegative amplitude function indicating the relative importance of the different orbits. The field estimate in Eq. (E1) is asymptotically insensitive (for large k) to the choice of the scale parameter Γ. However, for our current purposes, it is convenient to use Γ = P 0 /Q 0 (so that the Gaussian elements have the same width as the elementary Gaussian beam at the waist). The Jacobian in Eq. (E1) can then be found from Eqs. (2) and (3) to be where with v x and v y being the components of the Jones vector in Eq. (1). Equation (E1) can then be written as Here, I is the integral over the τ dependent parts, which can be solved analytically: where in the second step we used the change of variables u = exp(−iτ ) so that the integral is over the unit circle (u.c.) of the complex plane, in the third step we used the shorthand u 0 = 2x · v/[Q 0 cos θ], in the fourth step we applied residue theory, and in the last step we recognized the generating function of the N th order Hermite polynomial H N . With this, the field estimate can be written in the form in Eq. (9).
The ray parameters can then be expressed in terms of these coordinates as τ = arcsin 1 + r − 2q 2 2q 2 − 1 + r , (F7) η = ϕ − arctan (1 − r)(1 + r − 2q 2 ) (1 + r)(2q 2 − 1 + r) . (F8) From these expressions we see that the OPL in Eq. (F4) can be written as a sum of a part dependent purely on q and a part dependent purely on ϕ, pointing at the separability of the problem in polar coordinates. We now evaluate the field estimate in Eq. (9) by writing the normalized position vector x also in polar coordinates as x = (Q 0ρ , ϑ). By setting A(η) = 1 and defining = N − 2n, Eq. (9) simplifies to U (x)≈kP 0 exp i π 4 exp − N + 1 2 where η = η − ϑ and a ± = (N + 1)/r( √ 1 + r ± √ 1 − r)ρ/2. Notice that the range of integration was extended from [0, π] to [0, 2π] by inserting a factor of 1/2 and using the fact that has the same parity as N and hence the integrand is a combination of even powers of exp(iη ). This integral can then be solved by using the expansion for Hermite polynomials: where (ρ, ϑ) are polar coordinates for x and = N − 2n is the vorticity of the beam.

Appendix G: Hermite-Gaussian beams
For HG beams, the Poincaré path is a straight line, chosen here in the vertical direction for convenience. Since this line is the projection of a circle of radius r over the unit sphere's surface, we parametrize it as s = ( √ 1 − r 2 , r sin η). The angles determining the parametrized point in the Poincaré sphere are then θ = arccos 1 − r 2 cos 2 η , (G1) φ = arctan 1 − r 2 , r sin η .
While the general parametrization of the ray positions employed in this document so far (in which τ = 0 corresponds to a vertex of the ellipses) gives the correct result, for this case there is an equivalent but more convenient parametrization given by We can then change variables to τ x,y = τ ± η/2 (with unit Jacobian), so the position and momentum vectors simplify to Q = (Q 0 ν x cos τ x , Q 0 ν y cos τ y ) , (G10) P = (−P 0 ν x sin τ x , −P 0 ν y sin τ y ) , and the OPL is now easily found to be The fact that the total OPL is separable as a sum of a part that depends only on the x components and a part that depends only on the y components heralds the separability of the problem in Cartesian coordinates. This separability becomes evident when considering the wave field estimate in Eq. (9), where by assuming A(η) to be constant (set to unity for simplicity), one arrives at the exact result, which is of course separable. With Γ = P 0 /Q 0 , the Jacobian in Eq. (9) gives