Harmonic oscillators are universally recognized as one of the most important paradigms for exploring quantum-classical correspondence. Under the paraxial approximation, the transverse part of the wave equation for spherical cavities can be mathematically analogous to the Schrödinger equation for two-dimensional (2D) harmonic oscillators [1]. Accordingly, various high-order transverse modes can be generated with the specially designed spherical laser cavities to analogously manifest the quantum wave function. The eigenfunctions of 2D quantum harmonic oscillators can be solved as the Hermite–Gaussian (HG) functions in rectangular coordinates or the Laguerre–Gaussian (LG) functions in polar coordinates [2]. The selectively diode-end-pumped solid-state lasers have been widely employed to generate both HG and LG functions from ground order to very high order ^{[3][4][5][6][7]}[3,4,5,6,7]. Additionally, the same laser technology was exploited to generate the so-called geometric modes in the degenerate cavities, which clearly revealed the ray-wave duality in the spatial domain ^[8][9][8,9]. In mesoscopic quantum phenomena, the degeneracy of energy levels was found to play an important role in the connection between the conductance fluctuation and the classical periodic orbits (POs) [10]. Similarly, the emergence of geometric modes was verified to originate from the degeneracy of eigenfrequencies in laser resonators ^[11][12][11,12]. The Lissajous stationary modes are one of the most remarkable geometric modes generated from the astigmatic laser cavities. Theoretically, the Schrödinger coherent state for the one-dimensional (1D) harmonic oscillator can be straightforwardly extended to the 2D harmonic oscillator to obtain the stationary coherent states that exactly spatially correspond to the Lissajous figures.

Additionally, 2D quantum billiards are another pedagogical model for comprehending the connection between quantum and classical worlds. Various dynamical features can be straightforwardly studied from the model of quantum billiards by changing the geometry. One main branch of research on billiard systems is focused on nonintegrable and irregular shapes to understand the characteristics in the field of quantum chaos ^{[13][14][15][16][17][18][19][20]}[13,14,15,16,17,18,19,20]. Classically, the chaotic nature renders all the orbits in a chaotic system as being unstable. An interesting finding in nonintegrable quantum billiards is the emergence of eigenstates on unstable POs, called scarred modes ^[21][22][23][21,22,23]. Quantum scars have been searched and analyzed in mesoscopic systems ^{[24][25][26][27][28][29][30][31][32][33][34]}[24,25,26,27,28,29,30,31,32,33,34]. Due to the similarity of the equations for different types of waves, scars have been observed in microwaves ^{[35][36][37][38][39][40][41]}[35,36,37,38,39,40,41]. Quantum scars have not only been confirmed from the accumulation of spin–orbit-coupled atomic gases for specific energies [42] but also generated in the 2D harmonic oscillators with local impurities ^[43][44][45][43,44,45]. Furthermore, quantum many-body scars have been hypothesized to cause weak ergodicity breaking and the unexpectedly slow thermalization of cold atoms ^{[46][47][48][49][50][51][52][53][54][55]}[46,47,48,49,50,51,52,53,54,55]. The similar phenomenon dynamical scar has also been experimentally found in a driven fraction system [56]. Nevertheless, the overall number of scarred modes is quite few. The eigenstates in nonintegrable billiards are mostly widely distributed in the coordinate space [13], often exhibiting common features of quasi-linear ridge structures [57].

Compared with non-integrable billiards, stable POs are generally abundant in integrable billiards with symmetrical shapes ^{[11][16][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75]}[11,16,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75]. The quantum wavefunctions associated with stable POs have been found to play a key role in ballistic transport ^{[76][77][78][79]}[76,77,78,79], quantum pointer states and decoherence ^{[80][81][82][83][84][85][86][87][88][89][90][91][92][93][94][95]}[80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95], universal conductance fluctuations ^{[95][96][97][98][99]}[95,96,97,98,99], and chaos-assisted quantum tunneling ^{[100][101][102][103][104]}[100,101,102,103,104]. Ballistic transport means that the mean free path of the particle is significantly longer than the size of the medium through which the particle travels. In addition to microwave cavities, quantum billiards can be analogously explored with the wave systems including optical fibers ^[105][106][105,106], optical resonators ^{[107][108][109][110][111][112][113][114][115][116]}[107,108,109,110,111,112,113,114,115,116], vibrating plates and acoustic waves ^{[117][118][119][120][121][122][123]}[117,118,119,120,121,122,123], and liquid surface waves ^{[123][124][125][126][127][128]}[123,124,125,126,127,128]. Noticeably, it has been confirmed that the vertical-cavity surface-emitting lasers (VCSELs) with a unique longitudinal wave vector k_z and the lateral oxide confinements can be modeled as 2D wave billiards with hard walls. To be brief, theoretical research on quantum billiards was intensively performed in the last century, and later, some researchers’ interests shifted to applied fields such as laser resonators.

2. Quantum Billiards

The quantum analogy of a classical billiard is called a quantum billiard. For the classical billiard with the 2D region denoted by R, the corresponding potential in quantum mechanics is given as The time-independent Schrödinger equation for the potential defined in Equation (1) can be expressed as the Helmholtz equation: where 𝑘=(2𝜇𝐸/ℏ)1/2; E and 𝜇 are energy and mass of the particle, and ℏ is the reduced Planck’s constant. The homogeneous Dirichlet boundary condition is given due to the condition of 𝑉=∞. To be brief, quantum billiard is defined using the Helmholtz equation in R with the Dirichlet boundary condition.

3. Square Billiard

The spatial distributions of quantum wave functions corresponding to classical POs ^{[76][77][78][79]}[76,77,78,79] have been an intriguing phenomenon in open ballistic cavities. Semiclassical PO theory has been used to explain the scarred wave functions in chaotic billiards ^[21][22][23][21,22,23]. Nevertheless, it is pedagogically useful for comprehending the quantum-classical correspondence in mesoscopic physics to fully develop the connection between quantum eigenfunctions and classical POs in integrable systems. One of the simplest integrable billiards is the square billiard [11]. In a square billiard, each family of POs can be specified with three parameters (p, q, and ϕ), where p and q are two positive integers describing the number of reflections with the horizontal and vertical boundaries, and ϕ (−π < ϕ < π) is associated with the wall positions of specular reflection points ^[11][65][11,65]. Alternatively, the parameter ϕ may also be linked to the starting point of the classical particle. Figure 1 depicts some examples for POs in a square billiard. The trajectory can be seen to constitute a single, non-repeated orbit when p and q are co-prime. When p and q have a common factor m, the trajectory is an orbit family that corresponds to m primitive POs with indices of (p/m, q/m, and ϕ/m). For a square billiard with the region in 0≤𝑥,𝑦≤𝑎, the eigenfunctions are given as where the quantum numbers m and n are positive integers. The eigenvalues corresponding to the eigenfunctions 𝜓𝑚,𝑛(𝑥,𝑦) are given using 𝐸(𝑚,𝑛)=ℏ2𝑘2𝑚,𝑛/(2𝜇), where the wave numbers 𝑘𝑚,𝑛 are expressed as

shows the wave patterns for several sets of quantum numbers . From Bohr’s correspondence principle, the classical limit of a quantum system should be asymptotically obtained when the quantum numbers are sufficiently large. However, the conventional eigenfunctions of a square billiard cannot reveal the features of classical POs no matter how large the quantum numbers are. The quantum states relevant to the classical POs have been verified to be the superpositions of the nearly degenerate eigenstates. For a given central order , the nearly degenerate condition can be derived from the differential of the eigenvalue function given as

Setting 𝑑𝐸(𝑚,𝑛)|𝑚𝑜,𝑛𝑜=0 leads the tangent of the constant-energy contour as Since both quantum numbers (𝑚,𝑛) are positive integers, the slope −𝑑𝑛/𝑑𝑚 must be a rational number. From Equation (6), the nearly degenerate condition can be generalized as 𝑚𝑜/𝑛𝑜=𝑞/𝑝, where p and q are coprime positive integers. Figure 3 depicts the spectrum 𝑘𝑚,𝑛 as a function of the ratio 𝑚/𝑛 for a square billiard with 1≤𝑚,𝑛≤700. The spectrum conspicuously reveals that the eigenvalues are clustered in the vicinity of 𝑚/𝑛=𝑞/𝑝 to form valley structures. Clustering means that levels with very different quantum numbers have very similar energies. From the condition 𝑚𝑜/𝑛𝑜=𝑞/𝑝, the central eigenstate for the coherent superposition can be in terms of a single parameter N as 𝑚𝑜=𝑞𝑁 or 𝑛𝑜=𝑝𝑁. The slope −𝑑𝑛/𝑑𝑚=𝑞/𝑝 signifies that the quantum numbers for the nearly degenerate eigenstates around the central mode can be given using 𝑚=𝑞𝑁+𝑝𝐾 and 𝑛=𝑝𝑁−𝑞𝐾 with the integer index K in a small range. Consequently, the coherent superposition of the nearly degenerate eigenstates around the central mode can be generalized as

where 𝜙 is the phase factor in the range of −𝜋≤𝜙≤𝜋 and (2𝑀+1) means the total number of the superposed eigenstates. Note that the parameter ϕ corresponds to the starting point of the classical particle shown in Figure 1. Under the circumstance of 𝑁>>𝑀, the eigen-energies of the superposed eigenstates can be confirmed to be close to a constant energy. Figure 4 illustrates the wave patterns |Ψ(𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙)|2 calculated using Equation (7) with N = 100, M = 5, and 𝜙=𝜋/2 for eigenstates clustered around the indices (p,q) shown in Figure 3. The wave patterns of |Ψ(𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙)|2 are evidently localized on the classical POs. The velocity direction of the trajectory can be straightforwardly determined with the relation of 𝑘𝑥/𝑘𝑦=𝑚/𝑛=𝑞/𝑝. By way of explanation, the wave function in Equation (7) is not a strictly stationary state since the eigenstate components are not exactly degenerate for the Hamiltonian H. Nevertheless, Δ𝐻/<𝐻> will rapidly approach to zero as 𝑁→∞ for a small M. Therefore, the coherent state in Equation (7) can be regarded as a quasi-stationary state in the mesoscopic region.

The trajectorial equations for POs can be derived from the quantum coherent state in Equation (7) for the central maximum of the wave intensity. Using the identity sin𝜃=(𝑒𝑖𝜃−𝑒−𝑖𝜃)/(2𝑖), the representation of the coherent state in Equation (7) can be organized as

where Θ±𝑔(𝑥,𝑦)=𝜋(𝑞𝑥±𝑝𝑦)/𝑎 and Θ±𝑡(𝑥,𝑦)=𝜋(𝑝𝑥±𝑞𝑦)/𝑎. The function 𝐷𝑀(𝜃) in Equation (9) is the Dirichlet kernel that exhibits the periodic maxima of the intensity at 𝜃=2𝑛𝜋 for any integer n. Using the periodic maximal characteristic of the Dirichlet kernel, the parametric equations for the central maxima of the intensity |Ψ(𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙)|2 can be generalized as 𝐴 𝑥+𝐵 𝑦±𝜙=2𝑛𝜋, where −𝐴 /𝐵 =𝜂 represents the slope. From Equation (8) and Θ±𝑡(𝑥,𝑦)=𝜋(𝑝𝑥±𝑞𝑦)/𝑎, the slopes for all parametric equations can be found to be two cases of 𝜂=±𝑝/𝑞. Furthermore, all parametric equations can be confirmed to exactly correspond to the trajectorial lines of classical POs in a square billiard. The initial position (𝑥𝑜,𝑦𝑜) and the velocity (𝑣𝑥,𝑣𝑦) in classical dynamics can be linked to Equation (8) using the condition of 𝐴 𝑥𝑜+𝐵 𝑦𝑜±𝜙=2𝑛𝜋 as well as 𝑣𝑦/𝑣𝑥=(𝑑𝑦/𝑑𝑡)/(𝑑𝑥/𝑑𝑡)|𝑥𝑜,𝑦𝑜=𝜂. From the result of 𝜂=±𝑝/𝑞, the velocity (𝑣𝑥,𝑣𝑦) can be verified to be consistent with the classical dynamics 𝑣𝑦/𝑣𝑥=𝑝/𝑞. To sum up, the trajectorial equations for the classical POs can be analytically extracted from the quantum coherent states in Equation (7). It is worthwhile to mention that the lines of the phase functions Θ±𝑔(𝑥,𝑦) and the lines of the trajectorial functions Θ±𝑡(𝑥,𝑦) in Equation (8) are mutually orthogonal.

4. Equilateral Triangular Billiard

Square billiard is a classically separable and integrable system, whereas the equilateral triangle billiard is an integrable but non-separable system. In theory, the correlation between the quantum level clustering and the classical POs was deeply discussed from the representation of the quantum coherent states. In experiments, the oxide-confined VCSEL devices were fabricated in the shape of an equilateral triangle to analogously manifest the quantum level distribution and the spatial features of the wave functions. For an equilateral-triangular billiard with three vertices at (0,0), (𝑎/2,3−−√𝑎/2), and (−𝑎/2,3−−√𝑎/2), the eigenfunctions are given using ^[62][129][62,130]

where the quantum numbers m and n are nonnegative integers, and the superscripts (o) and (e) denote the two types of degenerate modes with odd and even symmetries, respectively. The eigenvalues corresponding to the eigenfunctions 𝜓(𝑒)𝑚,𝑛(𝑥,𝑦) and 𝜓(𝑜)𝑚,𝑛(𝑥,𝑦) are given using 𝐸(𝑚,𝑛)=ℏ2𝑘2𝑚,𝑛/(2𝜇), where the wave numbers 𝑘𝑚,𝑛 are expressed as

Figure 5 shows the wave patterns of |𝜓(𝑒)𝑚,𝑛(𝑥,𝑦)|2 for several sets of quantum numbers (𝑚,𝑛). Since the wave patterns for |𝜓(𝑒)𝑚,𝑛(𝑥,𝑦)|2 and |𝜓(𝑜)𝑚,𝑛(𝑥,𝑦)|2 are the same in the spatial patterns, only the case of |𝜓(𝑒)𝑚,𝑛(𝑥,𝑦)|2 is presented in Figure 5. Due to the setting of the equilateral triangle, all wave patterns can be found to be symmetric with respect to the y axis. Like the results for a square billiard, the conventional eigenstates for an equilateral-triangular billiard cannot exhibit the spatial properties of classical POs, even in the correspondence limit of large quantum numbers. The eigenfunctions 𝜓(𝑒)𝑚,𝑛(𝑥,𝑦) and 𝜓(𝑜)𝑚,𝑛(𝑥,𝑦) in Equations (10) and (11) are the standing-wave representation. The traveling-wave representation is more convenient for constructing the coherent states relevant to classical POs. In terms of 𝜓(𝑒)𝑚,𝑛(𝑥,𝑦) and 𝜓(𝑜)𝑚,𝑛(𝑥,𝑦), the traveling-wave representation is given using 𝜓(±)𝑚,𝑛(𝑥,𝑦)=𝜓(𝑒)𝑚,𝑛(𝑥,𝑦)±𝑖𝜓(𝑜)𝑚,𝑛(𝑥,𝑦) , where the symbols (+) and (−) denote the forward and backward characteristics, respectively. Consequently, the wave functions for 𝜓(±)𝑚,𝑛(𝑥,𝑦) can be given using Note that 𝜓(+)𝑚,𝑛 and 𝜓(−)𝑚,𝑛 form a conjugate pair with identical spatial patterns. Again, the quantum coherent states related to the classical POs can be formed via superposition of the nearly degenerate eigenstates. As discussed in the case of square billiard, the nearly degenerate condition for equilateral triangle billiards with the central order (𝑚𝑜,𝑛𝑜) can be given using

Since the slope −𝑑𝑛/𝑑𝑚 needs to be a rational number, the nearly degenerate condition in Equation (14) can be generalized as (2𝑚𝑜−𝑛𝑜)/(2𝑛𝑜−𝑚𝑜)=𝑞/𝑝 with p and q being coprime positive integers. Figure 6 depicts the spectrum 𝑘𝑚,𝑛 as a function of the ratio (2𝑚−𝑛)/(2𝑛−𝑚) for an equilateral triangle billiard with 1≤𝑚,𝑛≤700. The eigenvalues can be seen to be clustered in the vicinity of (2𝑚−𝑛)/(2𝑛−𝑚)=𝑞/𝑝 to display the valley structures. Obviously, the level clustering is certainly accompanied by the emergence of the gap. From the condition (2𝑚𝑜−𝑛𝑜)/(2𝑛𝑜−𝑚𝑜)=𝑞/𝑝, the central eigenstate for the coherent superposition can be given by 𝑚𝑜=(2𝑞+𝑝)𝑁 and 𝑛𝑜=(2𝑝+𝑞)𝑁 with a single parameter N. Combining with −𝑑𝑛/𝑑𝑚=𝑞/𝑝 from Equation (14), the coherent superposition of nearly degenerate eigenstates around the central mode can be expressed as

where 𝜙 is the phase factor in the range of −𝜋≤𝜙≤𝜋. Under the circumstance of 𝑁>>𝑀, the eigen-energies of the superposed eigenstates can be found to be nearly a constant energy of

Figure 7 illustrates the wave patterns of |Ψ(±,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙)|2 calculated using Equation (15) with N = 100, M = 5, and 𝜙=𝜋/2 for eigenstates with the indices (𝑝,𝑞) shown in Figure 6. The wave patterns of |Ψ(+,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙)|2 can be seen to be precisely concentrated on the classical POs. Since Ψ(+,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙) and Ψ(−,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙) form a conjugate pair, the spatial patterns are completely identical. calculated using Equation (15) with = 100, = 5, and for eigenstates clustered around the indices (, ).

The same as in the case of a square billiard, the trajectory equations for POs can be derived from Equation (15) from the central maximum of the wave intensity. Using sin𝜃=(𝑒𝑖𝜃−𝑒−𝑖𝜃)/(2𝑖), the representation in Equation (15) can be organized as

where The representation for Ψ(−,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙) can be given using the conjugate of Ψ(+,𝑝,𝑞)𝑁,𝑀(𝑥,𝑦;𝜙). Using the maximal feature of the Dirichlet kernel, the parametric equations for the central maxima of |Ψ(+,𝑝,𝑞)𝑁,𝑀| can be deduced as Θ±𝑡 𝑗(𝑥,𝑦,𝜙)=2𝑛𝜋 with 𝑗=1,2,3. Consequently, the classical POs of the equilateral-triangular billiard can be confirmed to be constituted by six independent line equations with different slopes. Using the form 𝐴𝑗𝑥+𝐵𝑗𝑦+𝜙=2𝑛𝜋 to express the trajectory equations, the slopes can be generalized as −𝐴𝑗/𝐵𝑗=±𝜂𝑗 with 𝜂1=(2𝑝+𝑞)/3−−√𝑞, 𝜂2=(𝑝+2𝑞)/3−−√𝑝, and 𝜂3=(𝑝−𝑞)/3−−√(𝑝+𝑞). Similar to the quantum coherent state for a square billiard, the lines of the phase functions Θ±𝑔 𝑗(𝑥,𝑦) and the lines of the trajectorial functions Θ±𝑡 𝑗(𝑥,𝑦) in Equation (17) are mutually orthogonal for 𝑗=1,2,3 , respectively.

5. Circular Billiard

Circular billiard is another pedagogical paradigm in classically separable and integrable systems. The azimuthal and radial components of the eigenfunctions of a circular billiard are the form of exp(imϕ) and the Bessel function of the first kind, respectively. Helically phased light beams with the azimuthal phase form of exp(imϕ) are well known to carry an orbital angular momentum (OAM) of mη per photon, where m is an integer ^[130][131][131,132]. In ray dynamics, the function of a circular billiard is the same as the transverse confinement of a cylindrical waveguide for light. Consequently, the propagation-invariant solutions of the Helmholtz equation in a cylindrical waveguide can be in terms of the Bessel beams with well-defined OAM ^[132][133]. The OAM or optical vortex (OV) of light has been widely used in numerous applications, such as generating OAM-entangled photon pairs ^[133][134][134,135], trapping and rotating micron and submicron objects ^{[135][136][137]}[136,137,138], generating astrophysical OAM light ^[138][139], assembling DNA biomolecules ^[139][140], OAM-based microscopy and imaging ^[140][141], super-diffraction limit imaging ^[141][142], and optical communication ^[142][143]. The eigenstates in polar coordinates for a circular billiard with radius R are given using where 𝑚∈𝛧, 𝑛∈𝛮, and 𝐽𝑚(•) are the Bessel functions of the first kind with order m. The quantum numbers m and n are the quantization of the azimuthal and radial oscillations, respectively. The eigenvalues for 𝜓𝑚,𝑛(𝑟,𝜃) are given using 𝑘𝑚,𝑛=𝑥𝑚,𝑛/𝑅 with 𝑅=𝑎/2, where 𝑥𝑚,𝑛 is the nth zero of 𝐽𝑚(𝑥) and a is the billiard diameter. Figure 8 shows the wave patterns for the function |Re[𝜓𝑚,𝑛(𝑟,𝜃)]|2 with different quantum numbers (m, n). Here, the real part of the eigenfunction 𝜓𝑚,𝑛(𝑟,𝜃) is purposely used for revealing the nodal structures in the radial and azimuthal directions associated with the indices n and m, respectively. Unlike square and equilateral triangular billiards, the nearly degenerate condition for a circular billiard cannot be straightforwardly derived from the eigenvalues 𝑘𝑚,𝑛 determined with 𝐽𝑚(𝑘𝑚,𝑛𝑅). The Wentzel–Kramers–Brillouin (WKB) method was nicely used to analytically obtain the nearly degenerate condition for a circular billiard. The eigenvalues 𝑘𝑚,𝑛 derived from the WKB method is given using ^[143][144]