Numerical verification of Percival's conjecture in a quantum billiard

In order to verify Percival's conjecture [J. Phys. B 6,L229 (1973)] we study a planar billiard in its classical and quantum versions. We provide an evaluation of the nearest-neighbor level-spacing distribution for the Cassini oval billiard, taking into account relations with classical results. The statistical behavior of integrable and ergodic systems has been extensively confirmed numerically, but that is not the case for the transition between these two extremes. Our system's classical dynamics undergoes a transition from integrability to chaos by varying a shape parameter. This feature allows us to investigate the spectral fluctuations, comparing numerical results with semiclassical predictions founded on Percival's conjecture. We obtain good $global$ agreement with those predictions, in clear contrast with similar comparisons for other systems found in the literature. The structure of some eigenfunctions, displayed in the quantum Poincar\'e section, provides a clear explanation of the conjecture.

In order to verify Percival's conjecture [J. Phys. B 6,L229 (1973)] we study a planar billiard in its classical and quantum versions. We provide an evaluation of the nearest-neighbor level-spacing distribution for the Cassini oval billiard, taking into account relations with classical results. The statistical behavior of integrable and ergodic systems has been extensively confirmed numerically, but that is not the case for the transition between these two extremes. Our system's classical dynamics undergoes a transition from integrability to chaos by varying a shape parameter. This feature allows us to investigate the spectral fluctuations, comparing numerical results with semiclassical predictions founded on Percival's conjecture. We obtain good global agreement with those predictions, in clear contrast with similar comparisons for other systems found in the literature. The structure of some eigenfunctions, displayed in the quantum Poincaré section, provides a clear explanation of the conjecture.

I. INTRODUCTION
In 1973 Percival conjectured that in the semiclassical limit, the spectrum of a generic dynamical system consists of two parts with strongly contrasting properties: a regular and an irregular part [1]. At the classical level such a system exhibits a mixed dynamics: regular regions dominated by tori and chaotic regions with mixing behavior coexist in the phase space.
In order to characterize a semiclassical spectrum it is advantageous to consider a sensitive fluctuation measure. The probability distribution p(s) of the spacing s between succesive levels is of particular interest because it contains information of the spectrum on its finest scale. In the special case of multidimensional integrable systems, Berry and Tabor [2] showed that the levels are uncorrelated and p(s) is governed by a Poisson distribution. The other special case corresponds to mixing systems where almost all orbits explore densely and chaoticaly the energy surface. In this case, Bohigas et al. [3] conjectured that the fluctuation properties of these spectra can be modeled by the ensemble of random real symmetric matrices [the Gaussian Orthogonal Ensemble (GOE)] [4]. For such systems, where the energy levels display repulsion, p(s) is closely approximated by the Wigner distribution; however, we have used the exact distribution [5] because the differences are meaningful, as we shall see below.
In the generic case, Berry and Robnik [6], based on Percival's conjecture, considered independent sequences of levels associated with each connected regular or irregular classical phase-space region. When only one chaotic region predominates (the situation considered in this article), they provided an expression for the distribution p(s) in terms of the classical fraction ρ cl of regular regions [referred to as the Berry-Robnik distribution (BRD)].
Although the special cases have been extensively confirmed numerically [7][8][9]3,10,11], with some well understood exceptions [12,14], for generic systems numerical calculations give rise to contradictory conclusions. Recently, there has been a number of numerical works [15][16][17] showing that the Brody distribution (BD) [18] gives quite a satisfactory fit globally. The Brody distribution is a one-parameter family of distributions that interpolates between Poisson and Wigner in a simple way; however, it has no semiclassical meaning. On the other hand, in 1994 Prosen and Robnik [19] confirmed numerically semiclassical predictions (the BRD) working on an abstract dynamical system : the standard map on a torus. To agree with this theory, they needed to compute extremely high excited states (around the 30 10 6 ). However, at not so excited states they found good global agreement with the Brody distribution both for the standard map and for the Limacon-Robnik billiard. In 1995 Prosen arrived at the same conclusion working on a two-dimensional semiseparable oscillator [20].
The Brody-like behavior at small spacing is understood in terms of tunneling between classically separated regions of phase space; however, the global agreement with the Brody distribution has no theoretical foundations. On the other hand, the very slow trasition to semiclassical predictions can be explained by the presence of partial barriers in the chaotic regions because the corresponding statistics is not a GOE for finiteh [10,13].
The goal of the present article is to verify that the classical support of eigenfunctions can be clearly identified as regular or chaotic and this classification is only affected by tunneling between classically separated regions of phase space (this effect decays exponentially whenh decreases, as it was pointed out in ref. [6]). We compute the spectral fluctuations of a one-parameter family of planar billiards : the Cassini ovals [21]. The classical dynamics in this billiard is mixed, going from integrability to chaos by varying a single parameter. We have chosen this parameter such that the classical dynamics does not show partial barriers inmerse in the chaotic sea to study the Percival's conjecture through the accuracy of the BRD. Moreover, we study qualitatively the eigenfunctions in phase space to provide additional support to the results.
Our work is organized in the following way. In Sec. II we introduce the classical system. Section III is devoted to the description of the quantum system (its energy spectrum and the corresponding eigenfunctions). In Sec. IV we study the resulting energy level statistic. Finally, Sec. V is devoted to conclusions.

II. THE CASSINI OVAL BILLIARD
Our billiard consists of a free-moving point particle inside a two-dimensional box that bounces off the boundary elastically. The boundary of our billiard system is given by a fourth-order curve, the Cassini oval: where r 1 and r 2 are distances from two foci located at x = ±c and y = 0. In cartesian coordinates it can be cast into the form: (2.1) We have two characteristic lengths. However, the shape of the boundary is defined by the ratio d = a/c (from now on, the shape parameter), which determines the following boundary types: √ 2 < d, the boundary is an oval; 1 < d < √ 2, the boundary is an oval with a neck; d < 1, the boundary becomes unlinked (two ovals separate).
In the present work we investigate d > 1 values. Decreasing the shape parameter, the classical behavior goes from the regular motion (when d → ∞, the boundary is a circle) to the chaotic one (when d → 1). Using the reflection symmetries of the boundary, we consider the motion in the region x > 0, y > 0 (desymmetrized billiard). That is, we study the quarter billiard defined by the boundary (eq. 2.1) for x ≥ 0, y ≥ 0 and the coordinate axes x and y. We study the Cassini oval billiard for two values of d. One of them is d = √ 2 (the value of the parameter for which the neck begins to appear; see Fig. 1(a)) whose form mimics the Bunimovich stadium billiard [22]. Figure 2(a) shows the Poincaré surface at the boundary using Birkhoff coordinates. The coordinate q is related to the arclength coordinate at the boundary where the bounce takes place by q = (arclength)/(perimeter); and p = p ·t/|p| is the fraction of tangential momentum at this point (t being the tangent unit vector to the boundary). Exploiting the time-reversal symmetry we show only the p ≥ 0, 0 ≤ q ≤ 1/4 region. The classical phase space has a bouncing-ball regular region dominated by invariant curves and a chaotic region with unstable short periodic orbits equivalent to those appearing in the Bunimovich stadium billiard. A resonance of winding number 6 defines the last great regular region before chaos begin to appear (an eigenfunction existing on the chain of islands defined by the resonance is shown in Fig.  3 (a)). We have found two very small stable regions corresponding to a stable bifurcation of the unstable bow tie periodic orbit (two small dark dots can be observed in the chaotic region). (Figure 3 (b) shows an eigenfunction existing in that region of the phase space). We have not detected any other regular region embedded in the chaotic sea.
The other shape we have studied (d = 2) is closer to an ellipse (see Fig. 1(b)). In this case the bouncing-ball region of the phase space is greater than before as it can be seen in Fig. 2(b). Moreover, a regular region appears as a thin band for p ∼ > 0.9 values, dominated by whispering gallery trajectories. Phase space is very mixed and many stable islands of very different sizes are interspersed with the chaotic trajectories.
By selecting two regions corresponding to chaotic motion in the phase space, we have calculated diffusion times between them. The results for d = √ 2 are independent of the chosen regions. In the other case (d = 2) this time is strongly dependent on them, and we have obtained diffussion times one order of magnitude greater than those of d = √ 2. This fact is related to partial barriers such as those shown in Fig. 5(a). We have determined ρ cl , the fraction of the phase space that corresponds to regular motion for both the values of the parameter. We have found that ρ cl = 0.172 when d = √ 2 while ρ cl = 0.394 for d = 2.

III. THE QUANTUM BILLIARD
To study the quantum billiard we solve the time-independent Schrödinger equation for one particle inside a twodimensional box D with Dirichlet boundary conditions at the impenetrable walls ∂D: and D corresponds to the surface of the desymmetrized billiard, that is the surface bounded by the Cassini oval (2.1) and the x > 0 and y > 0 semiaxis. So, the Dirichlet boundary condition implies that only odd-odd solutions of the full billiard will be found. We have employed a new technique, the scaling method [23]: this is a very efficient one dimensional method developed to compute eigenvalues and eigenfunctions of quite general planar billiards (for three dimensional billiards this is practically the only available method to obtain high excited states [24]). The great advantage of the method is that all eigenvalues and eigenfunctions in a narrow k interval are computed simultaneously with comparable accuracy, thus avoiding time consuming searches and the posibility of missing some state.
We have calculated the energy levels from the fundamental state up to the 25000th level for d = √ 2 and up to the 10000th for d = 2. Moreover, we took a sequence of 5000 levels between the 62210th and 67210th for both values of d. We obtained the eigenvalues with an average precision of 10 −6 of the mean level spacing for chaotic eigenfunctions. Regular states are only limited by the computer (double) precision.
We have studied the eigenfunctions of the billiard in different regions of the spectrum. In general, it is possible to identify each eigenfunction with a classical region. This identification of the eigenfunctions is more clearly seen in the stellar representation [25]. In it, the Husimi distribution of the normal derivative on the boundary represents the eigenfunctions in the Poincaré section in Birkhoff coordinates. The stellar representation of an eigenfunction is compared directly with the classical Poincaré section. As an example we show some eigenfunctions for d = √ 2 (we have taken the area of the desymmetrized billiard equal to π/4). Fig. 3 (a) shows a linear density plot of the square of a state existing on the chain of islands defined by a resonance of winding number 6 and Fig. 4 (a) shows the same state in stellar representation. Figs. 3 (b) and 4 (b) show a scar [26]  We stress that the states appearing in Figs. 3 (c) and (d) are quasi degenerate. The distance between them is a very small fraction (s = 0.00018) of the mean level spacing. It is clear from Figures 4 (c) and (d) that these states practically do not interact because they exist in different regions of phase space. On the other hand, this example shows us that it is necessary to evaluate the eigenvalues with high accuracy in order to calculate the spectral fluctuations of the system (see the next section).
For d = 2, the qualitative description of the eigenfunctions is equivalent to the previous one. However, there is a significant fraction of localized eigenfunctions existing in the chaotic region. One of them is shown in Fig. 5(b) and (c).

IV. THE ENERGY-LEVEL STATISTIC
In this section we analyze the level spacing distribution of the numerical data described previously. The counting function N (k) gives us the number of levels with wave number below k. Weyl's formula with border, angle and curvature corrections [8] provides a good estimate for the smooth part < N (k) > of the counting function.
In order to verify that no levels had been lost we have compared < N (k) > with a smoothed version of N (k). Defining a new sequence by K n ≡< N (k n ) > ( where k n belongs to the original sequence of wave numbers ) we take into account of the "unfolding" procedure by which a unit mean spacing is given to the series of levels [19].
Following ref. [17] in the analysis of the data, we use the cumulative level spacing distribution W (s) = s 0 p(y)dy rather than p(s) because it is numerically easier to evaluate. We fit numerical curves with the analytical expression for the Brody family of cumulative spacing distributions: with b = (Γ((β + 2)/(β + 1))) β+1 and the theoretical Berry-Robnik distribution: λ p GOE (y)dy, p GOE (y) being the exact GOE spacing distribution, and W GOE = 1 + dQ GOE /d((1 − ρ)s). Q GOE and W GOE are tabulated as Ψ and F (taking as mean density 1 − ρ) respectively in reference [5], for instance. This exact GOE evaluation allows us to distinguish deviations of numerical data from theory without including the difference between Wigner and exact GOE formulas (which are of approximately the same order; this fact can be verified in Fig. 6 (a) which shows the difference between the BRD using Wigner surmise and the exact GOE results).
Deviations of numerical data from best-fitting curves can be better seen with a transformation defined by: which has constant statistical error over all s. Furthermore, if we plot U (W ) versus W B β we will have an equally spaced distribution of points on the abscissa (see Fig. 6, for example).
We have evaluated )/N so that we could find the optimal values of β and ρ BR (ρ BR is the resulting best fitting value for the fraction ρ of regular levels employed as a free parameter in order to find the lowest χ 2 ). Results from fittings can be seen in Table I. It is clear that, for d = √ 2, BRD curve fits data much better than BD curve. Taking the sequence of 25000 first levels we found that χ 2 for the BRD was approximately five times lower than for the BD. Even for the small s range we could verify a better agreement between data and BRD formula than for the BD one, although deviations due to tunnelling can be clearly seen in Figure 6 (a). For large values of s, differences between data and BRD could evidence discrepancies with respect to GOE behaviour. Taking the sequence of 5000 levels between number 62210 and 67210 we got a χ 2 for BRD that was 25 times lower than for BD. This can be checked with Fig. 6 (b), where the region of small s still shows some deviations from best fitting BRD but the range where this happens has become very narrow. In fact, for all values of s, the agreement is excellent and we can say that BRD is working perfectly well at these not so high energy levels.
In the case of d = 2, for the first 10000 levels, a same order χ 2 was found for the two distributions although was better for the BRD than for the BD. In Fig. 7 (a) we can appreciate that BD fits numerical data in an acceptable way for small s values. This is due to tunnelling effects that persist in a wider s range than for the d = √ 2 case. A more complicated structure of classical phase space makes tunnelling processes more significant. However, for greater s BD no longer follows numerical data, so there is no global agreement with it. As in the previous case, we take 5000 levels between 62210 and 67210. For these energies, χ 2 for BD is twice the value for BRD. In Fig. 7 (b) we can see that discrepancies in the small s region have reduced, and, though not as clear as in d = √ 2 case, numerical data is well adjusted by BRD globally.
Finally, we investigated the behaviour of data while going from low to high energies. In order to do so we took first 3000 levels for both of the shapes; then three stripes of levels, the first, second and third 8000 levels for d = √ 2 and the first, second and third 4000 levels for d = 2. Results of fittings can be seen in Table II. It is clear from Figures 8  and 9, that there is no transition from Brody to Berry-Robnik regime, neither for d = √ 2, nor for d = 2.

V. SUMMARY AND CONCLUSIONS
To verify Percival's conjecture, we have studied the quantum version of a billiard, depending on one shape parameter d. This system shows mixed classical dynamics, going from integrability (d → ∞) to chaos (d → 1) as the parameter is varied.
On the one hand, the Husimi distribution of the eigenfuntions that we have obtained clearly displays the classical structure of the phase space. Though mixings among regular wavefunctions and irregular ones are expected (based on the abscence of degeneracies in this one-parameter system that has been desymmetrized) they seem to happen only when the energy difference is surprisingly small ( even for levels as low as N ≃ 2500 such as the ones exemplified in Fig. 3 (c) and (d)). So, we can say that Percival's conjecture is effectively working. Mixed functions are exceptions mainly originated in the states whose Husimi distributions localize on the last KAM tori.
On the other hand, the level spacing statistics was fitted by two distributions depending on one parameter: the semiclassical Berry-Robnik distribution (BRD), which is founded on Percival's conjecture, and the Brody distribution (BD). We have found, in all cases we have analyzed that the BRD is the best one. Moreover, as expected, fits of the BRD are better for increasing values of the wave number k; however, discrepancies due to tunnelling effects are found for small values of s, even for the best-fitting case (see Fig. 6 (b) -This subject is currently under study at the moment). For d = √ 2, the parameters ρ of the BRD corresponding to best fits, do not show significant differences from the fraction of the regular region of the classical phase space. In the case of d = 2, the best ρ's, are systematically greater than the classical value.
Our calculation of diffusion times shows that the cassini billiard with d = √ 2 has only one chaotic region. There are not partial barriers dividing chaotic regions of comparable sizes. So we expect that the chaotic component of the statistic is given by a single GOE distribution. In the case d = 2 we have determined partial barriers between a main chaotic region and a little one near the whispering gallery region. The size of the latter is so little that it enlarges (for low energies) the regular region of the phase space rather than contributes as an independent GOE. This fact can be seen in Fig. 5(b) and (c), where we show an eigenfunction localized in this small classical chaotic region but having regular characteristics. Moreover, there are small regular islands immerse in the chaotic sea that are more relevant for d = 2 (see Fig. 2 ) than for d = √ 2. They produce quantum localization in their chaotic neighborhood (for instance, see the Poincaré surface section Fig. 2 (a) and the Husimi distribution Fig. 4 (b) ) and consequently the resulting ρ of the fits, overestimates the fraction corresponding to regular motion in the classical phase space. We stress that, in all cases, the BD has significant differences from the numerically calculated distribution. So we do not observe the BD to BRD transition that was seen in other systems.

ACKNOWLEDGMENTS
We would like to thank Marcos Saraceno for useful suggestions on the original manuscript. This work was partially supported by CONICET and UBACyT. . We can observe two partial barriers limiting flux and defining two small chaotic regions (which are filled more densely than the major chaotic part of phase space). b) Husimi distribution of an eigenfunction localized over one of these small regions. c) Same eigenfunction in configuration space.
FIG. 6. Differences between numerical U and U BR taking the best-fitting value ρ BR are displayed with a solid (fluctuating) line for d = √ 2. Dotted lines that follow the solid one represent the δU uncertainty band. a) N between 1 and 25000 and b) N between 62210 and 67210. In a), difference between U BR calculated with Wigner surmise and with exact GOE results is also plotted with a dashed line. We can see that deviations of numerical data from this curve are of the same order.