# Geometric phases and coherent states for commensurate bichromatic polarization and astigmatic cavities

Miguel A Alonso

► To cite this version:

Miguel A Alonso. Geometric phases and coherent states for commensurate bichromatic polarization and astigmatic cavities. Nanophotonics, 2025, 10.1515/nanoph-2025-0347. hal-05297984

HAL Id: hal-05297984

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

Submitted on 5 Oct 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research Article

Miguel A. Alonso*

# Geometric phases and coherent states for commensurate bichromatic polarization and astigmatic cavities

https://doi.org/10.1515/nanoph-2025-0347

Received July 27, 2025; accepted September 5, 2025;

published online September 22, 2025

**Abstract:** Both the polarization state of coherent bichromatic fields produced by harmonic generation and a class of anisotropic paraxial optical cavities are examples of commensurate two-dimensional harmonic oscillators. The geometric phase for these systems is studied here, both in the classical/ray and quantum/wave regimes. The quantum geometric phase is described in terms of the coherent states of the system, for which recursive expressions are derived that yield the exact result and are numerically stable even for high modal orders.

**Keywords:** geometric phase; bichromatic polarization; coherent states; anisotropic harmonic oscillator

## 1 Introduction

A commensurate anisotropic harmonic oscillator is a system consisting of two or more degrees of freedom, where each is described by harmonic oscillations of different but commensurate frequencies. There are at least two physical situations in paraxial optics that correspond to commensurate two-dimensional anisotropic harmonic oscillators. The first is that of laser cavities with specific amounts of astigmatism [1], or equivalently to gradient index waveguides with appropriate anisotropic refractive index profiles. Here, the “classical” description corresponds to the behavior of rays under propagation, while the analogues of the quantum eigenstates are the transverse profiles of the wave modes. It has been shown that these cavities accept modes that take the shape of Lissajous curves [1]–[5], and that specific cases accept modes whose shapes happen to be separable

in parabolic coordinates [6], [7]. These systems have been implemented experimentally by introducing appropriate anisotropic focusing elements in a laser cavity [3], [4], [7]. The second optical realization is not in the spatial field profile but in its polarization state, and corresponds to coherent bichromatic fields with commensurate frequencies [8]–[11], obtainable for example through nonlinear harmonic generation [12]–[15]. Such fields have been shown theoretically and experimentally to exhibit many polarization distribution patterns and singularities [8]–[11], [16]. The classical trajectories are then simply the paths traced by the electric field vector as a function of time, while the quantum eigenstates correspond to exotic states of quantum polarization. Of course, these two types of realization are completely different physically, as are the experimental techniques needed to implement and measure them, and the type of information a measurement can obtain. Nevertheless, their shared underlying model reveals analogous topological features, potentially leading to cross-fertilization of theoretical and experimental ideas.

Commensurate two-dimensional anisotropic harmonic oscillators are a simple yet rich mechanical model that accepts closed-form solutions both in the classical and quantum regimes due to their factorizable nature. They are an example of what has been termed “accidental degeneracy” [17], [18], coming from different combinations between the eigenvalues of the factorizable parts that provide the same sum. They are described by the Hamiltonian

$$ \hat{H} = -\frac{1}{2k^2}(\partial_x^2 + \partial_y^2) + V(x, y), \quad (1) $$

$$ V(x, y) = V_2 \frac{m_x^2 x^2 + m_y^2 y^2}{2}, \quad (2) $$

where $x = (x, y)$ are two spatial variables, the constants $m_x$ and $m_y$ are nonnegative integers that do not share any prime factors, and $k = 1/\hbar$ (or $2\pi$ divided by the wavelength in paraxial optics) is a constant that could be set to unity by using appropriate units but that we retain for asymptotic bookkeeping. Henceforth we choose units for which the mass and the potential coefficient $V_2$ equal unity.

<footnote>*Corresponding author: Miguel A. Alonso, Aix Marseille Univ, CNRS, Centrale Med, Institut Fresnel, UMR 7249, 13397, Marseille, Cedex 20, France; and The Institute of Optics, University of Rochester, Rochester, NY, 14627, USA, E-mail: miguel.alonso@fresnel.fr.</footnote>

<footnote>https://orcid.org/0000-0001-7037-5383</footnote>
In this work we study how the concept of geometric (or Pancharatnam-Berry) phase [19] extends to these systems, both in the classical (or ray) and quantum (or wave) domains, and whether this phase corresponds to an area or solid angle over the manifold that describe the different Lissajous oscillation shapes. We also provide a method for calculating exactly and robustly the coherent states of the system associated with each classical oscillation shape. These coherent states serve precisely as the basis for the definition of the quantum geometric phase, which for anisotropic oscillators is shown to be approximately proportional to the classical one only in the limit of highly excited states.

## 2 Classical description

We consider first the classical solutions, which correspond to closed orbits of Lissajous shape for which the position vector $Q = (Q_x, Q_y)$ is composed of two sinusoids of different frequencies. These orbits can be parametrized as

$$Q = \left[ \begin{array}{l} \frac{Q_0}{m_x} \cos \frac{\theta}{2} \cos \left( m_x t + \frac{\phi}{2m_y} \right), \\ \frac{Q_0}{m_y} \sin \frac{\theta}{2} \cos \left( m_y t - \frac{\phi}{2m_x} \right) \end{array} \right], \quad (3)$$

where $Q_0$ is a constant, $\theta$ determines the relative amplitude of the oscillations in both directions, and $\phi$ determines the relative phase delay. Note that the orbits are fully covered in a period in $t$ of $2\pi$. In this respect, varying $\phi$ over an interval of size $2\pi$ is sufficient to cover the full range of orbit shapes for the given constant $\theta$. The momentum $\mathbf{P} = (P_x, P_y) = \partial_t Q$ is given by

$$\mathbf{P} = - \left[ \begin{array}{l} Q_0 \cos \frac{\theta}{2} \sin \left( m_x t + \frac{\phi}{2m_y} \right), \\ Q_0 \sin \frac{\theta}{2} \sin \left( m_y t - \frac{\phi}{2m_x} \right) \end{array} \right]. \quad (4)$$

There are several independent combinations of $Q$ and $\mathbf{P}$ that are constants of the motion, and correspond to generalizations of the Stokes parameters for monochromatic polarization. The first corresponds to the oscillator's energy:

$$T_0 = m_x^2 \frac{Q_x^2}{2} + m_y^2 \frac{Q_y^2}{2} + \frac{P_x^2}{2} + \frac{P_y^2}{2} = \frac{Q_0^2}{2}. \quad (5)$$

A second invariant is a measure of the mismatch in oscillation amplitudes in both directions/frequencies:

$$T_1 = m_x^2 \frac{Q_x^2}{2} - m_y^2 \frac{Q_y^2}{2} + \frac{P_x^2}{2} - \frac{P_y^2}{2} = \frac{Q_0^2}{2} \cos \theta. \quad (6)$$

The remaining two invariants are not quadratic forms, given the different frequencies of oscillation in the two

directions. They can be constructed as the real ($T_2$) and imaginary ($T_3$) parts of the combination

$$T_2 + i T_3 \propto (m_x Q_x - i P_x)^{m_y} (m_y Q_y + i P_y)^{m_x} \\ = Q_0^{m_x+m_y} \cos^{m_y} \frac{\theta}{2} \sin^{m_x} \frac{\theta}{2} \exp(i\phi). \quad (7)$$

We can use appropriate powers of $T_0$ to normalize the remaining parameters $T_n$ according to

$$\begin{aligned} \vec{\tau} &= (\tau_1, \tau_2, \tau_3) \\ &= \left( \cos \theta, C \cos^{m_y} \frac{\theta}{2} \sin^{m_x} \frac{\theta}{2} \cos \phi, C \cos^{m_y} \frac{\theta}{2} \sin^{m_x} \frac{\theta}{2} \sin \phi \right), \end{aligned} \quad (8)$$

where the constant $C$ can be chosen, say, so that the maximal value for $(\tau_2^2 + \tau_3^2)^{1/2}$ is unity:

$$C = \sqrt{\frac{(m_x + m_y)^{m_x+m_y}}{m_x^{m_x} m_y^{m_y}}}. \quad (9)$$

The vector $\vec{\tau}$ then defines a manifold known as a Kummer shape [20], which in the case $m_x = m_y = 1$ reduces to the Poincaré sphere, with the invariants $T_n$ then becoming the Stokes parameters. Each point over this manifold corresponds to an orbit shape, namely a Lissajous curve with a given sense. Figure 1 shows some examples of Kummer shapes, as well as some of the orbits for one of them. By construction, Kummer shapes have rotational symmetry around the $\tau_1$ axis, and one can see that they are widest at $\tau_1 = \epsilon$, where

$$\epsilon = \frac{m_y - m_x}{m_x + m_y}. \quad (10)$$

We refer to this circular cross-section as their equator. Note that (except for the degenerate case $m_x = m_y$ corresponding to the Poincaré sphere) the only geodesics for the Kummer shapes that are contained in planes are precisely this equator as well as the curves of constant $\phi$ (denoted as meridians).

It is worth stressing that only for the isotropic case ($m_x = m_y = 1$) does the invariant $T_3$ represent an angular momentum (which is orbital when we are describing cavity modes and spin when we are considering polarization). In fact, if we define the mean angular momentum of the orbit as $L = (2\pi)^{-1} \int_0^{2\pi} (Q_x P_y - Q_y P_x) dt$, one can easily show that this quantity vanishes exactly except in the isotropic case, for which it gives $L = (Q_0^2/2) \sin \theta \sin \phi$. For anisotropic oscillators, while $T_3$ does not correspond to an angular momentum, it does share some of its qualitative features: the orbit's shape depends on the magnitude of $T_3$ and not on its sign, but this sign determines the sense of circulation. The subset of orbits for which $T_3 = 0$ (such as those shown in Figure 1c as black curves) are those for which the orbit oscillates back and forth along the same path, so they are
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Kummer shapes (orange) for anisotropic oscillators with (a) $m_x = 1, m_y = 2$, (b) $m_x = 1, m_y = 3$, (c) $m_x = 2, m_y = 3$, (d) $m_x = 21, m_y = 23$. The equator of these shapes is shown as a black circle whose radius is unity, so the shapes can be inscribed in a cylinder of unit radius and height two. In (c) some of the orbits are shown, corresponding to the intersections of the axes with the Kummer shape. Note also that the pole $\vec{\tau} = (1, 0, 0)$ is smooth only for $m_x = 1$, and correspondingly the pole $\vec{\tau} = (-1, 0, 0)$ is smooth only for $m_y = 1$. In fact, the case $m_x = m_y = 1$ (not shown) corresponds to the standard Poincaré sphere.</figcaption></figure>

the analogues of linear polarizations for the isotropic case (although they are only strictly linear when $T_2$ also vanishes). If we were studying polarization textures over an extended field, these $T_3 = 0$ orbits would serve to define the analogues of L-lines or L-surfaces. Note, however, that much of the published work on polarization singularities for bichromatic fields [8]–[11] associates to each frequency not a linear but a circular polarization.

## 3 Jones vector and classical geometric phase

We can write the position vector as $Q = Q_0 \text{Re}\{\text{diag}[\exp(-im_x t), \exp(-im_y t)]\mathbf{v}\}$, where the normalized complex vector $\mathbf{v}$ (analogous to the Jones vector for optical polarization when $m_x = m_y = 1$) is defined as

$$ \mathbf{v}(\theta, \phi) = \left[ \cos \frac{\theta}{2} \exp\left(-i \frac{\phi}{2m_y}\right), \sin \frac{\theta}{2} \exp\left(i \frac{\phi}{2m_x}\right) \right]. \quad (11) $$

In analogy with the isotropic case, we can define a geometric phase associated with transformations of $\mathbf{v}$. Let us consider continuous changes corresponding to a parametrization of $\theta$ and $\phi$ in terms of a continuous variable $\eta$. The rate of accumulation of geometric phase is then defined as the phase difference between two infinitesimally separated states of polarization:

$$ d\Phi_c = \arg\{\mathbf{v}^*[\theta(\eta), \phi(\eta)] \cdot \mathbf{v}[\theta(\eta + d\eta), \phi(\eta + d\eta)]\}. \quad (12) $$

By substituting the expression for $\mathbf{v}$ and simplifying, one gets

$$ \begin{aligned} d\Phi_c &= \left[ -\frac{1}{2m_y} \cos^2 \frac{\theta}{2} + \frac{1}{2m_x} \sin^2 \frac{\theta}{2} \right] \partial_{\eta} \phi \, d\eta \\ &= -\frac{m_x + m_y}{4m_x m_y} (\cos \theta - \epsilon) \partial_{\eta} \phi \, d\eta. \end{aligned} \quad (13) $$

Note that the factor $(\cos \theta - \epsilon) \partial_{\eta} \phi \, d\eta$ can also be written as $(\tau_1 - \epsilon) d\phi$, suggesting the following geometric interpretation [shown in Figure 2a]: Consider the projection from the Kummer shape onto a cylinder of unit radius centered at the $\tau_1$ axis. The accumulation of geometric phase is then proportional to the area enclosed between the projection of the equator and the projection of the curve traced by $\tau_1$. In the degenerate case $m_x = m_y = 1$, this area over a cylinder equals the area over the sphere, namely a solid angle, yielding the standard geometric phase for optical polarization. In the general non-degenerate case, on the other hand, the geometric phase does not correspond to an area or solid angle over the Kummer shape.

Note that the rate of accumulation of geometric phase described above is associated with the specific definition of the Jones vector in Eq. (11); different definitions would result from multiplying this one by an arbitrary phase factor that depends on $\theta$ and $\phi$. However, when we consider the total phase accumulated over a closed loop, the geometric phase becomes independent of the specific definition and is of a more fundamental nature. Let us consider only loops over the Kummer shape that do not cross themselves. There are two types of such loops:

i) Loops that do not surround the $\tau_1$ axis [illustrated in Figure 2b], for which the initial and final values of $\phi$ coincide. For these, the total geometric phase is

$$ \Phi_c = \frac{m_x + m_y}{4m_x m_y} \Omega, \quad (14) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2:** A geometric interpretation for the classical geometric phase. (a) The integral in $\eta$ of $(\tau_1 - \epsilon)\partial_\eta\phi d\eta$ over a given interval corresponds to the area (magenta) accumulated between the equator and the projection of the curve from the Kummer shape onto a cylinder of unit radius centered at the $\tau_1$ axis. (b) For a closed loop that does not</figcaption></figure>

where $\Omega$ is the area enclosed over the cylinder (or equivalently, over the unit sphere with spherical angles $\theta, \phi$), with a sign determined by the handedness of the loop (positive for a right-handed loop in the outward sense).

ii) Loops that surround the $\tau_1$ axis [illustrated in Figure 2c], where the initial and final values of $\phi$ differ by $2\pi$. Let us first assume that $\phi$ increases by $2\pi$ as we complete the loop. The area between the equator and the loop equals $2\pi(1 - \epsilon)$ minus the area above the curve, which can be regarded as the enclosed area $\Omega_+$. Therefore,

$$ \begin{aligned} \bar{\Phi}_c &= \frac{m_x + m_y}{4m_x m_y} [\Omega_+ - 2\pi(1 - \epsilon)] \\ &= \frac{m_x + m_y}{4m_x m_y} \Omega_+ - \frac{\pi}{m_y}. \end{aligned} \quad (15a) $$

If, on the other hand, $\phi$ decreases by $2\pi$ as the loop is closed (i.e., the sense is reversed) one gets

$$ \bar{\Phi}_c = \frac{m_x + m_y}{4m_x m_y} \Omega_- - \frac{\pi}{m_x}, \quad (15b) $$

where $\Omega_-$ corresponds to the cylinder's area below the curve. Notice that the extra contributions in Eqs. (15) occur even in the isotropic case, where they reflect the nonperiodicity of the Jones vector, namely $v(\theta, \phi + 2\pi) \neq v(\theta, \phi)$. That is, there is a part of this accumulated phase that is not an intrinsic geometric phase proportional to an enclosed area but a consequence of the parametrization (the bar in $\bar{\Phi}_c$ indicating this). For the isotropic oscillator $m_x = m_y = 1$, $v^*(\theta, \phi) \cdot v(\theta, \phi \pm 2\pi) = -1$ for any $\theta$ and $\phi$, and this extra phase of $-\pi$ corresponds exactly to the extra contribution in Eqs. (15), so that the intrinsic part of the geometric phase is still given by Eq. (14). In the anisotropic case, this nonperiodicity is more complex, since $v^*(\theta, \phi) \cdot v(\theta, \phi + 2\pi)$ in general has a magnitude different from unity, and its phase depends on $\theta$. These phases take simple forms at the poles, and these forms coincide exactly with the extra terms in Eqs. (15): $\arg[v^*(0, \phi) \cdot v(0, \phi + 2\pi)] = -\pi/m_y$ and $\arg[v^*(\pi, \phi) \cdot v(\pi, \phi - 2\pi)] = -\pi/m_x$.

<footnote>surround the $\tau_1$ axis, $\Omega$ corresponds to the area (green) enclosed by the projection of the loop onto the cylinder, where the sign is determined by the handedness of the loop. (c) For a closed loop that surrounds the $\tau_1$ axis, $\Omega_+$ and $\Omega_-$ correspond to the cylinder's area above (green) and below (pale blue) the projection of the loop onto the cylinder.</footnote>
## 4 Semiclassical wave state construction

We now consider eigenstates for this system in the quantum/wave regime. Asymptotic estimates to these eigenstates can be constructed based on the Lissajous orbits. In order to do so, we must consider two-parameter families of these trajectories, given that the system is two-dimensional. First, we introduce a parameter $\xi$ to denote a temporal shift in the orbit according to $t \to t + \xi$. Since we are constructing modes whose shape is invariant in $t$, it is sufficient to consider $t = 0$, so effectively we just replace the physical variable $t$ with the delay parameter $\xi$. The second parameter, called here $\eta$, is used to parametrize both $\theta$ and $\phi$. The construction will also require the classical action $S$, which is linked to the positions and momenta by the following equations:

$$ \partial_{\xi} S = \mathbf{P} \cdot \partial_{\xi} \mathbf{Q}, \quad \partial_{\eta} S = \mathbf{P} \cdot \partial_{\eta} \mathbf{Q}. \qquad (16) $$

In the particular context of classical paraxial optics, where $\xi, \eta$ identify a ray, $S$ represents the optical path length from a normal to the rays, and Eqs. (16) guarantee that the family of rays constitutes a normal congruence, that is, that the rays are locally perpendicular to a surface element. As shown in Appendix A, these equations have a solution of the form

$$ S(\xi, \eta) = \Sigma(\xi, \theta, \phi) + S_{GP}(\eta), \qquad (17) $$

where $\Sigma(\xi, \theta, \phi)$ is given by a fairly simple closed-form expression in which the only dependence on $\eta$ is within the parameters $\theta, \phi$, and $S_{GP}$ is related to the geometric phase:

$$ S_{GP}(\eta) = \frac{Q_0^2}{4m_x m_y} \int^{\eta} \cos \theta \partial_{\eta} \phi \, d\eta. \qquad (18) $$

We now seek to construct the wave eigenstates based on the orbits. To do this, we use what is known as a semiclassical method, which is a procedure to construct asymptotic estimates of wave fields or quantum wavefunctions by using only the information provided by rays or classical trajectories. While several approaches exist [21], we use one that consists on dressing each orbit with a Gaussian contribution [22], [23] and that does not suffer from singularities at caustics. This approach has been applied to a variety of self-similar Gaussian beams [24]–[27] (which correspond to the modes of the isotropic 2D harmonic oscillator), for which it is able to actually reconstruct exact wave fields from the ray description. It is shown in what follows that this exact

ray-based wave (or classical-based quantum) construction extends to the anisotropic case, so the results found are not simply asymptotic estimates but true modes. The wave field is expressed as a superposition of Gaussians centered at the ray coordinates $\mathbf{Q}$, each presenting a linear phase ramp proportional to the ray momentum $\mathbf{P}$, that is,

$$ \Psi(\mathbf{x}) \propto \iint \sqrt{\frac{\partial(\mathcal{P}\mathbf{Q} + i\mathbf{P})}{\partial(\xi, \eta)}} \exp \left\{ -k \frac{(\mathbf{x} - \mathbf{Q})\mathcal{P}(\mathbf{x} - \mathbf{Q})}{2} + ik[S + (\mathbf{x} - \mathbf{Q}) \cdot \mathbf{P}] \right\} d\xi d\eta, \qquad (19) $$

where $\mathcal{P}$ is a constant matrix that determines the extent of the Gaussian contributions, whose eigenvalues have positive real parts.

The semiclassical estimate in Eq. (19) incorporates naturally Maslov phase shifts in the phase of the square root of the Jacobian if care is taken when evaluating the square root so that the resulting complex function is continuous. For sufficiently large $k$ (namely, for eigenstates of sufficiently high order) the estimate is asymptotically insensitive to the choice of $\mathcal{P}$. However, it turns out that analytic results for parts of the integral can be obtained in this case by using the diagonal matrix $\mathcal{D} = \text{diag}(m_x, m_y)$, for which $\exp(-k\mathbf{x}\mathcal{D}\mathbf{x}/2)$ corresponds to the ground state of the system. The determinant then reduces to a product of a function of $\eta$ and an exponential in $\xi$:

$$ \frac{\partial(\mathcal{P}\mathbf{Q} + i\mathbf{P})}{\partial(\xi, \eta)} = \mathcal{D} \exp \left[ i \left( \frac{1}{m_y} - \frac{1}{m_x} \right) \frac{\phi}{2} \right] \times \exp[-i(m_x + m_y)\xi], \qquad (20) $$

where

$$ \mathcal{D} = \frac{Q_0^2}{4} \{2 \sin \theta \partial_{\eta} \phi - i [m_x + m_y + (m_x - m_y) \cos \theta] \partial_{\eta} \theta\}. \qquad (21) $$

The solution to each of the two integrals in Eq. (19) is considered separately in the next two sections.

## 5 Coherent states

Given that both $\mathcal{D}$ and $\mathcal{S}$ have closed form expressions for their dependence on $\xi$, the integral in $\xi$ for the semiclassical estimate can be solved first. Let us then rewrite Eq. (19) as

$$ \Psi(\mathbf{x}) \propto \int \sqrt{\mathcal{D}} \exp(i\Theta_{GP}) \mathcal{V}_n(\mathbf{x}, \theta, \phi) \, d\eta, \qquad (22) $$
where $\Theta_{GP}(\eta) = kS_{GP} + (1/m_y - 1/m_x)\phi/4$, and $\mathcal{V}_n(\mathbf{x}, \theta, \phi)$ is the wavefunction contribution associated with the orbit specified by $\theta$ and $\phi$, given by the following integral in $\xi$ over an interval of length $2\pi$:

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi) = \int \exp \left\{ -k \frac{(\mathbf{x} - \mathbf{Q}) \mathcal{I} (\mathbf{x} - \mathbf{Q})}{2} - \frac{m_x + m_y}{2} i \xi + i k \left[ \Sigma + (\mathbf{x} - \mathbf{Q}) \cdot \mathbf{P} \right] \right\} d\xi. \quad (23) $$

It is shown in Appendix B that the periodicity of the integer imposes the quantization condition

$$ Q_0 = \sqrt{\frac{2n + m_x + m_y}{k}}, \quad (24) $$

where $n$ is a non-negative integer. With this condition, the integral can be calculated in closed form, as shown in Appendix B, giving a result expressible in terms of the components of the Jones vector $\mathbf{v}$:

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi) = -2\pi \exp \left[ -\frac{k}{2}(m_x x^2 + m_y y^2) - \frac{kQ_0^2}{4} \left( \frac{|\nu_x|^2}{m_x} + \frac{|\nu_y|^2}{m_y} \right) \right] \times F_n \left( kQ_0 \nu_x x, kQ_0 \nu_y y, -\frac{kQ_0^2 \nu_x^2}{4m_x}, -\frac{kQ_0^2 \nu_y^2}{4m_y} \right), \quad (25) $$

where $F_n(\alpha_x, \alpha_y, \beta_x, \beta_y)$ is an nth-order four-variable polynomial composed of terms proportional to $\alpha_x^{\mu_x} \alpha_y^{\mu_y} \beta_x^{\nu_x} \beta_y^{\nu_y}$, for all combinations of non-negative integer powers for which $m_x(\mu_x + 2\nu_x) + m_y(\mu_y + 2\nu_y) = n$. This polynomial satisfies the recurrence relations (for $\rho = x, y$):

$$ \partial_{\alpha_{\rho}} F_n = F_{n-m_{\rho}}, \quad \partial_{\beta_{\rho}} F_n = F_{n-2m_{\rho}}, \quad (26) $$

$$ F_n = \sum_{\rho=x,y} m_{\rho} \frac{\alpha_{\rho} F_{n-m_{\rho}} + 2\beta_{\rho} F_{n-2m_{\rho}}}{n}, \quad (27) $$

with $F_{n<0} = 0$ and $F_0 = 1$. The relation in Eq. (27) is particularly convenient for the fast and numerically robust computation of these polynomials when considering high orders [28], the computation time being proportional to $n$.

The solutions $\mathcal{V}_n$ just found are not simply asymptotic estimates but exact eigenstates of the Hamiltonian. This is shown by using the recurrence relations to find the constant eigenvalue:

$$ \begin{aligned} \frac{\hat{H} \mathcal{V}_n}{\mathcal{V}_n} &= \frac{[-k^{-2}(\partial_x^2 + \partial_y^2) + m_x^2 x^2 + m_y^2 y^2] \mathcal{V}_n}{2 \mathcal{V}_n} \\ &= \frac{m_x + m_y}{2k} + \sum_{\rho=x,y} m_{\rho} \frac{\alpha_{\rho} \partial_{\alpha_{\rho}} F_n + 2 \beta_{\rho} \partial_{\alpha_{\rho}}^2 F_n}{k F_n} \\ &= \frac{m_x + m_y}{2k} + \sum_{\rho=x,y} m_{\rho} \frac{\alpha_{\rho} F_{n-m_{\rho}} + 2 \beta_{\rho} F_{n-2m_{\rho}}}{k F_n} \\ &= \frac{2n + m_x + m_y}{2k}. \end{aligned} \quad (28) $$

This result holds for any $\theta$ and $\phi$, so all states with equal $n$ are degenerate and can be combined to construct other eigenstates with the same eigenvalue. The states $\mathcal{V}_n$ can be regarded as coherent states, and their wave profile mimics the Lissajous shape of the orbit for the corresponding values of $\theta$ and $\phi$, as illustrated in Figure 3. Some of these states have been produced and measured experimentally as cavity modes [3], [4]. Expressions for these coherent states have been given as specific combinations involving Wigner d coefficients [2], [29]–[31]. The different approach taken here where the orbit is dressed with appropriate Gaussians gives an expression in terms of a polynomial that can be computed directly in a numerically robust way by using the recurrence relation in Eq. (27).

The recurrence relations given above show that, following the ground state $n=0$, the first nonzero eigenstate corresponds to $n=\min(m_x, m_y)$, and eigenstates exist only for values of $n$ that are nonnegative integer combinations of $m_x$ and $m_y$. Incidentally, note that for the isotropic case $m_x = m_y (=1)$, this recurrence relation reduces to one with two (rather than four) terms on the right-hand side, and can be transformed into that for Hermite polynomials. This is not surprising since the coherent states for the isotropic oscillator are expressible in terms of Hermite polynomials of a complex expression involving the inner product of the position and Jones vectors [24], [32]. In the context of bichromatic polarization, $n$ is related to the number of photons in the state, but it is not necessarily equal to it: the state is composed of all possible numbers of photons $i_x + i_y$ for which $i_x m_x + i_y m_y = n$. In fact, for $\theta \neq 0, \pi$, the coherent states $\mathcal{V}_n$ are in general entangled in the two frequency components (i.e., non-separable in $x$ and $y$), except in special cases where $n$ is only composed as an integer multiple of one of the frequencies $m_\rho$ without involving the other.

Finally, note that the recurrence relations in Eqs. (26) and (27) can be combined to arrive at the following relation for the coherent states themselves:
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Coherent modes associated with a Lissajous orbit. (a) Modulus squared and (b) real part of a coherent mode with $n = 200$, $m_x = 2$, $m_y = 3$, $\theta = \pi/2$ and $\phi = 3\pi/8$. The classical orbit is overlaid on both parts as a black curve. In (b), the inset shows a detail of the mode's structure near a crossing of the orbit with itself.</figcaption></figure>

$$2nV_n = \sum_{\rho=x,y} Q_0 v_{\rho} (km_{\rho}\rho - \partial_{\rho}) V_{n-m_{\rho}}. \quad (29)$$

The operators acting on the coherent states on the right hand side are clearly creation operators on each of the Cartesian variables, which cause an increase in the total index $n$ not of unity but of $m_x$ or $m_y$. We can then see that, to get to a given order $n$, we cannot simply apply the same operator repeatedly, but we have to consider all the combinations of the operators for each of the cartesian variables/frequencies that cause the state to arrive at that order.

# 6 Semiclassical quantization of geometric phase

The calculation of a generic eigenstate in Eq. (22) involves a weighted superposition of coherent states. Of course, any

such superposition is an eigenstate, but this semiclassical formula is constructed such that neighboring coherent states are superposed asymptotically in phase along the caustics, so that the resulting waveform mimics the classic structure. If we consider a closed integral such that the integrand is a periodic function of $\eta$, this requirement of constructive interference of contiguous coherent states implies the accumulation of a phase that must be quantized to guarantee phase self-consistency when closing the integral.

We again consider the closed (non-self-crossing) path over the Kummer shape traced by $\theta(\eta)$, $\phi(\eta)$. As discussed in Section 3, we consider separately the two types of loop:

i) Loops that do not surround the $\tau_1$ axis, for which the initial and final values of $\phi$ coincide. Notice that $D$ accumulates a phase of $\pm 2\pi$ upon closing the loop (the sign depending on the loop's sense), and hence $\sqrt{D}$ accumulates a phase of $\pm \pi$. The periodicity of the integrand then requires that $\Theta_{GP}$ and therefore $kS_{GP}$ accumulate an odd multiple of $\pi$, so that the area $\Omega$ enclosed by the loop projected to the cylinder is discretized according to

$$\Omega = - \int \cos \theta \, \partial_{\eta} \phi \, d\eta = \frac{4m_x m_y \pi}{2n + m_x + m_y} (2\mu + 1), \quad (30)$$

for integer $\mu$.

ii) Loops that surround the $\tau_1$ axis, where the initial and final values of $\phi$ differ by $2\pi$. In this case, $D$ does not accumulate a phase. However, an increase of $2\pi$ in $\phi$ causes $V_n$ to accumulate a phase $\Phi_{V_n}$, shown in Appendix B to be given by

$$\Phi_{V_n} = \frac{i_y m_y - i_x m_x}{m_x m_y} \pi \quad (31)$$

with $i_x$ and $i_y$ being non-negative integers such that $i_x m_x + i_y m_y = n$. Note that the choice of $i_x$ and $i_y$ is typically not unique: if $n$ is sufficiently large, one can change $i_x \to i_x - \nu m_y$, $i_y \to i_y + \nu m_x$ for any integer $\nu$ that allows keeping these two numbers non-negative; however, this change simply leads to $\Phi_{V_n} \to \Phi_{V_n} + 2\pi\nu$, so the specific choice of the integers $i_x$ and $i_y$ (within their allowed values) is inconsequential. The total phase accumulated by the integrand in Eq. (22), given by

$$\Theta_{GP} + \Phi_{V_n} = \frac{2n + m_x + m_y}{4m_x m_y} \int \cos \theta \, \partial_{\eta} \phi \, d\eta + \left( \frac{1}{m_y} - \frac{1}{m_x} \right) \frac{\pi}{2} + \Phi_{V_n}, \quad (32)$$

must be an integer multiple of $2\pi$. Let us denote by $\Omega_+$ the cylinder area above the loop:
$$ \Omega_+ = \int (1 - \cos \theta) \partial_\eta \phi \, d\eta. \qquad (33) $$

This area is then quantized as

$$ \Omega_+ = \frac{4m_x m_y \pi}{2n + m_x + m_y} \left( 2\mu + \frac{2i_y}{m_x} + \frac{1}{m_y} \right), \quad (34a) $$

where in order to arrive to the last step we used the relation $n = i_x m_x + i_y m_y$. We can, without loss of generality, choose the smallest possible allowed value of $i_y$. Note that, had we chosen instead to use the complementary area below the curve $\Omega_- = 4\pi - \Omega_+$ we would get a similar result but with the subindices x and y swapped:

$$ \Omega_- = \frac{4m_x m_y \pi}{2n + m_x + m_y} \left( 2\mu + \frac{2i_x}{m_y} + \frac{1}{m_x} \right). \quad (34b) $$

The quantization conditions for the areas enclosed by the curves in the two cases, given, respectively, in Eqs. (30) and (34), resemble each other except for the offsets inside the parentheses. Notice that Eq. (34a) or Eq. (34b) can only fully coincide with Eq. (30) if $m_y$ or $m_x$ equal unity, respectively, and $n$ is an integer multiple of both $m_x$ and $m_y$ (since then $i_y$ or $i_x$ can be chosen to be zero). This suggests that, for cases where $n$ is an integer multiple of $m_x m_y$, the difference in the quantization condition has a topological origin and is related to whether the Kummer shape's pole associated with the area under consideration [namely, $\tau = (\pm 1, 0, 0)$ for $\Omega_\pm$] is smooth or not; if it is smooth, the quantization is not different from the one for a loop not surrounding the $\tau_1$ axis.

## 7 Coherent state overlap, norm, and Husimi representation

We now calculate the overlap between any two coherent states $\langle U_n^{(1)} | U_n^{(2)} \rangle$, where the superindices indicate different parameters $\theta_1, \phi_1$ and $\theta_2, \phi_2$ for each. By using the definition of the polynomials $F_n$ in terms of derivatives of an exponential in Eq. (57), the overlap integrals can be solved in closed form using standard Gaussian formulas:

$$ \begin{aligned} \langle U_n^{(1)} | U_n^{(2)} \rangle &= \iint U_n(\mathbf{x}, \theta_1, \phi_1)^* U_n(\mathbf{x}, \theta_2, \phi_2) \, dx \, dy \\ &= \frac{4\pi^3}{k\sqrt{m_x m_y}} \exp\left(-\frac{2n + m_x + m_y}{4}\right) \\ & \quad \times \sum_{\rho=x,y} \frac{|v_{1\rho}|^2 + |v_{2\rho}|^2}{m_\rho} \end{aligned} $$

$$ \times G_n \left( \frac{2n + m_x + m_y}{2} v_{1x}^* v_{2x}, \frac{2n + m_x + m_y}{2} v_{1y}^* v_{2y} \right), \quad (35) $$

where

$$ G_n(w_x, w_y) = \left. \left\{ \frac{\partial_{u_1}^n \partial_{u_2}^n}{(n!)^2} \exp \left[ \sum_{\rho=x,y} \frac{(u_1 u_2)^{m_\rho} w_\rho}{m_\rho} \right] \right\} \right|_{\substack{u_1=u_2=0 \\ u_1=u_2=1}} \quad (36) $$

is a polynomial with two arguments, which can be shown to satisfy the recurrence formula

$$ G_n = \frac{1}{n^2} \left( \sum_{\rho=x,y} m_\rho w_\rho G_{n-m_\rho} + \sum_{\rho_1=x,y} \sum_{\rho_2=x,y} w_{\rho_1} w_{\rho_2} G_{n-m_{\rho_1}-m_{\rho_2}} \right), \quad (37) $$

with $G_{n<0} = 0$ and $G_0 = 1$. Like $F_n$, this polynomial can be computed rapidly and in a numerically stable form using this type of recurrence relation.

The overlap formula in Eq. (35) allows calculating the norm of the coherent states $U_n$. By evaluating this equation for $U_n^{(1)} = U_n^{(2)} = U_n$ (corresponding to the angles $\theta, \phi$), one can verify that the result is independent of $\phi$, but that in general it depends on $\theta$ (except in the isotropic case $m_x = m_y = 1$). In fact, the norm of the coherent states tends to zero in the limit $\theta \to 0(\pi)$ unless $n$ is an integer multiple of $m_x(m_y)$. This is easy to understand, since in this limit the oscillator involves a single frequency, and $n$ must be an integer multiple of this frequency for the solution to exist.

The coherent states $U_n$ can be used to construct a Husimi distribution representation of arbitrary modes $\Psi$ with the same $n$ according to

$$ H(\theta, \phi) = \frac{|\langle U_n(\theta, \phi) | \Psi \rangle|^2}{\langle U_n(\theta, \phi) | U_n(\theta, \phi) \rangle}, \quad (38) $$

similar to that used for the isotropic oscillator, for characterizing quantum states of polarization [33] or structured Gaussian beams [34]. This distribution can be represented over the surface of the corresponding Kummer shape. Figure 4 shows this distribution for the case $m_x = 2, m_y = 3$ for a wavefunction $\Psi$ proportional to a coherent state $U_n$ with $\phi = 0$ and for both $\theta \approx 0$ (a,b) and $\theta = \pi/2$ (c,d). In all cases, the maximum is located around the point that represents the corresponding orbit. Note, however, that the coherent state $\theta = 0$ is not allowed except when $n$ is an integer multiple of $m_x$, which is the case in (a,c) corresponding to $n = 204$, but not in (b,d) corresponding to $n = 205$. For this reason, a state was used in (a,c) with $\theta = 10^{-5}$ rather than
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Husimi representation for $m_x = 2$, $m_y = 3$ for coherent states with $\phi = 0$, for (a,b) $\theta = 10^{-5} \approx 0$ and (c,d) $\theta = \pi/2$, and where (a,c) $n = 204$ (an integer multiple of $m_x$ and $m_y$), and (b,d) $n = 205$ (not an integer multiple of $m_x$ or $m_y$).</figcaption></figure>

one corresponding exactly to the pole. Similarly, when evaluating Eq. (38) when $n$ is not an integer multiple of $m_x(m_y)$, the Husimi distribution cannot be evaluated exactly at $\theta = 0(\pi)$. This irregular behavior at the pole makes the maximum of the Husimi function appreciably more extended in (b) than in (a), despite $n$ being marginally larger for the latter. In contrast, the maxima in (c) and (d) are essentially identical (the latter being in theory very slightly more localized due to the larger value of $n$).

When applied to more complex modes, the Husimi distributions become more extended and can present zeros at a number of points. In the symmetric case ($m_x = m_y = 1$ the location of these zeros is sufficient to characterize the state (up to a global factor, of course) and is known as the Majorana constellation [33], [34]. The number of zeros in the commensurate anisotropic case is $N$ if $n = Nm_xm_y$.

## 8 Quantum geometric phase

We now consider the rate of accumulation of a “quantum” geometric phase by the coherent states as their parameters trace a closed trajectory in $\theta$ and $\phi$. This geometric phase, defined as

$$ \Phi_Q = \oint \arg\{\langle V_n^*[\theta(\eta), \phi(\eta)] | V_n^*[\theta(\eta + d\eta), \phi(\eta + d\eta)] \rangle\}, \quad (39) $$

can be calculated using the overlap expression in Eq. (35). Note that the factors multiplying $G_n$ in Eq. (35) are real and positive, so they do not contribute to the phase. We start by considering the phase accumulation due to infinitesimal changes in $\theta$ and $\phi$ given by $d\theta = \partial_\eta \theta d\eta$ and $d\phi = \partial_\eta \phi d\eta$, respectively. We then evaluate the polynomials at $\omega_\rho = (kQ_0^2/2)\nu_\rho^*|_{\theta,\phi}\nu_\rho|_{\theta+d\theta,\phi+d\phi}$, i.e.,

$$ w_x = \frac{2n + m_x + m_y}{2} \cos^2 \frac{\theta}{2} \left( 1 - i \frac{d\phi}{2m_y} \right), \quad (40a) $$

$$ w_y = \frac{2n + m_x + m_y}{2} \sin^2 \frac{\theta}{2} \left( 1 + i \frac{d\phi}{2m_x} \right). \quad (40b) $$

That is, like for the classical geometric phase, only differentials in $\phi$ play a role, so it suffices to compare the rate of phase accumulation in $\phi$ as a function of $\theta$. Recall that, for the classical geometric phase, this rate equals $-(m_x + m_y)/(4m_xm_y)(\tau_1 - \epsilon)$.

It is useful to consider first the special case of the isotropic oscillator, corresponding to $m_x = m_y = 1$. One can see that the recurrence relation in Eq. (37) becomes then a two-term relation given by $G_n = (WG_{n-1} + W^2G_{n-2})/n^2$ with $W = w_x + w_y = (n+1)\exp(-i\tau_1 d\phi/2)$ for infinitesimal changes, following Eqs. (8). This recursion relation has a solution of the form $G_n \propto W^n$, so the rate of phase accumulation in $\phi$ by $G_n$ is exactly $-n/2 \cos \theta$, and hence the exact quantum geometric phase is simply $n$ times the classical one: $\Phi_Q = n\Phi_C$.

For anisotropic oscillators, on the other hand, the rate of phase accumulation takes a more complicated form that typically requires numerical calculation. There are, however, limiting situations where this rate can be calculated analytically: at the poles, one of the two expressions in Eqs. (8) vanishes, reducing the recurrence relation in Eq. (37) to $G_n = (m_\rho w_\rho G_{n-m_\rho} + w_\rho^2 G_{n-2m_\rho})/n^2$, where $\rho = x$ for $\theta = 0$ and $\rho = y$ for $\theta = \pi$. If $n$ is an integer multiple of $m_\rho$, we find $G_n \propto w_\rho^{n/m_\rho}$, giving a rate of phase accumulation of $\mp n/(2m_xm_y)$, the sign corresponding to $\theta$ being 0 or $\pi$. This suggests the following asymptotic relation for the rate of quantum geometric phase accumulation with $\phi$:

$$ \partial_{\phi} \Phi_Q \approx -\frac{n}{2m_x m_y} \tau_1. \quad (41) $$

That is, ignoring the offset of $\epsilon$ in the classical geometric phase (which is inconsequential for loops not surrounding the $\tau_1$ axis), the quantum geometric phase is approximately $(m_x + m_y)n/2$ times larger than the classical one for sufficiently large $n$. (Note that this approximate result indeed reduces to $\Phi_Q = n\Phi_C$ for the isotropic case $m_x = m_y = 1$.)

Figure 5 shows numerically calculated plots of $-2m_xm_y\partial_{\phi}\Phi_Q/n$ versus $\tau_1$, for $m_x = 2$, $m_y = 3$, and for
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Plots of $-2m_x m_y \partial_\phi \Phi_Q / n$ versus $\tau_1$ for $m_x = 2, m_y = 3$. The values of $n$ are (a) 12 (red), 102 (green) and 1002 (dashed blue), and (b) 10 (red), 100 (green) and 1000 (dashed blue).</figcaption></figure>

increasing values of $n$ that are (a) integer multiples of $m_x$ and $m_y$, and (b) integer multiples of $m_x$ but not of $m_y$. In both cases, the approximate result in Eq. (41) is validated for sufficiently large $n$, but in the latter case there is a noticeable departure from linearity within a region around $\tau_1 = -1$. This discrepancy can be understood by the fact that the relation $G_n \propto \omega_\rho^{n/m_\rho}$ mentioned earlier for $\theta = 0, \pi$ is actually not valid if $n/m_\rho$ is not an integer. That is, no

quantized state exists that involves only one frequency $m_\rho$ if $n$ is not an integer multiple of that frequency. To obtain the behavior in the close vicinity of this forbidden value, we must consider $\theta$ being close to its endpoints but not equal to them, so that one of the two expressions in Eqs. (8) is much smaller than the other but not zero. The dominant contribution to $G_n$ will then be approximately proportional to $\omega_\rho^i \omega_{\bar{\rho}}^{\bar{i}}$, with $\rho = x, y, \bar{\rho} = y, x$ and where the nonnegative integers $i_\rho$ and $i_{\bar{\rho}}$ satisfy the relation $i_\rho m_\rho + i_{\bar{\rho}} m_{\bar{\rho}} = n$ with $i_{\bar{\rho}}$ being as small as possible. Therefore, the correct result near $|\tau_1| = 1$ is $-2m_x m_y \partial_\phi \Phi_Q / n \approx m_x i_x - m y i_y$, for $i_x$ as small as possible if $\tau_1 \approx 1$ and $i_y$ as small as possible if $\tau_1 \approx -1$. This explains why in Figure 5b, for which $n$ is not an integer multiple of $m_y$, the curves start at $(2m_x - 2m_y)/n = -0.2$ for $n = 10$, $(2m_x - 32m_y)/n = -0.92$ for $n = 100$, and $(2m_x - 332m_y)/n = -0.992$ for $n = 1000$.

## 9 Geometric phase for the comparison of states

So far we have consider the accumulation of geometric phase due to continuous transformations of the orbit shapes. For the isotropic oscillator, these continuous transformations are easily realizable: one can consider, for example, the passage of polarized light through birefringent media with arbitrary eigenpolarizations as a way to transform continuously the polarization state of a monochromatic beam. For commensurate anisotropic oscillators, on the other hand, the situation is different. Uniform shifts in $\phi$ (corresponding to rotations around $\tau_1$, which is the only axis of rotational symmetry of the Kummer shape) are easy to implement, through birefringence for bichromatic polarization or by detuning the frequency ratio for the cavity. On the other hand, operations that change $\theta$ are not easy to envisage, especially if we want these operations to be unitary. It can therefore be challenging to implement experiments to observe the accumulation of geometric phase described so far.

There is, however, a simple manifestation of geometric phase that is not based on continuous transformations but on the comparison of the phases of at least three states. Suppose that we have three different states, labelled as A, B, and C. We start by measuring the interference between A and B, and we adjust the global phase of B so that it is “in phase” with A, meaning that the constructive interference of A and B is maximized. We then do the same between B and C, giving C the necessary phase to maximize its interference with B. Despite having optimized the phase agreement between A and B and between B and C, in general
the phase agreement between C and A is not optimal, in that an extra phase would be needed to make their interference as constructive as possible. This difference is due to the geometric phase. For the isotropic oscillator (e.g. for standard monochromatic polarization), this phase is given by half the area over the sphere enclosed by the geodesic triangle defined by the three states. We now consider the general anisotropic case.

Let us start by considering the classical geometric phase, defined as

$$ \Phi_C = \arg[(\mathbf{v}_A^* \cdot \mathbf{v}_B)(\mathbf{v}_B^* \cdot \mathbf{v}_C)(\mathbf{v}_C^* \cdot \mathbf{v}_A)]. \quad (42) $$

A geometric interpretation for this phase follows from rewriting the Jones vector in Eq. (11) as

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Quantum geometric phases between three states. (a) Quantum geometric phase between the states with angles $(\theta, \phi)$ given by $(\pi/2, 0)$, $(\pi/2, \pi/4)$, and $(\pi/4, 0)$, as a function of $n$ and for $m_x = 2, m_y = 3$ (blue), and $m_x = 3, m_y = 2$ (orange). The gray line corresponds to the corresponding estimate in Eq. (44). (b) Relative error between the quantum geometric phase and the asymptotic estimate for $n = 1000$, $m_x = 2, m_y = 3$, and for the same states as in (a) except that the third state is replaced by $(\theta_c, 0)$ with $\theta_c$ varying from 0 to $\pi$.</figcaption></figure>

$$ \begin{aligned} \mathbf{v}(\theta, \phi) = & \exp\left(i \frac{m_y - m_x}{4m_x m_y} \phi\right) \\ & \times \left[ \cos \frac{\theta}{2} \exp\left(-i \frac{m_y + m_x}{4m_x m_y} \phi\right), \right. \\ & \left. \sin \frac{\theta}{2} \exp\left(i \frac{m_y + m_x}{4m_x m_y} \phi\right) \right]. \end{aligned} \quad (43) $$

Note that the global phase factor at the front has no effect on $\Phi_C$, since it cancels for each state in Eq. (42). The second, vectorial part is mathematically identical to a standard Jones vector for a single frequency, except that $\phi$ is scaled by a rational factor. Therefore, a geometric interpretation for $\Phi_C$ corresponds to half the area enclosed inside the geodesic triangle (or polygon, for more than three states) over a sphere (or more precisely, a spherical section) whose azimuthal angle is not $\phi$ but $\frac{m_y+m_x}{4m_x m_y} \phi$. Note, however, that geodesics over this sphere do not correspond to geodesics over the projected cylinder nor over the Kummer shape. In particular, for a given set of states corresponding to three sets of angles $\theta$ and $\phi$, the classical geometric phase does not change if we exchange the values of $m_x$ and $m_y$.

We now consider the quantum geometric phase, given by

$$ \Phi_Q = \arg[\langle \mathcal{V}_n^{(A)} | \mathcal{V}_n^{(B)} \rangle \langle \mathcal{V}_n^{(B)} | \mathcal{V}_n^{(C)} \rangle \langle \mathcal{V}_n^{(C)} | \mathcal{V}_n^{(A)} \rangle]. \quad (44) $$

In general, this phase has no simple closed form expression or geometric interpretation, but it can be easily calculated by using Eq. (35). Like for the continuous case, for $n \gg m_x, m_y$ we expect the asymptotic relation

$$ \Phi_Q \approx \frac{m_x + m_y}{2} n \Phi_C. \quad (45) $$

Figure 6a considers the geometric phases for all $n$ up to 100 corresponding to three states with $(\theta, \phi)$ given by $(\pi/2, 0)$, $(\pi/2, \pi/4)$, and $(\pi/4, 0)$, for two anisotropic oscillators with $m_x = 2, m_y = 3$ and $m_x = 3, m_y = 2$. We observe that, for small $n$ the quantum geometric phases differ significantly from the asymptotic estimate and from each other, but later they grow approximately linearly in $n$, mimicking Eq. (44) but with a slope that is slightly higher or lower depending on whether $m_x$ or $m_y$ is larger, somewhat reflecting the Kummer shape's asymmetry in contrast to the classical geometric phase. To understand this asymmetry, Figure 6b shows the difference between $\Phi_Q$ and the estimate in Eq. (44), normalized by $\Phi_Q$ for the states with $(\theta, \phi)$ given by $(\pi/2, 0)$, $(\pi/2, \pi/4)$, and $(\theta_c, 0)$ as a function of $\theta_c$ for $n = 1000$. The maximum error in slope is of about
3–4 %, but for $n = 1000$ this means a phase difference of about $\pi/2$.

## 10 Concluding remarks

The semiclassical treatment provided here for the two-dimensional commensurate anisotropic harmonic oscillator leads to two main results. The first is a robust exact construction of the coherent states associated with the orbits, in which each “particle” is dressed by an appropriate Gaussian wave packet. Remarkably, the integral can be solved in closed form and the resulting polynomials satisfy simple recurrence relations that facilitate their computation. These coherent states can be used for constructing more complicated modes, as well as for defining Husimi representations and Majorana constellations. Appendix C provides a simple implementation for Wolfram’s Mathematica of this recurrence relation, as well as for the one that gives access to the overlap between two coherent states.

The second result is the definition of both classical and quantum geometric phases for these systems. In the semiclassical limit (corresponding to highly excited states, for which $n \gg m_x, m_y$), the quantum geometric phase is proportional to the classical one, the proportionality constant being $(m_x + m_y)n/2$. Like in the case of the isotropic oscillator, these phases are proportional to enclosed areas, but these do not correspond to the area of the curved manifold that naturally represents the constants of the motion for the system (the Kummer shape). When we are considering continuous transformations, the area can be associated with that of a projection from the $\tau_1$ axis onto a cylinder of unit radius. (For the isotropic oscillator this projected area happens to equal that of the sphere.) When we consider instead the phase between a discrete set of states, the geometric interpretation requires considering geodesic polygons over a sphere where the azimuthal angle is scaled.

These results can also be applied to commensurate oscillators in which the oscillations of each frequency are not along a line but trace a more general elliptic shape. One such example are bichromatic superpositions of counterrotating circular polarizations [3], [8]–[10], [16], [35]. A general class of these results can be obtained by applying a linear canonical transformation (or Collins’ ABCD formula) [36], such as an astigmatic fractional Fourier transform whose axes are not aligned with x and y. This transformation is unitary, so the results for the classical and quantum geometric phases still apply. Further, because these integral transformations involve kernels that are exponentials of imaginary quadratic polynomials, they can be performed in closed form when applied to the expressions in Appendix B.

Finally, note that the semiclassical approach used here can also be applied to anisotropic commensurate harmonic oscillators in three or more dimensions, probably leading also to closed forms for the coherent states and to asymptotic expressions for the geometric phase. It will be interesting to see if in the three-dimensional case this phase accepts a geometric interpretation similar to one proposed recently for isotropic 3D oscillators [37].

**Acknowledgment:** The author thanks Miguel A. Bandres, Luis L. Sánchez-Soto, Mark Dennis, Konstantin Bliokh, Rodrigo Gutierrez-Cuevas, Emilio Pisanty, Mike Raymer and Greg Forbes for useful discussions.

**Research funding:** This research received funding from the Agence Nationale de la Recherche (ANR) through the project 3DPol, ANR-21-CE24-0014-01.

**Author contributions:** The author confirms the sole responsibility for the conception of the study, presented results and manuscript preparation.

**Conflict of interest:** Author states no conflicts of interest.

**Data availability:** Data sharing does not apply to this article as no datasets were generated or analysed during the current study.

## Appendix A: Equations for the action

We can solve Eqs. (16) by first setting $\xi$ to a constant such as zero and integrating in $\eta$ the equation for $\partial_\eta S$ to find a function $S_0(\eta)$, and then adding to this solution the integration in $\xi$ (from zero) of the equation for $\partial_\xi S$ for fixed $\eta$ (and hence fixed $\theta$ and $\phi$). That is, we write $S(\xi, \eta) = S_0(\eta) + \bar{S}(\xi, \theta, \phi)$, where each of the two parts satisfies one of the following equations:

$$ \partial_{\eta} S_0 = \frac{Q_0^2}{8m_x m_y} \left[ \left( m_x \sin \frac{\phi}{m_x} + m_y \sin \frac{\phi}{m_y} \right) \sin \theta \partial_{\eta} \theta + 2 \left( \cos \theta + \sin^2 \frac{\theta}{2} \cos \frac{\phi}{m_x} - \cos^2 \frac{\theta}{2} \cos \frac{\phi}{m_y} \right) \partial_{\eta} \phi \right], \quad (46) $$

and

$$ \partial_{\xi} S = Q_0^2 \left[ \cos^2 \frac{\theta}{2} \sin^2 \left( m_x \xi + \frac{\phi}{2m_y} \right) + \sin^2 \frac{\theta}{2} \sin^2 \left( m_y \xi - \frac{\phi}{2m_x} \right) \right]. \quad (47) $$

This latter equation can be integrated to give a closed-form expression of the form
$$ \bar{S}(\xi, \theta, \phi) = \frac{Q_0^2}{2} \left[ \xi - \frac{\sin(m_x \xi)}{m_x} \cos\left(m_x \xi + \frac{\phi}{m_y}\right) \cos^2 \frac{\theta}{2} - \frac{\sin(m_y \xi)}{m_y} \cos\left(m_y \xi - \frac{\phi}{m_x}\right) \sin^2 \frac{\theta}{2} \right]. \quad (48) $$

Equation (46) does not have a generic closed-form primitive, but it can be separated into a part that is a total derivative and one proportional to the derivative of the classical geometric phase:

$$ \partial_{\eta} S_{0} = -\partial_{\eta} \left\{ \frac{Q_{0}^{2}}{4} \left[ \frac{\operatorname{Im}(v_{x}^{2})}{m_{x}} + \frac{\operatorname{Im}(v_{y}^{2})}{m_{y}} \right] \right\} \\ + \frac{Q_{0}^{2}}{4m_{x}m_{y}} \cos \theta \partial_{\eta} \phi. \quad (49) $$

These results lead to the expression in Eq. (17), where the integrable part is given by

$$ \Sigma(\xi, \theta, \phi) = \bar{S}(\xi, \theta, \phi) - \frac{Q_0^2}{4} \left[ \frac{\operatorname{Im}(v_x^2)}{m_x} + \frac{\operatorname{Im}(v_y^2)}{m_y} \right]. \quad (50) $$

## Appendix B: Derivation of the coherent states and their recurrence relations

The wavefunction contribution $\mathcal{V}_n(\mathbf{x}, \theta, \phi)$ defined in Eq. (23) can be rewritten as

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi) = \exp\left(-k \frac{m_x x^2 + m_y y^2}{2}\right) \\ \times \int \exp\left\{k\mathbf{x} \cdot (\mathbb{Q} + i\mathbf{P}) - \frac{m_x + m_y}{2}i\xi\right\} \\ + i\kappa[\Sigma - \mathbf{Q} \cdot \mathbf{P}] \, d\xi, \quad (51) $$

the integral being over an interval of length $2\pi$.

By substituting $u = \exp(-i\xi)$ (so that $d\xi = iu^{-1}du$), the integral can be changed to one along the unit circle over the complex $u$ plane:

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi) = i \exp\left[-\frac{k}{2}(m_x x^2 + m_y y^2) - \frac{kQ_0^2}{4}\left(\frac{|v_x|^2}{m_x} + \frac{|v_y|^2}{m_y}\right)\right] \\ \times \oint \frac{\exp(\alpha_x u^{m_x} + \alpha_y u^{m_y} + \beta_x u^{2m_x} + \beta_y u^{2m_y})}{u^{n+1}} \, du, \quad (52) $$

with $\mathbf{v} = (v_x, v_y)$ defined in Eq. (11), and where

$$ n = \frac{kQ_0^2 - (m_x + m_y)}{2}, \quad (53) $$

$$ \alpha_x = kQ_0 v_x x, \quad \alpha_y = kQ_0 v_y y, \quad (54) $$

$$ \beta_x = -k \frac{(Q_0 v_x)^2}{4m_x}, \quad \beta_y = -k \frac{(Q_0 v_y)^2}{4m_y}. \quad (55) $$

Note that the integral in $u$ over a unit circle in the complex plane is uniquely defined only if the integrand is fully periodic, that is, if the power $-n-1$ of $u$ in the first factor is an integer. This is true if the quantization condition in Eq. (24) is satisfied with $n$ a nonnegative integer. Further, the remarkable fact that the exponential includes only a few powers of $u$, all of them positive, means that the integral can be evaluated in closed form by using the residue theorem, since the only singularity is a pole of order $n+1$ at the origin. The result can then be written as

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi) = -2\pi \exp\left[-\frac{k}{2}(m_x x^2 + m_y y^2) - \frac{kQ_0^2}{4}\left(\frac{|v_x|^2}{m_x} + \frac{|v_y|^2}{m_y}\right)\right] \\ \times F_n\left(kQ_0 v_x x, kQ_0 v_y y, -\frac{kQ_0^2 v_x^2}{4m_x}, -\frac{kQ_0^2 v_y^2}{4m_y}\right), \quad (56) $$

where $F_n$ is an $n$th-order four-variable polynomial defined as

$$ \begin{aligned} F_n(\alpha_x, \alpha_y, \beta_x, \beta_y) \\ &= \frac{1}{2\pi i} \oint \frac{\exp(\alpha_x u^{m_x} + \alpha_y u^{m_y} + \beta_x u^{2m_x} + \beta_y u^{2m_y})}{u^{n+1}} \, du \\ &= \frac{1}{n!} \left[ \partial_u^n \exp(\alpha_x u^{m_x} + \alpha_y u^{m_y} + \beta_x u^{2m_x} + \beta_y u^{2m_y}) \right]_{u=0}. \end{aligned} \quad (57) $$

It is easy to show that $F_n$ satisfies the recurrence relations in Eqs. (26) and (27), the proof for the latter consisting on using integration by parts in the integral form in Eq. (57).

It turns out that the coherent state $\mathcal{V}_n$ is not periodic in $\phi$, but it accumulates a phase factor if $\phi$ increases by $2\pi$:

$$ \mathcal{V}_n(\mathbf{x}, \theta, \phi + 2\pi) = \mathcal{V}_n(\mathbf{x}, \theta, \phi) \exp(i\Phi_{\mathcal{V}_n}). \quad (58) $$

Note that the factors in front of $F_n$ in the definition of $\mathcal{V}_n$ are real, but $F_n$ itself is complex. As mentioned earlier, this polynomial is composed of terms proportional to $\alpha_x^{\mu_x} \alpha_y^{\mu_y} \beta_x^{\nu_x} \beta_y^{\nu_y}$, for non-negative integer combinations $m_x(\mu_x + 2\nu_x) + m_y(\mu_y + 2\nu_y) = n$. By using the forms in Eqs. (54) and (55) taken by the arguments of the polynomial, we see that the phase of each term is
$$ (\mu_y + 2\nu_y) \frac{\phi}{2m_x} - (\mu_x + 2\nu_x) \frac{\phi}{2m_y} = \frac{i_y m_y - i_x m_x}{2m_x m_y} \phi, \quad (59) $$

with $i_\rho = \mu_\rho + 2\nu_\rho$ a non-negative integer such that $i_x m_x + i_y m_y = n$. From here we deduce that the phase accumulated by the polynomial following an increase of $2\pi$ in $\phi$ is given by

$$ \Phi_{\nu} = \frac{i_y m_y - i_x m_x}{m_x m_y} \pi. \quad (60) $$

## Appendix C: Mathematica code for the recurrence relations

We give here simple code for implementing in Wolfram's Mathematica the recurrence relations for the coherent states and their overlap. The recurrence formula for the coherent state given in Eq. (27) can be implemented as

```math
Frecur[\[Alpha]x_, \[Alpha]y_, \[Beta]x_, \[Beta]y_, mx_, my_, n_] :=
Module[{list}, list = {1};
For[i = 1, i < n + 1, i++,
  list = Append[
    list, (mx (\[Alpha]x If[i > mx, list[[i - mx + 1]]],
      If[i == mx, 1, 0]] +
      2 \[Beta]x If[i > 2 mx, list[[i - 2 mx + 1]]],
      If[i == 2 mx, 1, 0])) +
    my (\[Alpha]y If[i > my, list[[i - my + 1]]],
      If[i == my, 1, 0]] +
      2 \[Beta]y If[i > 2 my, list[[i - 2 my + 1]]],
      If[i == 2 my, 1, 0]])]/i]]]; Last(list)]
```

while that for the overlap of coherent states given in Eq. (37) can be implemented as

```math
Grecur[mx_, my_, n_, wx_, wy_] :=
Module[{list}, list = {1};
For[i = 1, i < n + 1, i++,
  list = Append[
    list, ((mx wx If[i > mx, list[[i - mx + 1]], If[i == mx, 1, 0]] +
      my wy If[i > my, list[[i - my + 1]],
      If[i == my, 1, 0]]) + (wx^2 If[i > 2 mx,
        list[[i - 2 mx + 1]], If[i == 2 mx, 1, 0]] +
      2 wx wy If[i > mx + my, list[[i - mx - my + 1]]],
      If[i == mx + my, 1, 0]] +
      wy^2 If[i > 2 my, list[[i - 2 my + 1]]],
      If[i == 2 my, 1, 0]]))/(i^2)]]; Last(list)]
```

## References

[1] Y. C. Lin, T.-H. Lu, K.-F. Huang, and Y.-F. Chen, "Model of commensurate harmonic oscillators with SU(2) coupling interactions: analogous observation in laser transverse modes," Phys. Rev. E, vol. 85, no. 4, p. 046217, 2012.

[2] Y. Shen, "Rays, waves, SU(2) symmetry and geometry: toolkits for structured light," J. Opt., vol. 23, no. 12, p. 124004, 2021.

[3] P. H. Tuan, K. T. Cheng, and Y. Z. Cheng, "Generating high-power Lissajous structured modes and trochoidal vortex beams by an off-axis end-pumped Nd: YVO 4 laser with astigmatic transformation," Opt. Express, vol. 29, no. 15, pp. 22957–22965, 2021.

[4] P. H. Tuan, W. C. Tsai, and K. T. Cheng, "Beam stability improvement of high-power Lissajous modes by an off-axis

pumped YVO 4/Nd: YVO 4 laser," Opt. Continuum, vol. 1, no. 8, pp. 1696–1703, 2022.

[5] X.-L. Zheng, Y.-H. Fang, W.-C. Chung, C.-L. Hsieh, and Y.-F. Chen, "Exploring the origin of Lissajous geometric modes from the ray tracing model," Photonics, vol. 11, no. 5, p. 456, 2024.

[6] P. B. Charles and K. B. Wolf, "The 2:1 anisotropic oscillator, separation of variables and symmetry group in Bargmann space," J. Math. Phys., vol. 16, no. 11, pp. 2215–2223, 1975.

[7] K. Tschernig, D. Guacaneme, O. Mhibik, I. Divliansky, and M. A. Bandres, "Observation of Boyer-Wolf Gaussian modes," Nat. Commun., vol. 15, no. 1, p. 5301, 2024.

[8] D. A. Kessler and I. Freund, "Lissajous singularities," Opt. Lett., vol. 28, no. 2, pp. 111–113, 2003.

[9] I. Freund, "Bichromatic optical Lissajous fields," Opt. Commun., vol. 226, nos. 1-6, pp. 351–376, 2003.
[10] I. Freund, "Polarization critical points in polychromatic optical fields," *Opt. Commun.*, vol. 227, nos. 1-3, pp. 61–71, 2003.

[11] I. Freund, "Polychromatic polarization singularities," *Opt. Lett.*, vol. 28, no. 22, pp. 2150–2152, 2003.

[12] S. Kerbstadt, D. Timmer, L. Englert, T. Bayer, and M. Wollenhaupt, "Ultrashort polarization-tailored bichromatic fields from a CEP-stable white light supercontinuum," *Opt. Express*, vol. 25, no. 11, pp. 12518–12530, 2017.

[13] Á. Jiménez-Galán, et al., "Control of attosecond light polarization in two-color bicircular fields," *Phys. Rev. A*, vol. 97, no. 2, p. 023409, 2018.

[14] H. Yan and B. Lü, "Dynamical evolution of Lissajous singularities in free-space propagation," *Phys. Lett. A*, vol. 374, no. 35, pp. 3695–3700, 2010.

[15] A. Fleischer, O. Kfir, T. Diskin, P. Sidorenko, and O. Cohen, "Spin angular momentum and tunable polarization in high-harmonic generation," *Nat. Photonics*, vol. 8, no. 7, pp. 543–549, 2014.

[16] E. Pisanty, et al., "Knotting fractional-order knots with the polarization state of light," *Nat. Photonics*, vol. 13, no. 8, pp. 569–574, 2019.

[17] J. D. Louck, M. Moshinsky, and K. B. Wolf, "Canonical transformations and accidental degeneracy. I. The anisotropic oscillator," *J. Math. Phys.*, vol. 14, no. 6, pp. 692–695, 1973.

[18] M. Moshinsky, J. Patera, and P. Winternitz, "Canonical transformation and accidental degeneracy. III. A unified approach to the problem," *J. Math. Phys.*, vol. 16, no. 1, pp. 82–92, 1975.

[19] M. V. Berry, "Quantal phase factors accompanying adiabatic changes," *Proc. R. Soc. Lond. A*, vol. 392, pp. 45–58, 1984.

[20] D. D. Holm, *Geometric mechanics-Part I: Dynamics and Symmetry*, Singapore, World Scientific Publishing Company, 2011.

[21] M. A. Alonso, *Phase-Space Optics Fundamentals and Applications*, chapter 8, New York, McGraw-Hill Companies, 2010.

[22] G. W. Forbes and M. A. Alonso, "Using rays better. I. Theory for smoothly varying media," *J. Opt. Soc. Am. A*, vol. 18, no. 5, pp. 1132–1145, 2001.

[23] M. A. Alonso and G. W. Forbes, "Stable aggregates of flexible elements give a stronger link between rays and waves," *Opt. Express*, vol. 10, no. 16, pp. 728–739, 2002.

[24] M. A. Alonso and M. R. Dennis, "Ray-optical Poincaré sphere for structured Gaussian beams," *Optica*, vol. 4, no. 4, pp. 476–486, 2017.

[25] M. R. Dennis and M. A. Alonso, "Swings and roundabouts: optical Poincaré spheres for polarization and Gaussian beams," *Philos. Trans. R. Soc. A: Math., Phys. Eng. Sci.*, vol. 375, no. 2087, p. 20150441, 2017.

[26] M. R. Dennis and M. A. Alonso, "Gaussian mode families from systems of rays," *J. Phy.: Photon.*, vol. 1, no. 2, p. 025003, 2019.

[27] R. Gutiérrez-Cuevas, M. R. Dennis, and M. A. Alonso, "Ray and caustic structure of Ince-Gauss beams," *New J. Phys.*, vol. 26, no. 1, p. 013011, 2024.

[28] G. W. Forbes, "Robust and fast computation for the polynomials of optics," *Opt. Express*, vol. 18, no. 13, pp. 13851–13862, 2010.

[29] Y.-F. Chen and K.-F. Huang, "Vortex structure of quantum eigenstates and classical periodic orbits in two-dimensional harmonic oscillators," *J. Phys. A: Math. Gen.*, vol. 36, no. 28, p. 7751, 2003.

[30] M. Sanjay Kumar and B. Dutta-Roy, "Commensurate anisotropic oscillator, SU(2) coherent states and the classical limit," *J. Phys. A: Math. Theor.*, vol. 41, no. 7, p. 075306, 2008.

[31] J. Moran and V. Hussin, "Coherent states for the isotropic and anisotropic 2D harmonic oscillators," *Quantum Rep.*, vol. 1, no. 2, pp. 260–270, 2019.

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

[33] A. Z. Goldberg, et al., "Quantum concepts in optical polarization," *Adv. Opt. Photon.*, vol. 13, no. 1, pp. 1–73, 2021.

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

[35] T.-H. Lu, Y. C. Lin, Y.-F. Chen, and K.-F. Huang, "Three-dimensional coherent optical waves localized on trochoidal parametric surfaces," *Phys. Rev. Lett.*, vol. 101, no. 23, p. 233901, 2008.

[36] J. J. Healy, M. Alper Kutay, H. M. Ozaktas, and J. T. Sheridan, *Linear Canonical Transforms: Theory and Applications*, vol. 198, New York, Springer, 2015.

[37] K. Y. Bliokh, M. A. Alonso, and M. R. Dennis, "Geometric phases in 2D and 3D polarized fields: geometrical, dynamical, and topological aspects," *Rep. Prog. Phys.*, vol. 82, no. 12, p. 122401, 2019.