Quantum Mechanical Approach to Bifurcation Point Detection in Hamiltonian Dynamical Systems

Energy level statistics of a bounded quantum system, whose classical dynamical system exhibits bifurcations, is investigated using the two-point correlation function (TPCL), which at the bifurcation points exhibits periodic spike oscillations owing to the accumulation of levels called the shell effect. The spike oscillations of the TPCL is analyzed by the reduced chi-squared value which deduced to exhibit abrupt increases at bifurcation points, thereby yielding a novel detection approach. Using this method, we attempt to numerically detect the bifurcation points of a lemon-shaped billiard.


Introduction
Bifurcation point detection is a crucial technique in the research field of nonlinear dynamical systems, which has played a significant role in understanding the nonlinear phenomena in the actual world [1]. There are several conventional approaches based on the Newton-Raphson method [2,3], and more recently, the meta-heuristic approaches, such as the particle swarm optimization (PSO), the differential evolution (DE) and the evolution strategy (ES), have garnered considerable attention [4][5][6][7]. Because these methods rely on an objective function obtained from the eigenvalues of the monodromy matrix (Jacobi matrix of variational equations) to evaluate the stability of each periodic point(fixed point), it is difficult to determine the bifurcation point when the matrix cannot be easily obtained. The derivation of the objective function is generally difficult in most nonlinear dynamical systems that cannot be solved analytically, and it is not easy to determine the position of the fixed points, let alone analyze its stability. The objective of this study is to propose a novel approach to the bifurcation-point detection for Hamiltonian dynamical systems that have quantum mechanical counterparts. Instead of the objective function, the bifurcation point is detected by the quantum mechanical data of the eigenenergy levels. The core principle that makes this possible is a physical phenomenon called the shell effect [9,10].
The shell effect is a phenomenon in which the degeneracy of periodic orbits due to the bifurcation produces a periodic strong accumulation of eigenenergy levels in the corresponding quantum system. This phenomenon has been reported to exist in bounded quantum systems such as atomic nuclei, metallic clusters, and mesoscopic semiconductor systems such as quantum dots, and quantum billiards, as well as in partially open quantum systems; in addition the related research field is still expanding [11][12][13][14][15][16][17][18][19][20][21]. On the other hand, its mechanism has been elucidated by the semiclassical theory [22][23][24], where the quantum effects of bifurcation in the statistical properties of energy levels have also been elucidated using an extended version of the Gutzwiller's trace formula that addressed the divergence problem at the bifurcation point [25]. Numerical attempts to verify the impact of the shell effect on energy level statistics have been made. Berry, Keating, and Prado investigated the level number variance (LNV) of the perturbed cat map at a saddle-node (tangent) bifurcation and reported that there is an additional contribution to the long-range spectral correlation called the "lift-off," which causes the LNV to increase rapidly at a certain correlation length determined by the semiclassical theory [26]. The lift-off was observed also for the pitchfork bifurcation of the coupled quartic oscillators by Gutiérrez et al. [27]. For the short-range spectral correlation, Makino, Harayama, and Aizawa numerically investigated the nearestneighbor level-spacing distribution (NNLSD) of the quantum oval billiard, and then reported anomalous accumulation between adjacent levels at the bifurcation [28]. These phenomena were recently investigated by Makino in terms of the two-point correlation function (TPCF), where the close relationship to the shell effect that causes periodic spike accumulation of levels is elucidated [29].
The TPCF is defined as the probability density of identifying two levels at a specific energy distance. In the research field of quantum chaology, this function is typically derived from the eigenenergy levels in a sufficiently small energy range [E, E + ∆E]. This is because the phase space structure of Hamiltonian system is generally energy-dependent, and fixing the dynamics for all eigenstates considered in the statistics provides a clear correspondence between the quantum and classical aspects. For the bifurcation detection method proposed in this paper to work well, this point must be taken into account for the energy range to be analyzed. We will discussed this point again in section Section 6. Note that for the systems with f degrees-of-freedom, the number of levels in the interval ∆E diverges as O(∆E/ f ) in the semiclassical limit → 0, indicating that the small interval ∆E still contains a sufficient number of levels to ensure good statistical significance in the deep semiclassical regime.
The billiard systems have a convenient scaling law that makes the phase space structure invariant with respect to energy E, and hence, the entire energy spectrum can be used in the analysis. In this paper we will analyze the one-parameter family of lemon-shaped billiard whose bifurcation points are well understood [30][31][32][33], and will attempt to estimate the bifurcation points by the quantum mechanical data of eigenenergy levels characterized by the TPCF.
The remainder of this paper is organized as follows. The lemon-shaped billiard is introduced in Section 2, while bifurcation parameters are determined analytically in Section 3. In Section 4, the TPCF of the quantum lemon billiard that exhibits periodic spike oscillations is analyzed. These oscillations are evaluated by the Piason's χ 2 -test in Section 5, where the estimated values of the bifurcation parameters are determined by χ 2 values obtained from the quantum mechanical data of eigenenergy levels. The quantum mechanical estimates are compared with the bifurcation parameters of the classical dynamical system, followed by a summary and discussion in Section 6.
2 Lemon-shaped billiard Figure 1(a) presents a schematic diagram of the lemon-shaped billiard whose boundary wall ∂D comprises two symmetrical arcs, where its curvatures are determined by a single parameter δ ∈ (0, 1]. For δ = 1(circular wall), the dynamical system is integrable, and the motion of a particle on the domain D is regular for any choice of initial conditions. For 0 < δ < 1, the dynamical system is non-integrable and the motion of the particle on D is regular or chaotic, depending on the initial condition.  Fig.1(a)], and α denotes the angle between the inner normal and the orbit reflected from the wall [34]. When δ = 1, the entire surface of the section is filled with the invariant tori[ Fig.1(b)]. As δ is altered from 1 to 0, the system transitions from the integrable to the non-integrable with successive bifurcations, while the surface of the section is filled with tori and chaos, as shown in Fig.1(c). The bifurcation parameters are derived in the next section.

BIFURCATIONS
The periodic orbits observed in this system are a countability infinite set that can be divided roughly into two families: the bouncing mode family, as illustrated in the Fig.2(a), comprising back-and-forth trajectories between upper and lower walls, and the glancing mode family presented in Fig. 2(b), which comprises trajectories with repeated shallow reflections along the wall.
The periodic orbits n = 2, 3, 4, · · · belonging to the bouncing mode family form the pairs comprising a stable periodic orbit and an unstable periodic orbit of the same period 2n, where R(δ) = √ 1 + δ 2 /δ represents the radius of circular arc ∂D, and w(δ) In addition, the length of the bouncing periodic orbit n at δ B n is

Figures 4 and 5 represent the Poincaré surfaces of the section and bifurcation diagram around
the bifurcation points δ B 2 and δ B 3 , respectively. The red and blue marks or lines in each figure represent the unstable (hyperbolic) and stable (elliptic) periodic points, respectively, which merge at the bifurcation points. The bifurcation structures analyzed in this study are similar to the saddle-node type of one-dimensional systems.
Each periodic orbit with a period 2(n + 1) belonging to the glancing mode family n = 1, 2, 3, · · · , emerges at δ G n ≡ γ n − γ 2 n − 1, γ n = 1 + 1/ tan 2 α(n) from the positions φ = 0 and 0.5, with an angle α(n) = nπ/2(n + 1), and exists in the parameter region δ ∈ [δ G n , 1), where δ G n for n = 1 − 3 are described as The length of the glancing periodic orbit n at δ G n is l G n = 4nR(δ G n ) cos α(n).  where the the red and blue lines correspond to the unstable and stable points with period 6, respectively. Table 1 presents the parameters δ B n and δ G n , including orbit lengths l B n and l G n for periodic orbits n = 1-9. Because δ G n has a property δ G n < δ G n+1 for all n, there are no periodic orbits belonging to the glancing-mode family in the region δ ∈ [0, δ G 1 ). Hence, we will focus on the region δ ≤ 0.5 where the bifurcations of bouncing periodic orbits are mainly observed, and also explore the bifurcation points δ B n from the quantum mechanical data of eigenenergy levels. Table 1 Bifurcation parameters of the bouncing and glancing periodic orbits, including the orbit length at each bifurcation point. and ψ(x, ±y) = ±ψ(x, y), the eigenenergy levels are divided into mutual independent components belonging to these four classes; hence, the unfolding transformation needs to be carried out separately as =N (E )/4 for each of the four components, thereby yielding four unfolded sets of levels. The TPCF analyzed in this study represents the probability density of identifying two levels at spacing L, and is defined for each of the energy level components using the level density d(x) = δ(x − ) as R 2 (L) = d(x − L/2)d(x + L/2) , where the bracket · · · stands for an averaging over x [37]. This quantity is suitable for studying the quantum-mechanical effects of bifurcation, as its relationship with classical periodic orbits is well understood in the semiclassical theory. Based on Gutzwiller's trace formula [25], TPCF is expressed by the periodic-orbit sum as where j labels each primitive periodic orbit and its repeatations, C j ( ) represents the amplitude factor determined by the degeneracy and stability of the orbit, S j ( ) denotes the action integral along the orbit j, which is defined here to include the Maslov index, and T j = ∂S j ( )/∂ represents the time period of the periodic orbit. The second term in the RHS of Eq. (3) is expected to disappear via the smoothing procedure over , TPCF is approximated by a sum of periodic functions, whose periods with respect to L are described as p j = 2π /T j . For the billiard problem, this period is rewritten by the orbit length l j as It should be noted that the creation of a periodic orbit j * across the bifurcation generates an additional contribution |C j * ( )| 2 cos [T j * L/ ] to Eq.(4) and the creation of this term can exert a significant impact on the property of the TPCF, if the amplitude C j * is relatively nonnegligible in an infinite series j . Such a possibility is expected to emerge at the bifurcation point where the extended trace formula, which is derived from the improved stationaryphase-approximation, predicts a significantly large and non-divergent value of C j * (also refer to Ref. [22][23][24]). Furthermore, this should trigger periodic oscillations of the p j * period in the behavior of the TPCF.   this study is ultimately determined by the superposition of four TPCFs obtained respectively from the unfolded energy levels of the four parity-symmetry classes.
It is widely known that in time-reversal invariant quantum systems with classically fully chaotic counterpart, a universality proposed by Bohigas, Giannoni, and Schmit(BGS) exists [38], such that the unfolded energy levels in the semiclassical limit exhibit the same fluctuation properties as predicted by the gaussian orthogonal ensemble (GOE) statistics of the random matrix theory [37], which provides the TPCF in the following form where σ(L) = sin(πL)/(πL)[also refer to the blue curves in Figure.6]. While in quantum systems with a classically integrable counterpart, another universality proposed by Berry and Tabor exists [39], such that the unfolded energy levels in the semiclassical limit exhibit the same fluctuation properties as the random number from the Poisson point process, which gives R Poisson 2 (L) = 1[refer to the yellow lines of Figure 6]. The theoretical underpinnings of these two universalities remain a subject in the research field of quantum chaology [40][41][42][43][44]. For a quantum system whose classical dynamical system comprises regular and chaotic motions, the TPCF fits neither R Poisson 2 (L) nor R GOE 2 (L) as shown in Figs.6(b) and 6(c). In this case, it is useful to introduce their interpolation formula[refer to Fig.6(d)] whose physical meaning is supported by the Berry-Robnik level statistics [45,46]. In the Berry-Robnik level statistics, Eq. (7) provides the TPCF of eigenenergy levels, which is a product of the statistically independent superposition of two spectral components following the GOE and Poisson statistics, while the relative weight ρ ∈ [0, 1] is assumed to coincide with the relative phase volume in the Liouville measure of the chaotic component. In this research, we do not go into its physical meaning and deal with ρ as a fitting parameter. √ /l B n , which is determined by the orbit length (2) of the bifurcating orbit. It is quite interesting that the TPCF in each figure exhibits remarkable spike oscillations whose period is well approximated by the fundamental period p n of the series j * |C j * ( )| 2 cos [T j * L/ ] in Eq.(4); hence, the oscillation is indeed contributed by the bifurcating periodic orbits j * ∈ pair n. In the next section, we propose an effective method to determine the bifurcation points from the quantum mechanical data of eigenenergy levels characterized by the TPCF.

Bifurcation point detection
As shown in the previous section, the bifurcation of the periodic orbit j * ∈ pair n triggers spike oscillations in the TPCF, which are R 2 (L) > 1 at L = p n , = 0, 1, 2, · · · . Consequently, R 2 (L) can no longer be approximated by conventional functions derived from the GOE Pearson's χ 2 -test evaluates a measure χ 2 , which is a sum of differences between observed and expected outcome frequencies. For a given TPCF R 2 (L), the observed frequency for each class i = 1, 2, 3, · · · , K is calculated by using the class interval ∆L and total number N of eigenenery levels as o i = N R 2 (L i )∆L, while the expected frequency is calculated as e i = NR 2 (L i , ρ)∆L, where the interpolation parameter ρ is determined by the best-fitting curveR 2 (L, ρ) to the numerical data R 2 (L). Then, χ 2 and its degrees of freedom ν are obtained as /e i and ν = K − 1, where ∆L and K are determined to hold e i >> 1 for every i, and determined to hold e i > 5 in our analysis. If the reduced chisquared value defined as χ 2 ν ≡ χ 2 /ν is less than 2, the conventional criteria for statistical significance, the goodness of the fitting is ascertained to be sufficient, and the null hypothesis that the numerical data R 2 (L) fits the formulaR 2 (L) is not rejected. Conversely, in the case χ 2 ν ≥ 2, it is determined that R 2 (L) deviates significantly fromR 2 (L). Such a case occurs when R 2 (L) exhibits a spike oscillation; hence, maximal points significantly larger than 2 in the numerical plots of χ 2 ν (δ) provide candidates for the bifurcation points.  n , n = 2 − 9 determined in the classical dynamical system. Surprisingly, χ 2 ν (δ) in both energy ranges exhibits abrupt increases near the bifurcation points, each of which has its maximul value above the criterion of 2, as indicated by the red circle. It should also be noted that χ 2 ν (δ) has its maximals closer to the bifurcation points in the higher energy region than in the lower energy region.
Tables 2(a)-(b) present the maxima of χ 2 ν (δ) for each of the two energy ranges, in order from the largest value, including their position δ qm n in the parameter space. The analytical solution of the bifurcation parameter δ B n , which is closest to the position of each maximum, is also listed in the Table. In both energy ranges, the quantum mechanical data δ qm n obtained from large maximum value agree well with δ B n , n = 2, 3, 4, · · · , thereby providing an optimal estimate of the bifurcation points; in addition, the agreement is better in the higher energy region, which is in the deeper semiclassical regime. However, for δ qm n obtained from the small maxima, there are no counterparts of δ B n that exhibit erroneous estimates. This result imply that the larger the χ 2 ν (δ qm n ) value, the better the estimate δ qm n is for the bifurcation parameter δ B n . This study was dedicated to developing a novel method for obtaining estimates of the bifurcation points in Hamiltonian dynamical systems using the eigenenergy levels of its quantum systems characterized by TPCF. Owing to the strong accumulation of levels caused by the shell effect, the TPCF at the bifurcation point exhibited periodic spike oscillations, whose period was well approximated by the semi-classical theory. We introduced the Pearson's χ 2 test to capture the spike oscillations of TPCF and determined that the reduced chi-squared value χ 2 ν (δ) exhibits abrupt increases at the bifurcation points, thereby providing an optimal estimate of the bifurcation parameters. The accuracy of this estimation was determined to be better in the higher energy region, where the quantum system approaches the semiclassical domain.
Note that the method proposed in this study is solely effective for Hamiltonian dynamical systems that have a quantum mechanical counterpart. Because this method does not require the objective function or its derivative, it is very useful for analyzing systems whose bifurcation parameters are difficult to determine via the conventional analysis.
The proposed method described above has been applied to the lemon-shaped billiards with good results. And it could be effectively applied to more general Hamiltonian dynamical systems if the following condition is satisfied.
The phase space geometry of Hamiltonian dynamical system is generally energydependent, and to obtain a clear correspondence between the quantum and classical aspects, the eigenenergy levels must be obtained from a sufficiently small interval ∆ . However, when the interval is extremely small, our method may not detect bifurcations of short periodic orbits that produce level accumulations of long period. In our method, the interval ∆ should contain at least one accumulation point, which provides the condition for the lower limit ∆ > p( ). On the other hand, there is also the upper limit from another context. should be conducted in the future to collect specific examples and to determine the scope of effective application. In the following, the above argument will be confirmed for the TPCF of the lemon-shaped billiard we analyzed. The gray region in Figure 10 provides the energy range [ , + ∆ ] whose interval ∆ satisfies the inequality (8). We also plotted the three energy ranges (a)[0,4000], (b)[10000,12000], and (c)[160000,168000], where the ranges (a) and (c) were adopted to identify the bifurcation points in Section 4. Figures 11(a)-(c) present the TPCF at δ B 2 for the three energy ranges (a)-(c) of Fig. 10, respectively. The ranges (b) and (c) fully satisfy the condition (8), and the spikes of R 2 (L) are hence maintained over a region longer than the period p 2 ( ), as shown in Figs.11(b) and (c). While the range (a) has a short coherence length, and the spikes are collapsed within a range shorter than the period p 2 as shown in Fig.11(a).
The reduced chi-square statistics is also applicable to the level-spacing distribution that is more popular in the research field of level statistics. In the present case, where some of the spectral components contributed from the bifurcating orbits, are strongly accumulated, it can be expected that the level-spacing distribution is consistent with a sub-Poisson or  asymptotic Poisson distribution [47,48], and the reduced chi-square value is expected to exhibit abrupt increases in this statistics as well. This possibility will be studied elsewhere.