<figure><img src="image_1.png" /></figure>

HAL
open science

# Ray and caustic structure of Ince-Gauss beams

R Gutiérrez-Cuevas, M R Dennis, Miguel Alonso

► To cite this version:

R Gutiérrez-Cuevas, M R Dennis, Miguel Alonso. Ray and caustic structure of Ince-Gauss beams. New Journal of Physics, 2024, 26, 10.1088/1367-2630/ad17dc. hal-04440282

HAL Id: hal-04440282

https://hal.science/hal-04440282v1

Submitted on 5 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Ray and caustic structure of Ince-Gauss beams

To cite this article: R Gutiérrez-Cuevas et al 2024 *New J. Phys.* **26** 013011

View the [article online](#) for updates and enhancements.

You may also like

- Proteolytic remodeling of 3D bioprinted tumor microenvironments

Fatemeh Rasti Boroojeni, Sajjad Naeimipour, Philip Lifwergren et al.

- Preparation and application of low temperature protection materials by Na₂SO₄·10H₂O/PS phase change microcapsules

wenhan Xie, Yongqing Jiang, Zhengmei Liu et al.

- Neural Network Field Theories: Non-Gaussianity, Actions, and Locality

Mehmet Demirtas, James Halverson, Anindita Maiti et al.
<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

12 October 2023

REVISED

12 December 2023

ACCEPTED FOR PUBLICATION

21 December 2023

PUBLISHED

8 January 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_4.png" /></figure>

PAPER

# Ray and caustic structure of Ince-Gauss beams

R Gutiérrez-Cuevas¹,²,*, M R Dennis³,⁴ and M A Alonso²,⁵,⁶

¹ Institut Langevin, ESPCI Paris, Université PSL, CNRS, 75005 Paris, France

² Aix Marseille Univ, CNRS, Centrale Med, Institut Fresnel, UMR 7249, 13397 Marseille Cedex 20, France

³ School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, United Kingdom

⁴ H H Wills Physics Laboratory, University of Bristol, Bristol BS8 1TL, United Kingdom

⁵ The Institute of Optics, University of Rochester, Rochester, NY 14627, United States of America

⁶ Laboratory for Laser Energetics, University of Rochester, Rochester, NY 14627, United States of America

* Author to whom any correspondence should be addressed.

E-mail: rodrigo.gutierrez-cuevas@espci.fr

Keywords: structured light, Ince-Gauss beams, caustics, Bose-Hubbard model, ray optics, pendulum

Supplementary material for this article is available online

## Abstract

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.

## 1. Introduction

The paraxial wave equation accepts many types of solutions, some of which are separable in a given coordinate system [1–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–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–18].

The Ince-Gauss (IG) beams [9, 19–24], separable in elliptic coordinates, have received significant attention in the last couple of decades, finding also many applications [25–28]. Elliptical coordinates require specifying the separation $2f$ between its two foci. Let the radial $\xi \in [0, \infty)$ and angular $\nu \in [0, 2\pi)$ elliptic variables be defined as $\mathbf{r} = (x, y) = (f\cosh\xi \cos\nu, f\sinh\xi \sin\nu)$. Within a region of some length scale $w_0$ around the origin, these elliptic coordinates tend to the Cartesian coordinates as $f \to \infty$, while for $f \to 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 $\epsilon = 2f^2/w_0^2$, where $w_0$ is the fundemantal width of the Gaussian. For any fixed $\epsilon$, 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 $\mu$, 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 normalization factor) in terms of the Ince polynomials $K_{N,\mu}^{(p)}$ as

$$IG_{N,\mu}^{(p)}(\mathbf{r}; \epsilon) \propto K_{N,\mu}^{(p)}(i\xi, \epsilon) K_{N,\mu}^{(p)}(\nu, \epsilon) e^{-r^2/w_0^2}. \quad (1)$$
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Intensity distribution for the even IG beam with total order $N = 22$ and $\mu = 18$ as $\varepsilon$ varies from 0 to $\infty$. At the extremal values $\varepsilon = 0$ and $\infty$, 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 $\varepsilon$ is below or above a given boundary.</figcaption></figure>

The variation of an IG beam profile with $\varepsilon$ (for fixed values of the beam's indices) is shown in figure 1. Note that for $\varepsilon \to 0$ the beam reduces to a real LG beam with $(N - \mu)/2$ radial nodes and $2\mu$ azimuthal nodes (that is, $\mu$ azimuthal field oscillations). At the other limit, for $\varepsilon \to \infty$ the beam reduces to a HG beam with $\mu$ and $N - \mu$ nodes in each of the Cartesian directions. As $\varepsilon$ 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 $\varepsilon \to 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–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 $\varepsilon$ 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 $\varepsilon$ can be selected by adding small amounts of astigmatism and spherical aberration [24], their ratio being proportional to $\varepsilon$. 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–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 ray-based 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.
## 2. 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 $\tau$ and $\eta$, 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 $\tau$, so that at $z = 0$ they satisfy

$$ \frac{1}{w_0} \mathbf{Q}(\eta, \tau) + i \frac{k w_0}{2} \mathbf{P}(\eta, \tau) = \sqrt{N+1} \mathbf{v}(\eta) e^{-i\tau}, \qquad (2) $$

where $\mathbf{v}$ is a two-dimensional complex vector with unit norm. The normalization of $\mathbf{v}$ implies that

$$ \frac{1}{w_0^2} \|\mathbf{Q}\|^2 + \frac{k^2 w_0^2}{4} \|\mathbf{P}\|^2 = N+1. \qquad (3) $$

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 $\eta$, both $Q$ and $P$ describe ellipses as $\tau$ 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 figure 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 $\mathbf{v}$. Since the vector is normalized and its global phase is irrelevant, it can be parametrized as

$$ \mathbf{v}(\eta) = \cos[\theta(\eta)/2] e^{-i\frac{\phi(\eta)}{2}} \epsilon_+ + \sin[\theta(\eta)/2] e^{i\frac{\phi(\eta)}{2}} \epsilon_-, \qquad \epsilon_{\pm} = (\hat{\mathbf{x}} \pm i\hat{\mathbf{y}})/2^{1/2}, \qquad (4) $$

with $\theta$ and $\phi$ 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 figure 2. Like for polarization, the angle $\theta$ determines the ellipticity and handedness of the EFR, while $\phi$ determines its orientation. These angles are functions of the second parameter, $\eta$, 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]

$$ T_1 = \frac{1}{2w_0^2} (Q_x^2 - Q_y^2) + \frac{k^2 w_0^2}{8} (P_x^2 - P_y^2), \qquad (5a) $$

$$ T_2 = \frac{1}{w_0^2} Q_x Q_y + \frac{k^2 w_0^2}{4} P_x P_y, \qquad (5b) $$

$$ T_3 = \frac{k}{2} (Q_x P_y - Q_y P_x). \qquad (5c) $$

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 equation (2) for fixed $\eta$ and varying $\tau$. 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 equation (3) can be written as

$$ T_1^2 + T_2^2 + T_3^2 = \left( \frac{N+1}{2} \right)^2, \qquad (6) $$

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

$$ t_1^2 + t_2^2 + t_3^2 = 1. \qquad (7) $$
<figure><img src="image_6.png" /><figcaption>Figure 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.</figcaption></figure>

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 figure 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 figure 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.

### 3. Ray structure for the IG beams

#### 3.1. 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_3^2$, where the square is consistent with the fact that the aberration does not have an intrinsic skewness.

As shown in [24], IG modes result when both aberrations are present. These modes are therefore associated with ray families for which a linear combination of $T_3^2$ and $T_1$ equals some constant, where $\varepsilon$ determines the ratio between these aberrations:

$$T_3^2 + \frac{\varepsilon}{2} T_1 = \frac{a}{4}. \qquad (8)$$

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

$$t_3^2 + \frac{\varepsilon}{N+1} t_1 = \frac{a}{(N+1)^2}. \qquad (9)$$

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 equation (7) to eliminate $t_3$ from equation (9), giving

$$(t_1 - c)^2 + t_2^2 = R^2, \qquad (10)$$

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

$$c = \frac{\varepsilon}{2(N+1)}, \quad \text{and} \quad R = \sqrt{1 + \left[ \frac{\varepsilon}{2(N+1)} \right]^2 - \frac{a}{(N+1)^2}}. \qquad (11)$$
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Ray structure of an IG beam ($N = 22$, even and $\mu = 18$) as it transitions between the LG ($c = 0$) and HG ($c \to \infty$) 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 $\phi$ 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.</figcaption></figure>

The PP for IG modes then corresponds to the intersection of the ray-Poincaré sphere with this cylinder, as shown in figure 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 figure 3 correspond to the different types of generalized Viviani curves for the PP:

1. For $c + R < 1$ (that is, for $\varepsilon < 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 $\varepsilon = 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 $\varepsilon > 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 $\varepsilon$ increases these two families deform until they share an EFR for $\varepsilon = 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 $\epsilon$ 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 figure 3.

## 3.2. Caustic structure in terms of Keplerian and harmonic oscillator trajectories

As we can see from the bottom row in figure 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 large density of rays in their vicinity (on one side), caustics are associated with the brightest regions of an optical field [39–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 figure 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

$$ (t_1, t_2) = (c + R \cos \varphi, R \sin \varphi), \qquad (12) $$

where $\varphi$ is the angle measured from the positive $t_1$ axis (see figures 4 and 5). In the LG-like regime $\varphi$ traces a complete cycle while in the HG-like regime it is constrained to $[\varphi_0, 2\pi - \varphi_0]$ with

$$ \varphi_0 = \arccos h \quad \text{with} \quad h = \frac{1 - R^2 - c^2}{2Rc}. \qquad (13) $$

We now describe the geometric procedure for finding the caustics from the PP projection. As discussed in [17], the first step is to find what is known as the *medial axes* or *topological skeleton* [51] between the PP projection in equation (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 figure 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 $\varphi$), there are two bitangent circles, one of radius $r_o$ and whose center $\mathbf{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 $\mathbf{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

$$ \|\mathbf{m}_o - (c, 0)\| = r_o + R, \quad \|\mathbf{m}_s - (c, 0)\| = |R - r_s|. \qquad (14) $$

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,

$$ \|\mathbf{m}_o\| = 1 - r_o, \quad \|\mathbf{m}_s\| = 1 - r_s. \qquad (15) $$

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 $\mathbf{m}_o$ the elimination of $r_o$ leads in all cases to the equation

$$ \|\mathbf{m}_o\| + \|\mathbf{m}_o - (c, 0)\| = 1 + R, \qquad (16) $$

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 figure 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 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 $\mathbf{m}_o$, given the fact that the PP projection is an incomplete circle.

The equation for the second medial axis, traced by $\mathbf{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
<figure><img src="image_8.png" /><figcaption>Figure 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 $\phi$ by the point $m_0$ which is the center of the bitangent circle (blue) of radius $r_0$ 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 $\phi$. Note that the PP projection's center and $m_0$ are on opposite sides of the tangent point to the PP projection. (bottom) Same construction (with orange instead of blue) for the medial point $m_s$, 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.</figcaption></figure>

of the second expression in equation (14) equals $r_s - R$. The elimination of $r_s$ from equations (14) and (15) then gives

$$ \|\mathbf{m}_s\| + \|\mathbf{m}_s - (c, 0)\| = 1 - R, \qquad (17) $$

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 figure 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 equation (14) equals $R - r_s$, and the elimination of $r_s$ gives instead

$$ \|\mathbf{m}_s\| - \|\mathbf{m}_s - (c, 0)\| = 1 - R, \qquad (18) $$

which describes a branch of a hyperbola with the same foci, as illustrated in the bottom-right panel of figure 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 figure 4 and in Media 3), it is a straight line for $R = 1$, it opens towards the origin for $R > 1$ (as shown in Media 4), and it becomes a parabola with focus at the origin in the limit $R \to \infty$ (corresponding to the HG case). In fact, in this latter limit, the ellipse traced by $m_0$ 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 figure 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 2.

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_0$ in all cases (although $m_0$ 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),
<figure><img src="image_9.png" /></figure>

repulsive for $R < 1$, and there is no force for $R = 1$. In the limit $c \to \infty$ (and hence $R \to \infty$) corresponding to HG beams, both $m_0$ and $m_s$ trace parabolas corresponding to escape orbits, while in the opposite limit of $c \to 0$ corresponding to LG beams, both $m_0$ 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:

$$Q_x^{(c)} + iQ_y^{(c)} = \pm w_0 \sqrt{(N+1)(t_1 + it_2)} \quad (19)$$

The resulting point $Q^{(c)} = (Q_x^{(c)}, Q_y^{(c)})$ 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, $(\pm 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.

# 4. Semiclassical analysis

## 4.1. 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
<figure><img src="image_10.png" /><figcaption>Figure 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 µ.</figcaption></figure>

some interesting aspects of the topological transition between the LG-like and the HG-like regimes, and elucidates 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 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 $\tau$ and $\eta$. When applied to SG beams, the integral in $\tau$ over a period can be solved analytically, and the constant factors in equation (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]:

$$U(\mathbf{r}) \approx \oint A(\eta) \sqrt{\phi'(\eta) \sin\theta - i\theta'(\eta)} U_N(\theta, \phi; \mathbf{r}) e^{i(N+1)\Upsilon(\eta)} d\eta, \quad (20)$$

where $A(\eta)$ is a weight factor for the contributions from each EFR, primes denote derivatives with respect to $\eta$, $U_N$ is a mathematical analog of a spin coherent state [32, 59–61] defined as

$$U_N(\theta, \phi; \mathbf{r}) = \frac{(\mathbf{v} \cdot \mathbf{v})^{N/2}}{w\sqrt{2^{N-1}\pi N!}} e^{-\frac{\tau^2}{w^2}} H_N\left(\frac{\sqrt{2}\mathbf{v} \cdot \mathbf{r}}{w\sqrt{\mathbf{v} \cdot \mathbf{v}}}\right), \quad (21)$$

with $H_N$ representing an Nth order Hermite polynomial, and

$$\Upsilon(\eta) = \frac{1}{2} \int_{0}^{\eta} \phi'(\bar{\eta}) \cos\theta(\bar{\eta}) d\bar{\eta} \quad (22)$$

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 $\tau$ is guaranteed by the factors in equation (2), the periodicity of the integrand in $\eta$ imposes conditions on the phase contribution due to the geometric phase caused by the PP enclosing a solid angle $\Omega$ over the ray-Poincaré sphere. Because the square-root factor in equation (20) accumulates a phase of $\pm\pi$ on closing the PP loop, the solid angle must then be an odd multiple of $2\pi/(N+1)$, [17, 18]

$$\Omega = \frac{2\pi}{N+1} (2s+1), \quad (23)$$

for integer s. This condition restricts the valid values for R in equation (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\pi/(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\pi/(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 figure 6 there are $N+1$ allowed values of R for the HG-like regime, but only $\lfloor(N+1)/2\rfloor$ for the LG-like regime.
## 4.2. 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 equation (8):

$$ \left( \hat{T}_3^2 + \frac{\varepsilon}{2} \hat{T}_1 \right) \mathrm{IG}_{N,\mu}^{(p)} = \frac{a}{4} \mathrm{IG}_{N,\mu}^{(p)}, \qquad (24) $$

where the operators $\hat{T}_j$ are obtained from the expressions for $T_j$ in equations (5) by substituting $\mathbf{P} \to -ik^{-1}(\partial_x, \partial_y)$. Note that $a$ corresponds now to an eigenvalue. The dotted lines in figure 6 were calculated by using the values of $a$ found through the numerical solution of equation (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 figure 6 for the semiclassical estimate and the result 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 equation (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$.

## 4.3. The semiclassical amplitude and a connection with the simple pendulum

The remaining element for evaluating the semiclassical estimate in equation (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 resembling that of IG beams. However, the form for this function that makes the estimate approach closely the IG beams can be calculated from substituting equation (20) into equation (24), leading to an equation that can be separated into different orders of $N$ (see appendix B for more details). Unsurprisingly, the dominant term gives equation (8), namely the PP. The next term in the series leads to an equation for the amplitude $A(\eta)$, 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(\eta)$:

$$ A^2 = \frac{B}{c} \frac{t_3'}{t_2}, \qquad (25) $$

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:

$$ t_1' = -t_2 t_3, \quad t_2' = t_1 t_3 - ct_3, \quad t_3' = ct_2. \qquad (26) $$

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 $\varphi$ used to parametrize the cylinder in equation (12):

$$ \varphi'' - cR \sin \varphi = 0. \qquad (27) $$

This is the equation for a simple pendulum where the parameter $\eta$ 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
<figure><img src="image_11.png" /><figcaption>Figure 7. Mapping between the phase space curves of the pendulum and the PPs of IG beams on the ray-Poincaré sphere.</figcaption></figure>

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 equation (12) to write $t_3$ in terms of $\varphi$ as $t_3 = \pm\sqrt{1-t_1^2-t_2^2} = \pm\sqrt{1-c^2-R^2-2cR\cos\varphi}$, and substitute this into the third relation in equations (26), which leads to the simple expression

$$ \varphi'(\eta) = t_3. \qquad (28) $$

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 $(\varphi, t_3)$ for the pendulum, in which its motion is described by a curve. However, given the possible periodicity of $\varphi$ 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 figure 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 equation (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 $(\varphi, t_3)$ phase space, becomes $t_3^2 = 2cR(h - \cos\varphi)$, 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| \le 1$.

### 4.4. 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 equation (20) has no closed-form solution, we approximate it as a discrete sum of $M$ terms:

$$ U(\mathbf{r}) \approx \sum_{j=0}^{M-1} \sqrt{\phi'(\eta_j) \sin\theta(\eta_j) - i\theta'(\eta_j)} U_N[\theta(\eta_j), \phi(\eta_j); \mathbf{r}] e^{i(N+1)\Upsilon(\eta_j)}, \qquad (29) $$
<figure><img src="image_12.png" /><figcaption>Figure 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$, $\mu = 4$, $c = 0.1$ and (fourth and fifth columns) $N = 22$, $\mu = 18$, $c = 0.3$, and (second row) even IG beams in the HG-like regime with (second and third columns) $N = 6$, $\mu = 4$, $c = 0.4$ and (fourth and fifth columns) $N = 22$, $\mu = 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 $\eta$. Note that their projection onto the $(t_1, t_2)$ plane is equivalent to snapshots at uniform time intervals of the evolution of a pendulum.</figcaption></figure>

where we use the parametrization defined by equations (26) and (27) that makes $A = 1$. The dependence of the angles $\theta$ and $\phi$ on the parameter $\eta$ is determined by solving these differential equations with appropriate initial conditions determined by the PP. Setting $\eta_0 = 0$ then $\eta_j = j\eta^{(p)}/M$ where the pendulum's period is given by

$$ \eta^{(p)} = \begin{cases} \int_{0}^{2\pi} \frac{d\eta}{\sqrt{1-c^2-R^2-2cR\cos\eta}} = 4\frac{K(h)}{\sqrt{1-(R+c)^2}}, & \text{if } c+R < 1, \\ \int_{\varphi_0}^{2\pi-\varphi_0} \frac{2d\eta}{\sqrt{1-c^2-R^2-2cR\cos\eta}} = 8\frac{F(\varphi_0/2|h)-K(h)}{\sqrt{1-(R+c)^2}}, & \text{if } c+R > 1, \end{cases} \quad (30) $$

with $h = -4Rc/[1 - (R+c)^2]$, 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 $\mu$ as $IG_{N,\mu}^{(\pm)} = IG_{N,\mu}^{(e)} \pm iIG_{N,\mu}^{(o)}$ [21]. As shown in figure 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 equation (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 figure 8 and for others with a lower value for $\mu$ for different values of $c$ (and hence of $\varepsilon$). 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$.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Rms error $\epsilon_{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.</figcaption></figure>

**5. 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 [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. Cavities play a central role not only in the generation of structured laser beams, but also in advanced metrology applications such as gravitational wave detection [63, 64]. Understanding the modes of these cavities, especially in the presence of the most basic aberrations, is of central importance. Since some amount of these aberrations are typically present, the modes produced or supported by realistic cavities are more likely approximations to IG beams that might be near the HG or the LG limits. The ray treatment provided here sheds light on the structure of the beams from the ray-optical perspective, as we now discuss. Imagine a subset of rays corresponding to a single EFR like those shown in figure 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 $\tau$. 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 $\tau$ and a slow pendulum-like precession in $\eta$. From this point of view, we see that the connection between IG modes and the pendulum does not arise in the semiclassical 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. This type of ray-based calculation is not only intuitive but also typically much faster [65–67]. However, it has been so far limited to simple focused beams, while calculations of trapping forces and torques for structured beams [68–70], including IG beams [26], have relied on the wave model.

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.

## Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/rodguti90/InceGauss. [71].

## Acknowledgment

The authors dedicate this work to the memory of Kurt Bernardo Wolf and Willard Miller, Jr., both of whom passed recently. MRD and MAA first presented aspects of this work almost exactly a decade ago at a workshop to celebrate Wolf's 70th birthday and where Miller was also present. Discussions with both of them led us to a deeper appreciation of these results. The authors also acknowledge J M Fellows and D H J O'Dell for useful discussions. R G-C and M A A acknowledge the Excellence Initiative of Aix-Marseille University—A*MIDEX, a French 'Investissements d'Avenir' programme for funding this research. R G C also acknowledges funding from the Labex WIFI (ANR-10-LABX-24, ANR-10-IDEX-0001-02 PSL*). MRD acknowledges support from the EPSRC Centre for Doctoral Training in Topological Design (EP/S02297X/1).

## Conflict of interest

The authors declare no conflicts of interest.

## Appendix A. Analytic formula for equation (22)

As mentioned in the main text, it is possible to perform the integral defining the phase $\Upsilon$ analytically for the ray families describing IG fields. The result is given by the following expression in terms of elliptic functions:

$$ \Upsilon(\eta) = \mathrm{Re} \left\{ \tan \left[ \frac{\phi(\eta)}{2} \right] \sqrt{1 - c^2 - R^2 - 2cR \cos \phi(\eta)} + \frac{i \operatorname{sgn}(C_+)}{\sqrt{C_+}} \left[ C_+ E \left( \gamma \left| \frac{C_-}{C_+} \right. \right) + \frac{2R(1+c^2-R^2)}{c-R} F \left( \gamma \left| \frac{C_-}{C_+} \right. \right) - \frac{4cR}{c^2-R^2} \Pi \left( \frac{1-C_-}{1-C_+}; \gamma \left| \frac{C_-}{C_+} \right. \right) \right] \right\}, \quad (A1) $$

where $C_{\pm} = 1 - (c \pm R)^2$, $\gamma = i \arcsinh\{\tan[\phi(\eta)/2]\}$, E is the elliptic integral of the second kind, F elliptic integral of the first kind, and $\Pi$ 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 $\phi$ within the interval $[-\pi, \pi]$ such that it is zero for $\phi = 0$ in the LG-like regime and for $\phi = \pm \phi_0$ in the HG-like regime. There is therefore a discontinuity at $\phi = \pm \pi$. 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\pi$. For the HG-like regime, on the other hand, the discontinuity must be removed by hand to make the function continuous at $\phi = \pm \pi$. Unfortunately the limit $\phi \to \pm \pi$ 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(\mathbf{r}; \mathbf{Q}, \mathbf{P}) = \sqrt{\frac{k\Gamma}{2\pi}} \exp \left[ -\frac{k\Gamma}{2} \|\mathbf{r} - \mathbf{Q}\|^2 + ik(\mathbf{r} - \mathbf{Q}) \cdot \mathbf{P} \right] e^{ikL(\eta, \tau)}. \quad (B1)$$

The center **Q** and direction **P** of the Gaussians are the same as those of the corresponding rays. $\Gamma$ 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, $\Gamma = 2/kw_0^2$. 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(\mathbf{r}) = \iint A(\eta) \sqrt{J(\tau, \eta)} g(\mathbf{r}; \mathbf{Q}, \mathbf{P}) \, d\tau \, d\eta \quad (B2)$$

where $L(\eta, \tau)$ is the optical path length that plays the role of the action in mechanics. The weight $A$ only depends on the parameter $\eta$ 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,

$$\iint A(\eta) \sqrt{J(\tau, \eta)} \left[ \hat{T}_3^2 g(\mathbf{r}; \mathbf{Q}, \mathbf{P}) + \frac{\epsilon}{2} \hat{T}_1 g(\mathbf{r}; \mathbf{Q}, \mathbf{P}) - \frac{a}{4} g(\mathbf{r}; \mathbf{Q}, \mathbf{P}) \right] d\tau d\eta = 0. \quad (B3)$$

Since $g$ is the only part that has a spatial dependence, we start by computing the following relations:

$$\hat{T}_1 g(\mathbf{r}, \eta, \tau) = \left[ \frac{N+1}{2} (|z_x|^2 - |z_y|^2) + \frac{\sqrt{N+1}}{w} (z_x \Delta_x - z_y \Delta_y) \right] g(\mathbf{r}, \eta, \tau) \quad (B4a)$$

$$\hat{T}_3^2 g(\mathbf{r}, \eta, \tau) = \left[ (N+1)^2 (p_x q_y - p_y q_x)^2 + \frac{N+1}{2} \mathbf{q} \cdot \mathbf{z} - \frac{2i(N+1)^{3/2}}{w} (p_x q_y - p_y q_x) (z_x \Delta_y - z_y \Delta_x) \right. \\ \left. + \frac{N+1}{w^2} (2z_x z_y \Delta_x \Delta_y - z_x^2 \Delta_y^2 - z_y^2 \Delta_x^2) + \frac{\sqrt{N+1}}{2w} \mathbf{z} \cdot \mathbf{\Delta} \right] g(\mathbf{r}, \eta, \tau) \quad (B4b)$$

where we introduced $\mathbf{\Delta} = (\Delta_x, \Delta_y) = \mathbf{r} - \mathbf{Q}$,

$$\mathbf{Z} = \Gamma \mathbf{Q} + i \mathbf{P} = \frac{2}{wk} \sqrt{N+1} \mathbf{v}(\eta) e^{-i\tau} = \frac{2}{wk} \sqrt{N+1} (\mathbf{q} + i\mathbf{p}) = \frac{2}{wk} \sqrt{N+1} \mathbf{z} \quad (B5)$$

and wrote explicitly the dependence of $g$ on $\eta$ and $\tau$.

The next step is to move all the spatial dependence to $g$, in order to be able to identify part of the integrand as being equal to zero. This is done by using the following identities,

$$\Delta_x g = \frac{w}{2(N+1)^{1/2}} \left( \frac{iz_y}{j} \frac{\partial g}{\partial \eta} + \frac{z'_y}{j} \frac{\partial g}{\partial \tau} \right), \quad \Delta_y g = - \frac{w}{2(N+1)^{1/2}} \left( \frac{iz_x}{j} \frac{\partial g}{\partial \eta} + \frac{z'_x}{j} \frac{\partial g}{\partial \tau} \right) \quad (B6)$$

where primes indicate derivatives with respect to $\eta$ (all derivatives with respect to $\tau$ can be taken explicitly because the dependence of **Q** and **P** on $\tau$ is known) and

$$J(\tau, \eta) = \frac{\partial \mathbf{Z}}{\partial (\tau, \eta)} = \frac{\partial Z_x}{\partial \tau} \frac{\partial Z_y}{\partial \eta} - \frac{\partial Z_y}{\partial \tau} \frac{\partial Z_x}{\partial \eta} = \frac{4i(N+1)}{w^2 k^2} \left( z_y \frac{\partial z_x}{\partial \eta} - z_x \frac{\partial z_y}{\partial \eta} \right) = \frac{4(N+1)}{w^2 k^2} j(\tau, \eta), \quad (B7)$$
and then integrating by parts using the fact that the integrations are periodic. This leads to the following identities which apply to all the terms appearing in equation (B4)

$$ 
\begin{aligned}
\iint AF\sqrt{j}g d\eta d\tau &= \frac{2(N+1)^{1/2}}{wk} \iint AF\sqrt{j}g d\eta d\tau \\
\iint AF\sqrt{j}\Delta_x g d\eta d\tau &= -\frac{1}{k} \iint \left[ iF'z_y \frac{A}{\sqrt{j}} + iFz_y \left( \frac{A}{\sqrt{j}} \right)' \right] g d\eta d\tau \\
\iint AF\sqrt{j}\Delta_y g d\eta d\tau &= \frac{1}{k} \iint \left[ iF'z_x \frac{A}{\sqrt{j}} + iFz_x \left( \frac{A}{\sqrt{j}} \right)' \right] g d\eta d\tau
\end{aligned}
$$

where $F(\eta, \tau) = f(\eta) \exp(-i\tau)$, with $f$ being an arbitrary function of $\eta$, and for the rest of the terms we can be more specific,

$$ 
\begin{aligned}
\iint A\sqrt{j}z \cdot \Delta g d\eta d\tau &= -\frac{1}{k} \iint j \frac{A}{\sqrt{j}} g d\eta d\tau \\
\iint Az_y^2 \sqrt{j} \Delta_x^2 g d\eta d\tau &= \frac{w(N+1)^{1/2}}{k} \iint \left\{ z_y^2 \frac{A}{\sqrt{j}} \left[ iq_x'z_y + p_x z_y' \right] \right\} g d\eta d\tau \\
&\quad - \frac{w}{2k(N+1)^{1/2}} \iint \left\{ \frac{z_y}{3} \left[ \left( \frac{z_y^3}{j} \right)' \right]' \left( \frac{A}{\sqrt{j}} \right) + \left( \frac{z_y^4}{j} \right)' \left( \frac{A}{\sqrt{j}} \right)' \right. \\
&\quad \left. + \left( \frac{z_y^4}{j} \right) \left( \frac{A}{\sqrt{j}} \right)'' \right\} g d\eta d\tau \\
\iint Az_x^2 \sqrt{j} \Delta_y^2 g d\eta d\tau &= -\frac{w(N+1)^{1/2}}{k} \iint \left\{ z_x^2 \frac{A}{\sqrt{j}} \left[ iq_x'z_x + p_y z_x' \right] \right\} g d\eta d\tau \\
&\quad - \frac{w}{2k(N+1)^{1/2}} \iint \left\{ \frac{z_x}{3} \left[ \left( \frac{z_x^3}{j} \right)' \right]' \left( \frac{A}{\sqrt{j}} \right) + \left( \frac{z_x^4}{j} \right)' \left( \frac{A}{\sqrt{j}} \right)' \right. \\
&\quad \left. + \left( \frac{z_x^4}{j} \right) \left( \frac{A}{\sqrt{j}} \right)'' \right\} g d\eta d\tau
\end{aligned}
$$

$$ 
\begin{aligned}
\iint Az_x z_y \sqrt{j} \Delta_x \Delta_y g d\eta d\tau &= \frac{w(N+1)^{1/2}}{k} \iint \frac{i}{2} z_x z_y (|z_y|^2)' \frac{A}{\sqrt{j}} g d\eta d\tau \\
&\quad + \frac{w}{2k(N+1)^{1/2}} \iint \left\{ \frac{z_x z_y}{2j^2} \left\{ j \left( z_x'' z_y + z_x z_y'' \right) + 4jz_x' z_y' - j' \left( z_x z_y \right)' \right\} \frac{A}{\sqrt{j}} \right. \\
&\quad \left. + \left( \frac{z_x^2 z_y^2}{j} \right)' \left( \frac{A}{\sqrt{j}} \right)' + \left( \frac{z_x^2 z_y^2}{j} \right)'' \right\} g d\eta d\tau.
\end{aligned}
$$

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 + \frac{\varepsilon N+1}{2} (|z_x|^2 - |z_y|^2) = \frac{a}{4}. \qquad (B8) $$

Noting that

$$ T_3^2 = (N+1)^2 (p_x q_y - p_y q_x)^2 \qquad T_1 = \frac{N+1}{2} (|z_x|^2 - |z_y|^2) \qquad (B9) $$
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 = \left\{ \left[ (p_x q_y - p_y q_x) (z_x^2 + z_y^2) \right]' \frac{A}{\sqrt{j}} + 2 \left[ (p_x q_y - p_y q_x) (z_x^2 + z_y^2) \right] \left( \frac{A}{\sqrt{j}} \right)' \right\} - \frac{\varepsilon i}{2(N+1)} \left[ (z_x z_y)' \frac{A}{\sqrt{j}} + 2z_x z_y \left( \frac{A}{\sqrt{j}} \right)' \right].$$

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 = \frac{B_j}{(p_x q_y - p_y q_x) (z_x^2 + z_y^2) - \frac{\varepsilon 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 $\eta$, as can be appreciated when writing them in terms of $\theta$ and $\phi$,

$$(N+1)^2 \cos^2 \theta + \varepsilon (N+1) \cos \phi \sin \theta = a, \qquad (B10a)$$

$$A^2 = \frac{-\sin\theta\phi' + i\theta'}{\cos\theta\sin\theta + \frac{\varepsilon i}{2(N+1)}(i\cos\phi\cos\theta + \sin\phi)}. \qquad (B10b)$$

Taking the derivative of the first equation

$$\cos\theta\sin\theta = \frac{\varepsilon}{2(N+1)} \left( -\frac{\phi'}{\theta'} \sin\phi \sin\theta + \cos\phi \cos\theta \right) \qquad (B11)$$

and substituting it into the second one we get

$$A^2 = \frac{2(N+1)}{\varepsilon} \frac{\theta'}{\sin\phi}. \qquad (B12)$$

Finally, by writing this equation in terms of $\phi$, we have that

$$A^2 = \pm \frac{\phi'}{\sqrt{1 - c^2 - R^2 - 2cR \cos \phi}}, \qquad (B13)$$

where we used the relation $c = \varepsilon/2(N+1)$. The denominator is proportional to the angular velocity of a simple pendulum, therefore if we set $\phi = \eta$ 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 $\phi$ follows the dynamics of the pendulum where $\eta$ plays the role of time. The connection to the simple pendulum becomes clearer after taking the second derivative of the previous equation, leading to

$$\phi'' - cR \sin \phi = 0. \qquad (B14)$$

This is the defining equation of a simple pendulum

## ORCID iDs

R Gutiérrez-Cuevas https://orcid.org/0000-0002-3451-6684

M R Dennis https://orcid.org/0000-0003-1147-1804

M A Alonso https://orcid.org/0000-0001-7037-5383
## References

[1] Boyer C P, Kalnins E G and Miller W 1975 Lie theory and separation of variables. 6. The equation $i\dot{u}_t + \delta_2 u = 0$. *J. Math. Phys.* **16** 499

[2] Siegman A E 1986 *Lasers* (University Science Books)

[3] Siviloglou G A, Broky J, Dogariu A and Christodoulides D N 2007 Observation of accelerating airy beams *Phys. Rev. Lett.* **99** 213901

[4] Bandres M A 2008 Accelerating parabolic beams *Opt. Lett.* **33** 1678

[5] Levy U, Derevyanko S and Silberberg Y 2016 Light modes of free space *Prog. Opt.* **61** 237–81

[6] Dennis M R and Ring J D 2013 Propagation-invariant beams with quantum pendulum spectra: from bessel beams to gaussian beam-beams *Opt. Lett.* **38** 3325

[7] Gutiérrez-Cuevas R and Alonso M A 2017 Polynomials of Gaussians and vortex-Gaussian beams as complete, transversely confined bases *Opt. Lett.* **42** 2205

[8] Gutiérrez-Cuevas R and Alonso M A 2017 Complete confined bases for beam propagation in Cartesian coordinates *J. Opt. Soc. Am. A* **34** 1697

[9] Boyer C P, Kalnins E G and Miller W 1975 Lie theory and separation of variables. 7. The harmonic oscillator in elliptic coordinates and ince polynomials *J. Math. Phys.* **16** 512

[10] Andrews D L 2008 *Structured Light and Its Applications: An Introduction to Phase-Structured Beams and Nanoscale Optical Forces* (Academic PR Inc)

[11] Berkhout G C G, Lavery M P J, Courtial J, Beijersbergen M W and Padgett M J 2010 Efficient sorting of orbital angular momentum states of light *Phys. Rev. Lett.* **105** 153601

[12] Zhou Y, Mirhosseini M, Fu D, Zhao J, Rafsanjani S M H, Willner A E and Boyd R W 2017 Sorting photons by radial quantum number *Phys. Rev. Lett.* **119** 263602

[13] Gu X, Krenn M, Erhard M and Zeilinger A 2018 Gouy phase radial mode sorter for light: Concepts and experiments *Phys. Rev. Lett.* **120** 103601

[14] Tsang M, Nair R and Lu X-M 2016 Quantum theory of superresolution for two incoherent optical point sources *Phys. Rev. X* **6** 031033

[15] Yao A M and Padgett M J 2011 Orbital angular momentum: origins, behavior and applications *Adv. Opt. Photon.* **3** 161

[16] Berry M V and McDonald K T 2008 Exact and geometrical optics energy trajectories in twisted beams *J. Opt. A: Pure Appl. Opt.* **10** 035005

[17] Alonso M A and Dennis M R 2017 Ray-optical Poincaré sphere for structured Gaussian beams *Optica* **4** 476

[18] Dennis M R and Alonso M A 2019 Gaussian mode families from systems of rays *J. Phys. Photon.* **1** 025003

[19] Arscott F M 1964 *Periodic Differential Equations* (Pergamon)

[20] Bandres M A and Gutiérrez-Vega J C 2004 Ince–Gaussian beams *Opt. Lett.* **29** 144

[21] Bandres M A and Gutiérrez-Vega J C 2004 Ince–Gaussian modes of the paraxial wave equation and stable resonators *J. Opt. Soc. Am. A* **21** 873

[22] Yao-Li Y L, Hu X-B, Perez-Garcia B, Bo-Zhao B Z, Gao W, Zhu Z-H and Rosales-Guzmán C 2020 Classically entangled Ince–Gaussian modes *Appl. Phys. Lett.* **116** 221105

[23] Sakpal S, Milione G, Li M-J, Nouri M, Shahoei H, LaFave T, Ashrafi S and MacFarlane D 2018 Stability of ince-gaussian beams in elliptical core few-mode fibers *Opt. Lett.* **43** 2656

[24] Gutiérrez-Cuevas R, O'Dell D H J, Dennis M R and Alonso M A 2023 Exactly solvable model behind Bose-Hubbard dimers, Ince-Gauss beams and aberrated optical cavities *Phys. Rev. A* **107** L031502

[25] Dennis M R 2006 Rows of optical vortices from elliptically perturbing a high-order beam *Opt. Lett.* **31** 1325

[26] Woerdemann M, Alpmann C and Denz C 2011 Optical assembly of microparticles into highly ordered structures using Ince-Gaussian beams *Appl. Phys. Lett.* **98** 111101

[27] Krenn M, Fickler R, Huber M, Lapkiewicz R, Plick W, Ramelow S and Zeilinger A 2013 Entangled singularity patterns of photons in ince-gauss modes *Phys. Rev. A* **87** 012326

[28] Plick W N, Krenn M, Fickler R, Ramelow S and Zeilinger A 2013 Quantum orbital angular momentum of elliptically symmetric light *Phys. Rev. A* **87** 033806

[29] Dennis M R and Alonso M A 2017 Swings and roundabouts: optical Poincaré spheres for polarization and Gaussian beams *Phil. Trans. R. Soc. A* **375** 20150441

[30] Malhotra T, Gutiérrez-Cuevas R, Hassett J, Dennis M, Vamivakas A and Alonso M 2018 Measuring geometric phase without interferometry *Phys. Rev. Lett.* **120** 233602

[31] Gutiérrez-Cuevas R, Dennis M R and Alonso M A 2019 Generalized Gaussian beams in terms of Jones vectors *J. Opt.* **21** 084001

[32] Gutiérrez-Cuevas R, Wadood S, Vamivakas A and Alonso M 2020 Modal Majorana sphere and hidden symmetries of structured-Gaussian beams *Phys. Rev. Lett.* **125** 123903

[33] Gutiérrez-Cuevas R and Alonso M A 2020 Platonic Gaussian beams: wave and ray treatment *Opt. Lett.* **45** 6759

[34] Schwarz U T, Bandres M A and Gutiérrez-Vega J C 2004 Observation of Ince-Gaussian modes in stable resonators *Opt. Lett.* **29** 1870

[35] Sakurai J J and Napolitano J J 2010 *Modern Quantum Mechanics* 2nd edn (Pearson)

[36] Padgett M J and Courtial J 1999 Poincaré-sphere equivalent for light beams containing orbital angular momentum *Opt. Lett.* **24** 430

[37] Calvo G F 2005 Wigner representation and geometric transformations of optical orbital angular momentum spatial modes *Opt. Lett.* **30** 1207

[38] Gati R and Oberthaler M K 2007 A bosonic Josephson junction *J. Phys. B: At. Mol. Opt. Phys.* **40** R61

[39] Berry M and Upstill C 1980 IV catastrophe optics: Morphologies of caustics and their diffraction patterns *Prog. Opt.* **18** 257–346

[40] Nye J 1999 *Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations* (Institute of Physics Publishing)

[41] Kravtsov Y A and Orlov Y I 1999 *Caustics, Catastrophes and Wave Fields* (Springer)

[42] Alonso M and Forbes G 2002 Stable aggregates of flexible elements give a stronger link between rays and waves *Opt. Express* **10** 728

[43] Alonso M A 2009 Rays and waves *Phase Space Optics: Fundamentals and Applications* (McGraw Hill Professional) pp 237–77

[44] Shen Y, Yang X, Naidoo D, Fu X and Forbes A 2020 Structured ray-wave vector vortex beams in multiple degrees of freedom from a laser *Optica* **7** 820

[45] Abramochkin E G and Volostnikov V G 2010 Generalized hermite-laguerre-gauss beams *Phys. Wave Phenom.* **18** 14

[46] lin Chao S 1974 High order transverse modes in an astigmatic cavity *PhD Thesis* University of Rochester

[47] Abbena E, Salamon S and Gray A 2017 *Modern Differential Geometry of Curves and Surfaces With Mathematica* (CRC Press)
[48] Graefe E-M, Korsch H J and Strzys M P 2014 Bose-Hubbard dimers, Viviani's windows and pendulum dynamics *J. Phys. A: Math. Theor.* **47** 085304

[49] Morales J D H, Rodríguez-Lara B M and Malomed B A 2017 Polarization dynamics in twisted fiber amplifiers: a non-hermitian nonlinear dimer model *Opt. Lett.* **42** 4402

[50] Azzam R M A 2000 Poincaré sphere representation of the fixed-polarizer rotating-retarder optical system *J. Opt. Soc. Am. A* **17** 2105

[51] Blum H 1967 A transformation for extracting new descriptors of shape *Models for the Perception of Speech and Visual Form* ed W Wathen-Dunn (MIT Press) pp 362–80

[52] Ogilvy C S 1990 *Excursions in Geometry* (Courier Corporation)

[53] Mittag L and Stephen M J 1992 Conformal transformations and the application of complex variables in mechanics and quantum mechanics *Am. J. Phys* **60** 207

[54] Alemán-Casta neda L A and Alonso M A 2019 Study of reflectors for illumination via conformal maps *Opt. Lett.* **44** 3809

[55] Zor D and Kay K G 1996 Globally uniform semiclassical expressions for time-independent wave functions *Phys. Rev. Lett.* **76** 1990

[56] Alonso M A and Forbes G W 2001 Using rays better II ray families to match prescribed wave fields *J. Opt. Soc. Am. A* **18** 1146

[57] Kay K G 1999 Exact wave functions for the coulomb problem from classical orbits *Phys. Rev. Lett.* **83** 5190

[58] Kay K G 2001 Exact wave functions from classical orbits: the isotropic harmonic oscillator and semiclassical applications *Phys. Rev. A* **63** 042110

[59] Radcliffe J M 1971 Some properties of coherent spin states *J. Phys. A: Gen. Phys.* **4** 313

[60] Arecchi F T, Courtens E, Gilmore R and Thomas H 1972 Atomic coherent states in quantum optics *Phys. Rev. A* **6** 2211

[61] Perelomov A M 1972 Coherent states for arbitrary lie group *Commun. Math. Phys.* **26** 222

[62] Forbes G W and Alonso M A 2023 Asymptotic expansions for field moments of bound states *J. Phys. A: Math. Theor.* **56** 244001

[63] Abbott B P et al (LIGO Scientific Collaboration and Virgo Collaboration) 2016 Observation of gravitational waves from a binary black hole merger *Phys. Rev. Lett.* **116** 061102

[64] Tao L, Green A and Fulda P 2020 Higher-order Hermite-Gauss modes as a robust flat beam in interferometric gravitational wave detectors *Phys. Rev. D* **102** 122002

[65] Ashkin A 1992 Forces of a single-beam gradient laser trap on a dielectric sphere in the ray optics regime *Biophys. J.* **61** 569

[66] Ashkin A 2000 History of optical trapping and manipulation of small-neutral particle, atoms and molecules *IEEE J. Sel. Top. Quantum Electron.* **6** 841

[67] Volke-Sepúlveda K, Chávez-Cerda S, Garcés-Chávez V and Dholakia K 2004 Three-dimensional optical forces and transfer of orbital angular momentum from multiringed light beams to spherical microparticles *J. Opt. Soc. Am. B* **21** 1749

[68] Barton J P, Alexander D R and Schaub S A 1989 Theoretical determination of net radiation force and torque for a spherical particle illuminated by a focused laser beam *J. Appl. Phys.* **66** 4594

[69] Gutiérrez-Cuevas R, Moore N J and Alonso M A 2018 Lorenz-Mie scattering of focused light via complex focus fields: an analytic treatment *Phys. Rev. A* **97** 053848

[70] Gutiérrez-Cuevas R and Alonso M A 2017 Scalar and electromagnetic nonparaxial bases composed as superpositions of simple vortex fields with complex foci *Opt. Express* **25** 14856

[71] Gutiérrez-Cuevas R InceGauss (available at: https://github.com/rodguti90/InceGauss)