Ray and caustic structure of Ince-Gauss beams

The Ince-Gauss beams, separable in elliptic coordinates, are studied through a ray-optical approach. Their ray structure can be represented over a Poincaré sphere by generalized Viviani curves (intersections of a cylinder and a sphere). This representation shows two topologically different regimes, in which the curve is composed of one or two loops. The overall beam shape is described by the ray caustics that delimit the beams’ bright regions. These caustics are inferred from the generalized Viviani curve through a geometric procedure that reveals connections with other physical systems and geometrical constructions. Depending on the regime, the caustics are composed either of two confocal ellipses or of segments of an ellipse and a hyperbola that are confocal. The weighting of the rays is shown to follow the two-mode meanfield Gross–Pitaevskii equations, which can be mapped to the equation of a simple pendulum. Finally, it is shown that the wave field can be accurately estimated from the ray description.


I. INTRODUCTION
The paraxial wave equation accepts many types of solutions, some of which are separable in a given coordinate system [1][2][3][4][5][6][7][8].A particular case is that of beams that take the form of a Gaussian times a polynomial involving the corresponding separation variables [1,2,9].The Hermite-Gauss (HG) and Laguerre-Gauss (LG) beams, for example, are separable in Cartesian and polar coordinates, respectively, and are used in a range of applications including quantum information, telecommunications and imaging [10][11][12][13][14][15].Given the simplicity of the closed-form expressions for these beams in the wave regime, only a handful of studies consider them in terms of a ray-optical picture [16][17][18].
The Ince-Gauss (IG) beams [9,[19][20][21][22][23][24], separable in elliptic coordinates, have received significant attention in the last couple of decades, finding also many applications [25][26][27][28].Elliptical coordinates require specifying the separation 2f between its two foci.Let the radial ξ ∈ [0, ∞) and angular ν ∈ [0, 2π) elliptic variables be defined as r = (x, y) = (f cosh ξ cos ν, f sinh ξ sin ν).Within a region of some length scale w 0 around the origin, these elliptic coordinates tend to the Cartesian coordinates as f → ∞, while for f → 0 the contours of this system reduce to radial lines and circles as in polar coordinates.IG beams are then a broader family of beams that include both HG and LG beams as its two limiting cases.The specific shape of an IG beam depends on the dimensionless parameter ε = 2f 2 /w 2 0 , where w 0 is the fundemantal width of the Gaussian.For any fixed ε, IG beams form a complete orthonormal set in terms of their total order N , parity p (= e for even modes and o for odd modes), and index µ, which runs from 0 (1) to N in steps of two for p = e (o).At the waist plane z = 0 these beams are defined (to within an unimportant LG LG-like Boundary HG-like HG N,µ (ν, ε)e −r 2 /w0 2 . ( The variation of an IG beam profile with ε (for fixed values of the beam's indices) is shown in Fig. 1.Note that for ε → 0 the beam reduces to a real LG beam with (N − µ)/2 radial nodes and 2µ azimuthal nodes (that is, µ azimuthal field oscillations).At the other limit, for ε → ∞ the beam reduces to a HG beam with µ and N − µ nodes in each of the Cartesian directions.As ε varies between these limits the profile undergoes a transition in which two regimes can be recognized: the LG-like regime where the beam resembles a deformed LG beam since it is confined between two concentric loops, and the HG-like regime where the beam resembles a deformed HG beam since it is confined within a quadrangle with curved sides.The distinction between these two regimes becomes more marked for modes with higher N .Note that, in all cases, the beams carry no net intrinsic orbital angular momentum (OAM), and hence the limit ε → 0 corresponds to real LG beams and not ones presenting a phase vortex.
The separable Gaussian solutions discussed so far are part of a larger set of paraxial beams known as structured Gaussian (SG) beams [17,18,[29][30][31][32][33].SG beams preserve their intensity profile under propagation (up to a hyperbolic scaling), and are the modes of stable resonant cavities with spherical mirrors [2,21,34].These cavities happen to be mathematically analogous to a two-dimensional rotationally symmetric quantum harmonic oscillator, and thus the functional form of SG beams at their waist plane also describes the eigenmodes of this oscillator.The quantization of the cavity modes is analogous to the energy quantization of the oscillator modes, and given the systems' symmetry they both depend on only one index: the total order N .All modes with equal N are then degenerate, explaining why a range of mode shapes (including HG, LG, IG for any ε and many others, all with any orientation) are possible.The degenerate subspace of modes with given N can be mapped onto a spin system with SU(2) symmetry by using Schwinger's coupled oscillator model [35].This connection has allowed the use of various spherical representations to describe SG beams and their changes under specific modal transformations, as well as the phases they can accumulate [17,32,36,37].
The degeneracy of the mirror cavity can be broken by inducing perturbations, such as small amounts of well-chosen aberrations.It was shown recently that IG modes with a given ε can be selected by adding small amounts of astigmatism and spherical aberration [24], their ratio being proportional to ε.In particular, a small amount of only spherical aberration selects LG modes, while a small amount of astigmatism alone leads to HG modes aligned with the symmetry axes of the aberration.Furthermore, this aberrated cavity system was shown to be mathematically analogous to a Bose-Hubbard dimer, a widely used model in condensed matter physics [24,38].
In this work, we take a fresh look at IG beams from a ray-optical perspective.Rays, and their envelopes known as caustics, have been used to explain intuitively the behavior of wave fields [17,30,33,[39][40][41][42][43].When applied to SG beams, the ray model can explain both their "self-healing" behavior, as well as the Gouy phase accumulated around the focal region and the geometric phase accumulated during a sequence of mode transformations [17,30].Rays can also be used to design new SG beams [17,44], and to describe even SG beams that can be deemed as "the least ray-like" [33].Here we show that a raybased description provides important insights into the shape of IG beams, revealing a surprising amount of geometry that connects them to several other physical phenomena.In particular, this ray treatment is naturally compatible with the cavity aberrations mentioned earlier, and hence clarifies the two regimes (HG-like and LG-like) from the shape of the resulting caustics, as well as the boundary between these regimes, which can be understood as a ray-optical topological transition.We also show that the ray description is essentially complete, since it is sufficient to calculate accurate wave field estimates.

II. RAY FORMALISM FOR SG BEAMS
In the paraxial ray formalism, at a given propagation distance z, say z = 0, each ray can be identified by its transverse coordinates Q = (Q x , Q y ) and its transverse direction cosine vector P = (P x , P y ).These quantities can be used to constitute a four-dimensional phase space (Q x , Q y , P x , P y ), where each point represents a possible ray.An optical beam is represented by a two-parameter family of rays that defines a surface (referred to as the Lagrange manifold) within this phase space [41,43].Typically, for well-localized beams, this manifold is topologically a torus, so that the two parameters, say τ and η, are periodic.For SG beams, the requirement of self similarity means that the cross-section of the ray family must maintain its overall shape up to a global scaling, and as shown elsewhere [17,18] this imposes a specific dependence for Q and P on one of the two parameters, say τ , so that at z = 0 they satisfy where v is a two-dimensional complex vector with unit norm.The normalization of v implies that This expression is equivalent to the equation for the conservation of energy of a two-dimensional harmonic oscillator, and it restricts the region of the 4D phase space occupied by the Lagrange manifold.Note that, for fixed η, both Q and P describe ellipses as τ varies, just like the time evolution for the position and momentum of a harmonic oscillator.Upon propagation, each of these elliptic families of rays (EFR) describes a ruled hyperboloid with an elliptic cross-section that satisfies the property of self similarity, as shown in Fig. 2: the elliptical cross-sections get larger upon propagation, but they retain their ellipticity and orientation.
An EFR is fully determined by the normalized complex two-vector v. Since the vector is normalized and its global phase is irrelevant, it can be parametrized as with θ and ϕ being the colatitude and longitude angles of a sphere, respectively.This vector is then mathematically equivalent to the Jones vectors used to represent polarization states as points on the surface of a Poincaré sphere.It is then possible to represent any EFR as a point on the surface of ray-Poincaré sphere as shown in Fig. 2. Like for polarization, the angle θ determines the ellipticity and handedness of the EFR, while ϕ determines its orientation.These angles are functions of the second parameter, η, so that an SG beam is described by a path over the surface of the sphere, referred to here as the Poincaré path (PP).Because this second parameter is typically periodic, the PP is typically a closed path.
The analogy between the ray structure for SG beams and polarization can be extended by introducing analogs of the Stokes parameters, which define a Cartesian coordinate system for the ray-Poincaré sphere.These are the Fradkin-Stokes parameters, defined as [24,29] For each ray, these parameters are invariant under propagation in free space, with T 3 being proportional to the skew invariant, which is a measure of the OAM of each ray.Further, it can be shown that these parameters are constant also over all rays within a given EFR as defined by the expression in Eq. ( 2) for fixed η and varying τ .In fact, an alternative definition for the EFR is as the set of all rays with given values for all three T n .It can be shown that Eq. ( 3) can be written as so these parameters are indeed constrained to a spherical surface in the space (T 1 , T 2 , T 3 ).Because the Poincaré sphere is often defined as a unit sphere, it is convenient to use the normalized quantities t j = 2T j /(N + 1), so that The ray families of the most common SG beams are simply related to the parameters T j (and hence t j ).For example, LG beams have definite OAM, and hence their ray-optical representation is composed of only rays with given t 3 .Their PP then corresponds to the intersection of the sphere with a plane of constant t 3 , that is, to a horizontal circle over the sphere whose height determines the OAM as shown in Fig. 2. A similar relation holds for HG beams oriented along the x and y axes, which are composed instead of rays with given t 1 , and whose PP are circles normal to the t 1 axis (see Fig. 2).Ray families for which a linear combination of t 1 and t 2 is constant correspond to HG beams with different orientations.More generally, ray structures for which a linear combination of all three parameters is constant, that is, where the PP is a circle with any orientation, correspond to the so-called generalized Hermite-Laguerre-Gauss (HLG) beams [18,29,31,45], which include the HG and LG beams as special cases.Note that, except for these two special cases, HLG beams are not separable in any coordinate system.

III. RAY STRUCTURE FOR THE IG BEAMS
A. Poincaré path as a generalized Vivani curve IG laser modes are known to result from aberrations or misalignment in a cavity [24,34].In fact, before IG beams were proposed, modes in aberrated cavities resembling IG modes were observed experimentally and studied theoretically in terms of ellipsoidal coordinates (beyond the paraxial regime) [46].We showed recently [24] that T 1 can be associated with the effect of a small amount of astigmatism in the cavity mirrors, whose effect is precisely to break the symmetry in the x and y directions, while a small amount of spherical aberration (which is nonlinear in Q and P) preserves rotational symmetry and can be associated instead with T 2 3 , where the square is consistent with the fact that the aberration does not have an intrinsic skewness.
As shown in Ref. [24], IG modes result when both aberrations are present.These modes are therefore associated with ray families for which a linear combination of T 2 3 and T 1 equals some constant, where ε determines the ratio between these aberrations: That is, unlike HLG beams, these modes involve quadratic combinations of the Stokes-Fradkin parameters.This relation can be written in terms of the normalized parameters as This equation defines a parabolic cylinder parallel to the t 2 axis, whose intersection with the unit sphere corresponds to the PP.However, the resulting curve is easier to understand by using Eq. ( 7) to eliminate t 3 from Eq. ( 9), giving which is the equation of a circular cylinder of radius R, with an axis parallel to the t 3 axis and that crosses the t 1 axis at c, with , and The PP for IG modes then corresponds to the intersection of the ray-Poincaré sphere with this cylinder, as shown in Fig. 3.Note that the condition −1 < c − R < 1 must be satisfied for an intersection between the cylinder and the sphere to exist.For the specific case of c = R = 1/2, this intersection corresponds to the Viviani curve (or Viviani window), proposed by seventeenth century mathematician Vincenzo Viviani [47].More generally, intersections of a cylinder with a sphere are referred to as generalized Viviani curves, or also as Euclidic spherical ellipses since the sum of the distances to two focal points on the surface of the sphere equals a constant [48].Generalized Viviani curves appear naturally in the description of a range of physical phenomena, including the Bose-Hubbard dimer [48], which as discussed earlier presents mathematical analogies with the aberrated cavity [24].In optics, generalized Viviani curves describe polarization evolution in nonlinear birefringent fibers (also analogous to the Bose-Hubbard dimer) [49], and the subset of the generalized Viviani curves that are eight-shaped (with c + R = 1) describe the polarization evolution for linearly polarized light traversing a rotating waveplate, or conversely the polarization states explored by a rotating-waveplate polarimeter [50].
We can now see that the regimes for the shape of IG beams in Fig. 3 correspond to the different types of generalized Viviani curves for the PP: 1.For c + R < 1 (that is, for ε < a/(N + 1)) the PP consists of two disjointed loops, and this defines the LG-like regime, since the PP resembles a deformed version of that for a real LG beam (two horizontal circles that are mirror images of each other with respect to the t 3 = 0 plane).
2. For c + R = 1 (that is, for ε = a/(N + 1)) the PP is an eight-shaped curve, which defines the boundary between the two regimes.
3. For c + R > 1 (that is, for ε > a/(N + 1)) the PP consists of a single loop, and this defines the HG-like regime since the PP resembles a deformed version of that for HG beams (a single vertical circle).
That is, in this ray-optical description, the transition between the LG-like and the HG-like regimes for IG beams corresponds to a topological transition for the ray family: in the LG-regime, the ray family is actually composed of two separate ray families (two phase-space tori) that are mirror images of each other.As ε increases these two families deform until they share an EFR for ε = a/(N + 1), and they then merge into a single continuous ray family (one phase-space torus) in the HG-like regime.Note that while the boundary is sharp in the ray regime, it becomes blurry when we take into account the wave nature of light.That is, for values of ε near the transition there can be appreciable tunneling effects between the two separate parts of the PP that approach each other.The shape correspondence between the ray and wave pictures is shown in Fig. 3.

B. Caustic structure in terms of Keplerian and harmonic oscillator trajectories
As we can see from the bottom row in Fig. 3, all rays are restricted to a region of space confined by an outer surface, and also by an inner surface for the case of LG-like beams.These surfaces correspond to what is known as caustics, which are the envelopes of the rays.Due to the larger density of rays in their vicinity (on one side), caustics are associated with the brightest regions of an optical field [39][40][41].Since the PP encodes the two-parameter ray family, the caustics can be inferred from it.The procedure for doing so turns out to be of a geometrical nature [17], particularly for the case of IG beams considered here.The coordinates (t 1 , t 2 ) of the PP are sufficient for the determination of the caustics, since the only extra information provided by t 3 is the sign of the rays' skewness, which does not affect where they form an envelope.As shown in the top row of Fig. 3, the ray-Poincaré sphere can then be projected onto a unit disk, referred to here as the equatorial Poincaré disk (EPD).Since for IG beams the PP is the intersection of a vertical circular cylinder with the sphere, the PP projects onto a circle (or a circular segment) of radius R centered at (t 1 , t 2 ) = (c, 0).This circle can be parametrized as where φ is the angle measured from the positive t 1 axis (see Figs. 4 and 5).In the LG-like regime φ traces a complete cycle while in the HG-like regime it is constrained to [φ 0 , 2π − φ 0 ] with We now describe the geometric procedure for finding the caustics from the PP projection.As discussed in Ref. [17], the first step is to find what is known as the medial axes or topological skeleton [51] between the PP projection in Eq. ( 12) and the unit circle (the edge of the EPD).Medial axes are curves that are equidistant to two other curves, and are defined as the loci of centers of what we call here bitangent circles, which are the circles that are tangent to the two original curves.This is shown in Fig. 4 for the case of interest here, in which the two original curves are the PP projection and the unit circle.For each point of the PP projection (corresponding to a given value of φ), there are two bitangent circles, one of radius r o and whose center m o is on the opposite side of the PP projection as its center, (c, 0), and one of radius r s and whose center m s is on the same side of the curve as (c, 0).Because the PP projection is also a circle, the distance between its center (c, 0) and each of the two medial points is just the sum or difference of R and the radius of the corresponding bitangent circle, namely Similarly, because both bitangent circles touch the unit circle, the distance between the origin and the medial points is the difference between unity (the radius of the unit circle) and the corresponding radius, that is, The equations for the medial axes are obtained by eliminating the radii r o and r s from these relations, and we now consider each medial axis separately.
For the medial axis traced by m o the elimination of r o leads in all cases to the equation which describes an ellipse with major axis 1+R and whose foci are the origin and (c, 0).This construction is shown in the top row of Fig. 4 and on the left side of Media 1-4.Note that in the LG-like regime the medial axis traces the complete ellipse, and the geometry of this construction is related of what is

Same side
Opposite side LG-like Boundary HG-like 4. Construction of the medial axes from the PP projection for the two regimes and their boundary.(top) Medial axis (thick blue curve) traced for varying φ by the point mo which is the center of the bitangent circle (blue) of radius ro that touches both the edge of the unit disk and the PP projection (thick black circle with radius R) at the point (red dot) corresponding to a given angle φ.Note that the PP projection's center and mo are on opposite sides of the tangent point to the PP projection.(bottom) Same construction (with orange instead of blue) for the medial point ms, which is on the same side of the tangent point of the PP projection as this curve's center.A more detailed view of this construction is shown in Media 1-4 for the different regimes.
known as a Steiner chain [52], a remarkable geometric result by 19th-century Swiss mathematician Jakob Steiner.The ellipse is also traced completely in the boundary case, for which the construction resembles instead a Pappus chain [52], described by Pappus of Alexandria in the third century.In the HG-like regime, on the other hand, only a fraction of the ellipse is traced by m o , given the fact that the PP projection is an incomplete circle.
The equation for the second medial axis, traced by m s , requires a different treatment for the two regimes, as well as for their boundary.In the LG-like regime we have that r s > R and R < 1, so that the right-hand side of the second expression in Eqs.(14) equals r s − R. The elimination of r s from Eqs. ( 14) and ( 15) then gives which also describes a complete ellipse with foci at the origin and (c, 0) but with major axis 1 − R, as shown in the bottom-left panel of Fig. 4) and the right panel of Media 1.In the HG-like regime, on the other hand, because R > r s the right-hand side of the second expression in Eqs. ( 14) equals R − r s , and the elimination of r s gives instead which describes a branch of a hyperbola with the same foci, as illustrated in the bottom-right panel of Fig. 4) and the right panel of Media 3,4.This hyperbolic branch presents different behavior for different values of R: it opens away from the origin for R < 1 (as in the case shown in the bottom-right part of Fig. 4 and in Media 3), it is a straight line for R = 1, it opens towards the origin for R > 1 (as shown in in Media 4), and it becomes a parabola with focus at the origin in the limit R → ∞ (corresponding to the HG case).In fact, in this latter limit, the ellipse traced by m o also becomes a parabola.At the boundary of the two regimes, for which c+R = 1, the medial axis traced by m s undergoes a transition that makes it occupy the complete unit radial line corresponding to the positive t 1 axis, as shown in the bottom-center panel of Fig. 4, because any circle centered on this line segment that is tangent to the unit disc at t 1 = 1 is also tangent to the PP projection and hence is a valid bitangent circle.Note, though, that two points along this straight medial axis are special, namely the origin and the point at (c, 0).When m s coincides with one of these points, r s equals either unity or R, meaning that the bitangent circle is identical either to the edge of the EPD or the PP projection and hence is tangent to the corresponding curve not only at one point but all along its perimeter.This transition case is illustrated in the right-hand panel of Media 3.
In summary, for both regimes the medial axes are conic sections with one focus at the origin of the EPD and the other at (c, 0).Therefore, geometrically, their shape coincides with that of Keplerian orbits for a potential centered at the origin.These orbits are closed (elliptic) and correspond to an attractive force for m o in all cases (although m o does not trace the complete orbit for c + R > 1), and also for m s in the LG-like regime.For m s in the HG-like regime, the force is attractive if R > 1 (although the orbit is not periodic), repulsive for R < 1, and there is no force for R = 1.In the limit c → ∞ (and hence R → ∞) corresponding to HG beams, both m o and m s trace parabolas corresponding to escape orbits, while in the opposite limit of c → 0 corresponding to LG beams, both m o and m s trace circular orbits.
The medial axes can now be mapped directly onto the caustics by using a square-root conformal map [17].This is done by associating a medial point m = (t 1 , t 2 ) with the complex number t 1 + it 2 and then taking its square root: The resulting point y ) corresponds to a caustic point.Given the sign ambiguity of the square root, each medial axis point maps onto two caustic points, and therefore each of the two medial axes maps onto two segments of the caustic curves.Remarkably, the square-root conformal map is known to have the property of mapping Keplerian orbits onto the orbits of a particle in a (repulsive or attractive) two-dimensional harmonic oscillator, which correspond to conic sections centered at the origin [53].(Incidentally, the properties of this mapping when applied to conics make it also useful for studying the shape of light reflectors in non-imaging applications [54].)Since the medial axes of IG beams correspond precisely to Keplerian orbits, their caustics take the shape of confocal conic sections with foci that coincide with those for the elliptical coordinate system, (±f, 0).For beams in the LG-like regime, the caustics are two complete confocal ellipses, which in the limit of LG beams become two concentric circles.In the HG-like regime, on the other hand, the caustic structure is composed of two segments of an ellipse and two segments of a hyperbola, forming a deformed quadrangle, and in the limit of HG beams these segments become straight lines (since the square-root map transforms parabolas to straight lines).In the boundary case, the outer caustic is an ellipse and the inner one occupies the major axis of this ellipse, with two main points that coincide with the foci.In all cases, all rays are contained within these caustics.

IV. SEMICLASSICAL ESTIMATES A. Wave reconstruction
The wave behavior of an optical field can be reproduced faithfully from the ray description if an appropriate semiclassical method is used.Since there is a closed-form expression for the IG beams, it is not necessary to rely on a semiclassical approximation to calculate them.However, it is still interesting to show that the ray picture is essentially sufficient to reconstruct these fields.Further, this asymptotic analysis helps illustrate some interesting aspects of the topological transition between the LG-like and the HG-like regimes, and helps elucidate a connection with yet another physical system.
The semiclassical approach we use is a Gaussian summation method [42,55,56], since it does not diverge at the caustics and was shown to reproduce the exact paraxial solutions for HG, LG and HLG beams [17,18,57,58].In this method, each ray is dressed with a Gaussian contribution, and the estimate takes the form of an integral over τ and η.When applied to SG beams, the integral in τ over a period can be solved analytically, and the constant factors in Eq. ( 2) guarantee that the integrand is consistent at the endpoints of integration.The semiclassical estimate then takes the form of an integral of contributions associated with each EFR that allows estimating a SG beam from its ray structure [17,18]: where A(η) is a weight factor for the contributions from each EFR, primes denote derivatives with respect to η, U N is a mathematical analog of a spin coherent state [32,[59][60][61] defined as with H N representing an N th order Hermite polynomial, and gives a phase to each contribution that can be interpreted as a geometric phase, since it corresponds to half an area on the ray-Poincaré sphere.This integral can be performed analytically in terms of elliptic integral functions.(See appendix A for more information.) While the periodicity of the integrand in τ is guaranteed by the factors in Eq. ( 2), the periodicity of the integrand in η imposes conditions on the phase contribution due to the geometric phase caused by the PP enclosing a solid angle Ω over the ray-Poincaré sphere.Because the square-root factor in Eq. ( 20) accumulates a phase of ±π on closing the PP loop, the solid angle must then be an odd multiple of 2π/(N + 1), [17,18] for integer s.This condition restricts the valid values for R in Eq. ( 10) needed to define the PP.However, recall that the number of closed loops composing the PP differs for the two regimes.For the LG-like regime there are two loops and each must satisfy the equation above, so the total subtended solid angle must be an odd multiple of 4π/(N + 1) (since the two loops enclose the same solid angle).On the other hand, in the HG-like regime the PP consists of a single loop so the solid angle must be an odd multiple of 2π/(N + 1).This leads to a discontinuity at the ray-optical boundary on the allowed values for R for a given c; as shown in Fig. 6 there are N + 1 allowed values of R for the HG-like regime, but only ⌊(N + 1)/2⌋ for the LG-like regime.

B. The Ince equation and its eigenvalues
In the wave domain, IG beams are defined as eigenstates of what we call here the Ince equation, which corresponds to the operator version of Eq. ( 8): where the operators T j are obtained from the expressions for T j in Eqs. ( 5) by substituting P → −ik −1 (∂ x , ∂ y ).Note that a corresponds now to an eigenvalue.The dotted lines in Fig. 6 were calculated by using the values of a found through the numerical solution of Eq. ( 24).These curves are continuous at the boundary, and they merge by pairs as one transitions into the LG-like regime.This is because LG-like IG beams with equal indices but different parity are nearly degenerate, and semiclassically they are represented by the same pair of PP loops, the only difference being the relative phase between the two contributions.The discrepancy between the curves in Fig. 6 for the semiclassical estimate and the result FIG. 7. Mapping between the phase space curves of the pendulum and the PPs of IG beams on the ray-Poincaré sphere.
using the exact eigenvalue is due to the fact that the semiclassical approach neglects tunneling effects that take place near the boundary.A more accurate semiclassical estimate would require considering asymptotic corrections to the amplitude term derived below, which would modify the phase consistency condition that leads to Eq. ( 23).This type of correction was used recently in the asymptotic estimation of eigenvalues for anharmonic potentials [62].Nonetheless, away from the boundary between the two regimes, the curves produced by the semiclassical estimate are in good agreement with those calculated using the exact eigenvalues, even for moderate N .

C. The semiclassical amplitude and a connection with the simple pendulum
The remaining element for evaluating the semiclassical estimate in Eq. ( 20) is to find the amplitude A. The choice of this function, as long as its phase is constant and it does not vary abruptly, does not affect the fact that the resulting beam is self similar or that it has a shape similar to that of IG beams.However, the form for this function that makes the estimate approach closely the IG beams can be calculated from substituting Eq. ( 20) into Eq.( 24), leading to an equation that can be separated into different orders of N (see appendix B for more detials).Unsurprisingly, the dominant term gives Eq. ( 8), namely the PP.The next term in the series leads to an equation for the amplitude A(η), and hence to the appropriate weighting for each EFR.The resulting equation for the weight takes a simple form when written in terms of the Cartesian parameters t j (η): where B is a constant of integration.If we assume that the parametrization is such that A is constant, then t ′ 3 = A 2 t 2 /B.For simplicity we set A 2 /B equal to unity.By combining this equation with the PP and the fact that the parameters t j lie on a unit sphere, we are led to the following system of first order differential equations: These equations are analogous to the two-mode Gross-Pitaevskii equations describing the meanfield evolution of a Bose-Hubbard dimer [38,48].They lead to a particularly simple form when written in terms of the angle φ used to parametrize the cylinder in Eq. ( 12): This is the equation for a simple pendulum where the parameter η plays the role of time.This equation implies that the EFR contributions become more important the closest their point at the PP is to the equator (which correspond to points where the pendulum lingers).The two regimes for the IG beams then correspond to the two topologically different regimes for the pendulum: full rotation corresponds to the LG-like regime and libration corresponds to the HG-like regime.The analogous connection between the pendulum and the Bose-Hubbard dimer was pointed out by Graefe et al. [48].Further, the connection with the generalized Viviani curves becomes more intuitive if we use Eq. ( 12) to write t 3 in terms of and substitute this into the third relation in Eqs.(26), which leads to the simple expression That is, t 3 corresponds to the pendulum's velocity (or momentum, for a unit mass).One can then define a two-dimensional phase space (φ, t 3 ) for the pendulum, in which its motion is described by a curve.However, given the possible periodicity of φ one can roll up this phase space into a cylinder of radius R. The pendulum's phase space curve would then coincide with the PP for the IG beams, that is, with the intersection of the cylinder with a sphere as shown in Fig. 7.
The fact that the rolling-up of the phase space of the pendulum falls on a sphere (or ellipsoid) is all the more surprising because the shape of the intersections of the cylinder with the sphere depends on two parameters, c and R, while for a pendulum one considers typically that the orbits in phase space depend on a single parameter.Let us use as this parameter h as defined in Eq. (13).The libration and rotation regimes for the pendulum then correspond to h being smaller or larger than unity, respectively.It is easy to show that the equation for the PP, when unrolled onto the pendulum's (φ, t 3 ) phase space, becomes t 2 3 = 2cR(h − cos φ), so that the shape of the curve depends only on h, and the factor of 2cR just provides a constant scaling for the pendulum's momentum that, amongst other things, guarantees that |t 3 | ≤ 1.

D. Ray-based estimates
We conclude by testing the ray-based estimate for the IG beams, to show that the ray families obtained recreate faithfully the wavefields.Since the integral in Eq. ( 20) has no closed-form solution, we approximate it as a discrete sum of M terms: where we use the parametrization defined by Eqs. ( 26) and ( 27) that makes A = 1.The dependence of the angles θ and ϕ on the parameter η is determined by solving these differential equations with appropriate initial conditions determined by the PP.Setting η 0 = 0 then η j = jη (p) /M where the pendulum's period is given by , and where K and F are, respectively, the complete and incomplete elliptic integrals of the first kind.
In the LG-like regime, each of the two loops provides a contribution, one being the complex conjugate of the other.However, the semiclassical theory does not constrain how these contributions are combined.Figure 8 shows the field due to the northern PP loop, resulting in a beam with the same elliptical shape as the corresponding IG beam, but with nonzero net OAM.The southern loop would then have the opposite OAM.These are then estimates (when sufficiently far from the boundary) for what has been referred to as helical IG beams, which are obtained by combining even and odd IG beams with the same N and µ as IG [21].As shown in Fig. 3, for both the upper and lower parts, the projections on the EPD of the sampling points used to compute the estimate provide us with snapshots of the evolution of a pendulum in full rotation.Real IG beams can be obtained by adding both contributions either in phase or out of phase (that is, by taking the real or imaginary part of the contribution from one loop).In the HG-like regime, on the other hand, the PP is given by a single loop so that Eq. ( 29) directly provides a ray-based estimate for the corresponding IG beam.Here, the projections on the EPD of the sampling points used to compute the estimate correspond to snapshots at equal temporal intervals of the libration evolution of the pendulum.Note that it is not possible to obtain an estimate at the boundary between the two regimes because the quantization condition becomes ill-defined.
Figure 8 also shows the excellent agreement between the ray-based estimates and the exact fields for IG beams within the two regimes, despite the low values of M being used.Figure 9 shows the rms error in terms of the number of contributions M , for beams with the same indices as in Fig. 8 and for others with a lower value for µ for different values of c (and hence of ε).The plots show that in all cases the sums converge rapidly (typically for M < N ), reaching a level of error that is due to the semiclassical approach itself, which is more accurate as the difference between c and the boundary between regimes increases.This means that this approach might be computationally convenient for large N .

V. CONCLUDING REMARKS
In this work we explored the structure of IG beams from the ray-optics perspective.This ray structure embodies the main features of these beams, in particular the two separate regimes for their overall shape, as well as the boundary between them, which corresponds to a ray-optical topological transition.The analysis also reveals connections with many geometric constructions such as the Steiner and Pappus chains, Viviani curves, Keplerian orbits and medial axes, as well as physical systems like the harmonic oscillator, the pendulum and the Bose-Hubbard dimer.
As shown in Ref. [24] and discussed in this work, IG beams are modes of resonant cavities whose curved mirrors contain a small amount of spherical aberration and astigmatism.Let us discuss how this fact sheds light on the structure of the beams from the ray-optical perspective.Imagine a subset of rays corresponding to a single EFR like those shown in Fig. 2. When confined inside an unaberrated cavity with appropriate mirror separation and curvatures, this group of rays retains its configuration after a round trip.Note though that each ray within the EFR changes position after each round trip, cycling within the EFR and exploring different values of the parameter τ .In fact, except for very specific cases in which certain rational relations hold, the complete EFR can be populated by letting a single ray bounce back and forth.Now consider the effect of small amounts of astigmatism and spherical aberration in the cavity.These aberrations cause the EFR to very slightly modify its profile after each round trip.If the aberrations are in the perturbative regime, the deformed EFR will still have an elliptic profile and therefore will still constitute an EFR, one corresponding to a nearby point over the ray-Poincaré sphere.Multiple round trips then correspond to a slow displacement of the point over the ray-Poincaré sphere, which traces the PP (a generalized Viviani curve if the aberrations are astigmatism and spherical).Remarkably, the temporal rate of this slow evolution corresponds, when projected onto the (t 1 , t 2 ) plane, to the rotation or libration of a simple pendulum.That is, a ray bouncing back and forth in the aberrated cavity presents a fast cyclic evolution in τ and a slow pendulum-like precession in η.From this point of view, we see that the connection between IG modes and the pendulum does not arise in the semiclassic regime (that is, for the construction of wave modes) but in the ray-optical one.
The ray-optical picture of IG beams provides intuition not only on their behavior within a cavity but also under free propagation.For example, rays give a simple description of the transport of energy (and therefore information) over significant distances, especially for beams with sufficiently large N for with the transverse profile is more structured.If one were to block a sector of a beam's profile, the change in this profile some distance away from the obstacle can be understood intuitively in terms of the rays that were blocked.Amongst other things, this explains why IG beams, like most other structured self-similar beams, present the property of "self-healing", in which the intensity features that were blocked reappear after some propagation distance (at the cost of others disappearing).This effect is due to the cycling under propagation of the rays within each EFR: the blocked rays associated with the missing feature are replaced by non-blocked rays, making the feature reappear.Finally, note that the ray representation of a structured beam allows the use of fast ray-based computations of the trapping forces and torques resulting from the interaction of the beam with an object whose spatial features are larger than the wavelength.
Finally, let us remark that the description given here can be used to predict the (ray and wave) shape of self-similar beams that are modes of cavities presenting other types of aberration.By associating such aberrations with powers of the parameters t j , one can easily find the shape of the PP, and from it that for the caustics by using the medial axis construction.Nevertheless, it is unlikely that any other set of aberrations would lead to beams with such a simple yet rich geometric description as the one shown here.For IG beams it is remarkable that, despite the nonlinear (that is, nonplanar) nature of the PP and its resulting different topological regimes, the caustics end up being simple confocal conic sections.functions: where }, E is the elliptic integral of the second kind, F elliptic integral of the first kind, and Π is the complete integral of the third kind.However, there is a caveat for the use of this expression: it is defined as a continuous function for ϕ within the interval [−π, π] such that it is zero for ϕ = 0 in the LG-like regime and for ϕ = ±φ 0 in the HG-like regime.
There is therefore a discontinuity at ϕ = ±π.This discontinuity poses no problem in the LG-like regime since the quantization of the solid angle ensures that, when placed in the exponent of the semiclassical estimate, the resulting phase discontinuity is an integer multiple of 2π.For the HG-like regime, on the other hand, the discontinuity must be removed by hand to make the function continuous at ϕ = ±π.
Unfortunately the limit ϕ → ±π in this analytic form is complicated to evaluate, so the constant that must be added/subtracted must be evaluated numerically by using a value very near this limit.

Appendix B: Derivation of the amplitude equation
To compute the ray-based estimate, we employ a Gaussian summation method where each ray is dressed with a Gaussian contribution of the form g(r; Q, P) = kΓ 2π exp − kΓ 2 ∥r − Q∥ 2 + ik(r − Q) • P e ikL(η,τ ) .(B1) The center Q and direction P of the Gaussians are the same as those of the corresponding rays.Γ is a free parameter that can be chosen to optimize the field estimate.For the study of SG beams, its optimal value is the inverse Rayleigh range, Γ = 2/kw 2 0 .Given this choice, these Gaussian contributions correspond to displaced coherent states of the two-dimensional isotropic harmonic oscillator.As opposed to what is done in the main text, here we keep the dependence on both parameters of the ray family so that the estimate is given by U (r) = A(η) J(τ, η) g(r; Q, P)dτ dη (B2) where L(η, τ ) is the optical path length that plays the role of the action in mechanics.The weight A only depends on the parameter η so that all the rays in the same ellipse are weighted equally, which is a necessary condition to obtain a self-similar beam.
To determine the appropriate equation for the amplitude, we substitute the semiclassical estimate into the Ince equation, T 2 3 g(r, η, τ ) = (N + 1) 2 (p x q y − p y q x ) 2 + N + 1 2 q • z − 2i(N + 1) 3/2 w (p x q y − p y q x )(z x ∆ y − z y ∆ x ) Using the previously derived identities, all the spatial dependence can be passed onto g so that it can be factored out along with A and j 1/2 .The rest of the integrand can then be set equal to zero.Moreover, since the classical (or ray) limit is attained for large N , the contributions of different orders can be treated separately.The first order leads to the following equation, (N + 1) 2 (p x q y − p y q x ) 2 + ε 2 Noting that T 2 3 = (N + 1) 2 (p x q y − p y q x ) 2 it is clear that the previous equation is nothing more than the PP path.The following term requires more work, but after some tedious mathematical manipulations we get 0 = (p x q y − p y q x )(z 2 x + z 2 y ) ′ A √ j + 2 (p x q y − p y q x )(z 2 x + z 2 y ) This equation says how the weight function changes depending on the parametrization of the PP.This equation can be integrated after we multiply it by A/j 1/2 , thus leading to A 2 = Bj (p x q y − p y q x )(z 2 x + z 2 y ) − εi 2(N +1) z x z y , with B being a constant of integration, and for the remainder of the derivation we will set B = i without loss of generality.Both the equation for the PP and that for the amplitude A depend only on η, as can be appreciated when writing them in terms of θ and ϕ, (N + 1) 2 cos 2 θ + ε(N + 1) cos ϕ sin θ = a (B10a) and substituting it into the second one we get Finally, by writing this equation in terms of φ, we have that where we used the relation c = ε/2(N + 1).The denominator is proportional to the angular velocity of a simple pendulum, therefore if we set φ = η then the amplitude of the ray-ellipses is inversely proportional to the square root of the angular velocity of the pendulum.Another way to look at it is to assume that the amplitude is constant, say A = 1, then the solution of the angular variable φ follows the dynamics of the pendulum where η plays the role of time.The connection to the simple pendulum becomes clearer after taking the second derivative of the previous equation, leading to φ ′′ − cR sin φ = 0. (B14) This is the defining equation of a simple pendulum.

FIG. 1 .
FIG. 1. Intensity distribution for the even IG beam with total order N = 22 and µ = 18 as ε varies from 0 to ∞.At the extremal values ε = 0 and ∞, the IG beams reduce to real LG and HG beams, respectively.At intermediate values they take the form of either deformed LG-like or HG-like beams when ε is below or above a given boundary.

FIG. 2 .
FIG. 2. (left) Ray-Poincaré sphere representation for elliptic families of rays (EFRs).Several EFRs are shown corresponding to different points over the ray-Poincaré sphere.(right) Examples of Poincaré paths for vortex LG and HG beams.

FIG. 3 .
FIG. 3. Ray structure of an IG beam (N = 22, even and µ = 18) as it transitions between the LG (c = 0) and HG (c → ∞) limits.(top) The intersection of the RPS with the surfaces defining the Poincaré path (PP) as their intersection along with the projection of the PP on the EPD where the angle φ is defined.(middle) Elliptic families of rays color coded as the PP on the first row overlaid on top of the amplitude distribution for the corresponding beam.(bottom) Three-dimensional propagation of the families of rays.

FIG. 5 .
FIG. 5. Construction of the caustics from the PP as they transition between the LG (c = 0) and HG (c → ∞) limits.(top) Equatorial Poincaré disk (EPD) with the projection of the PP (inner black curve) and the two medial axes (blue and orange curves) with examples of the bitangent circles used to construct them, (middle) the caustics obtained via the conformal mapping of the medial axes and (bottom) the caustics overlaid on top of amplitude distribution of the corresponding beam.

FIG. 6 .
FIG. 6. Semiclassical and exact values for the radius R for the IG beams with N = 6, N = 11 and N = 22 as a function of the center coordinate c.The different curves or pairs of curves correspond to different values of µ.

FIG. 8 .
FIG. 8. Examples of sampling and comparison between the exact mode and its ray-based estimate for (first row) helical IG beams in the LG-like regime with (second and third columns) N = 6, µ = 4, c = 0.1 and (fourth and fifth columns) N = 22, µ = 18, c = 0.3, and (second row) even IG beams in the HG-like regime with (second and third columns) N = 6, µ = 4, c = 0.4 and (fourth and fifth columns) N = 22, µ = 18, c = 0.7.The number of terms M was chosen as the smallest for which the estimate is visually indistinguishable from the exact distribution.For the helical IG beams the plots show the complex field with the phase encoded as hue, while for the HG-like IG beams the plots show the amplitude distribution.The figures in the first column show an example of sampling for eight points at equal integrals in η.Note that their projection onto the (t1, t2) plane is equivalent to snapshots at uniform time intervals of the evolution of a pendulum.

FIG. 9 .
FIG.9.Rms error ϵrms between the exact IG beams and the ray-based estimates as a function of the number of terms M used for the estimate for various values of the center c for the LG and HG-like regime and for four different IG beams.