Fluctuating Number of Energy Levels in Mixed-Type Lemon Billiards

: In this paper, the ﬂuctuation properties of the number of energy levels (mode ﬂuctuation) are studied in the mixed-type lemon billiards at high lying energies. The boundary of the lemon billiards is deﬁned by the intersection of two circles of equal unit radius with the distance 2 B between the centers, as introduced by Heller and Tomsovic. In this paper, the case of two billiards, deﬁned by B = 0.1953, 0.083, is studied. It is shown that the ﬂuctuation of the number of energy levels follows the Gaussian distribution quite accurately, even though the relative fraction of the chaotic part of the phase space is only 0.28 and 0.16, respectively. The theoretical description of spectral ﬂuctuations in the Berry–Robnik picture is discussed. Also, the (golden mean) integrable rectangular billiard is studied and an almost Gaussian distribution is obtained, in contrast to theory expectations. However, the variance as a function of energy, E , behaves as √ E , in agreement with the theoretical prediction by Steiner.


Introduction
The purpose of this paper is to analyze the energy spectra of two characteristic complex mixed-type lemon billiards within the scope of quantum chaos. The boundary of the lemon billiards is defined by the intersection of two circles of equal unit radius with the distance 2B between the centers, as introduced by Heller and Tomsovic in 1993 [1,2]. The present study represents a continuation of our recent paper [3] on the classical and quantum ergodic billiard (B = 0.5) with strong stickiness effects, from the family of lemon billiards, as well as on three simple mixed-type lemon billiards with only one regular region, surrounded by a uniform chaotic sea without stickiness regions, namely, with the shape parameters B = 0.42, 0.55, 0.6 [4].
In the present paper, two lemon billiards with the shape parameters B = 0.1953, 0.083 are studied. These lemon billiards are mixed-type billiards with several independent regular regions embedded in a chaotic sea with no significant stickiness regions, which serve as examples of systems with a complex divided phase space. These lemon billiards were selected by the criterion of the maximally complex chaotic component generated by a single chaotic orbit. The discovery of the present and past physically and dynamically different and interesting lemon billiards has only been made possible thanks to the recent extensive analysis of Lozej [2]. The entire family of classical lemon billiards for a dense set of about 4000 values of B ∈ [0.01, 0.99975] (in steps of dB = 0.00025) was systematically analyzed as for the corresponding phase space structure and stickiness effects. It must be emphasized that although all the lemon billiards belong to the same family of billiards as for the mathematical definition, individually, the lemon billiards have quite different, in fact, very rich, dynamical properties, important in the classical and quantum context.
A general introduction to the subjects in quantum chaos, related to this study, can be found in [3]. Let us also mention the books by Stöckmann [5] and Haake [6] on general quantum chaos and the recent reviews [7,8] on stationary quantum chaos in generic (mixedtype) systems.
The main purpose of the present paper is to analyze the two selected quantum lemon billiards of B = 0.1953 and 0.083, with the goal to study the energy spectra, while the structure of the Poincaré-Husimi functions in the phase space, the separation of the regular eigenstates and chaotic eigenstates, as well as the localization properties of the chaotic eigenstates and their statistics will be treated separately [9].
The main results are as follows. The energy spectra are calculated by the scaling method of Vergini-Saraceno [10] in two versions, one based on the plane waves and the other one based on the circular waves (Bessel functions for the radial part and trigonometric functions for the angular part). This is done for high-lying eigenstates with the wavenumber k (in the specific units) k = 2880 for up to 300,000 consecutive levels for each of the four symmetry classes (odd-odd, odd-even, even-odd, and even-even). The energy of the level at k is E = k 2 . It turns out that about 0.1% of the levels are lost for technical reasons, which is a known and experienced fact in numerically calculating billiard spectra with the scaling method, while the accuracy of individual levels is better than 1% of the mean energy level spacing. The spectral statistics are found to be stable with respect to these losses. The cumulative (integrated) energy level density (spectral staircase function) N (k) is well described by the Weyl formula (with the leading area term and the perimeter term) W(k) if there are no missing levels. Then the fluctuations of the actual staircase function N (k) around the Weyl function, the difference R(k) = N (k) − W(k), and the R(k) distribution are studied. In the literature, R(k) is called mode fluctuation [11][12][13][14]. Due to the lost levels, this difference has a drift to negative values and fluctuates around the mean value. In order to separate the drift and the fluctuations, the quantity w(k) = R − m(R), where m(R) is the local average of R(k) over 100 consecutive levels, is investigated. Then, the distribution of w(k) for about 300,000 levels of each symmetry class separately, starting at k = k 0 = 2880, within the interval approximately k ∈ [2880, ≈ 3700] is studied. One finds that the distribution of w follows a Gaussian distribution fairly well, which is a surprising result as the two billiards have the relative fraction of the chaotic phase space of only 0.28 and 0.16, respectively. For comparison, also the entirely regular, integrable, case of the maximally irrational rectangular billiard is investigated and it is surprisingly found an almost Gaussian distribution for w, but with the variance rising linearly with k, as generally predicted by the theory of Steiner [11][12][13][14], while the distribution itself in integrable systems is predicted to be nonuniversal, varying from case to case, which here is not confirmed. The validity of the theoretical predictions has also been checked in experiments with superconducting microwave billiards [15]. Finally, the level spacing distribution for the two lemon billiards is presented which provides good agreement with the Berry-Robnik-Brody distribution.
The paper is organized as follows. In Section 2, the definition are given and classical dynamical properties of the lemon billiards are examined. In Section 3, the analysis of the fluctuation quantity w(k) is performed. In Section 4, the statistical analysis of the entire energy spectra is made by calculating the level spacing distribution, P(S). In Section 5, the results are discussed and conclusions are presented.

The Lemon Billiards and Their Phase Portraits
The family of lemon billiards was introduced by Heller and Tomsovic in 1993 [1] and has been studied in a number of papers [16][17][18][19][20], most recently, in [2][3][4]21]. and in the recent studies [3,4] by us. The lemon billiard boundary is defined by the intersection of two circles of equal unit radius with the distance between the centers, 2B, being less than the diameters and B ∈ (0, 1), and is given by the following implicit equations in Cartesian coordinates: As usual, the canonical variables are used to specify the location, s, and the momentum component, p, on the boundary at the collision point. Namely, the arc length, s, is counted in the mathematical positive sense (counterclockwise) from the point (x, y) = (0, − √ 1 − B 2 ) as the origin, while p is equal to the sine of the reflection angle θ; thus p = sin θ ∈ [−1, 1], as θ ∈ [−π/2, π/2]. The bounce map (s, p) ⇒ (s , p ) is area preserving as in all billiard systems [22]. Due to the two kinks, the Lazutkin invariant tori (related to the boundary glancing orbits) do not exist. The period-2 orbit connecting the centers of the two circular arcs at the positions (1 − B, 0) and (−1 + B, 0) is always stable (and therefore surrounded by a regular island) except for the case B = 1/2, where it is a marginally unstable periodic orbit, the case being ergodic and treated by us earlier [3].
The circumference of the entire billiard, L, is given by: The area A of the billiard is: The structure of the phase space is shown in Figure 1 for the lemon billiard of B = 0.1953, with L = 5.4969. The relative fraction of the area of the chaotic component of the bounce map is χ c = 0.3585, while the relative fraction of the phase space volume of the same chaotic component is ρ 2 = 0.2804, the Berry-Robnik parameter. Three independent regular island chains are clearly visible, the largest one around the period-2 orbit, which is densely covered by the invariant tori, with no visible thin chaotic layers inside. Let us denote the largest island chain by L, the second largest one by M, and the smallest one by S. The relative phase space volume of all three regular regions taken together is ρ 1 = 1 − ρ 2 = 0.7196. The chaotic sea is quite uniform, with no significant stickiness regions, and is generated by a single chaotic orbit.
The structure of the phase space, as shown in Figure 2 for the lemon billiard of B = 0.083, is more complex. Here, L = 5.9508. The relative fraction of the area of the chaotic component of the bounce map is χ c = 0.2168, while the relative fraction of the phase space volume of the same chaotic component is ρ 2 = 0.1617. Thus, the relative phase space volume fraction of the complementary regular regions is ρ 1 = 1 − ρ 2 = 0.8383. In this case also, the chaotic sea is rather uniform, with no significant stickiness regions, and is generated by a single chaotic orbit, creating a very complex structure, perhaps mimicking stickiness in some thin regions. Both billiards are clearly of the Kolmogorov-Arnold-Moser (KAM) type, generic systems; examples of stochastic transition were studied already in [23].
One can conclude that the two cases B = 0.1953, 0.083 are interesting to verify the Berry-Robnik picture of quantum billiards [24], including the possible quantum localization of chaotic eigenstates, leading to the Berry-Robnik-Brody level spacing distribution and the universal statistical properties of the localization measures [4,9], as there are no significant stickiness effects, based on the results of the analysis of the recurrence time statistics in [2], unlike in the ergodic case, B = 0.5, studied in [3].

The Energy Spectra and the Fluctuation of the Number of Energy Levels
Let us now turn to the quantum billiard B described by the stationary Schrödinger equation in the chosen units (h 2 /2m = 1) given by the Helmholtz equation: with the Dirichlet boundary conditions ψ| ∂B = 0, and the energy is E = k 2 . Hereh is the reduced Planck constant and m is the mass of the particle.
The mean number of energy levels W(E) below E = k 2 is determined quite accurately, especially at large energies, asymptotically exact, by the celebrated Weyl formula (with perimeter corrections) using the Dirichlet boundary conditions, namely: where "c.c." stays for small constants, determined by the corners and the curvature of the billiard boundary, which differentially play no role. Thus, the mean density of levels d(E) = dW/dE is equal to: The numerical method used here to solve the Helmholtz equation is based on the Heller's plane wave decomposition method and the Vergini-Saraceno scaling method [10,25]. Both versions of the Vergini-Saraceno method are implemented, namely, the one, based on plane waves, and the other one, based on the circular waves, and the same results are obtained within an accuracy of 0.1% of the mean level spacing. As mentioned, typically at most about 0.1% of the levels are lost. The numerical accuracy was checked by the convergence test, by varying the method's parameters (the number of basis waves and the number of nodes on the boundary).
The billiard considered has two reflection symmetries; thus, the eigenstates have four symmetry classes: odd-odd, odd-even, even-odd, and even-even. For the purpose of analyzing the spectral statistics and the wave functions, let us consider only the quarter billiard. In this case, the Weyl formula for the four symmetry classes reads: whereL is defined for each of the above-defined symmetry classes as follows: Note that summing up the four contributions in Equation (7) with Equation (8), one obtains the Weyl formula for the entire spectrum Equation (5). For each symmetry class for each of the two billiards, about 300,000 levels were calculated and the difference, R(k) = N (k) −W(k), was studied between the staircase function, N (k), and the Weyl function,W(k). Figure 3 shows the results for the B = 0.1953 billiard. It is clear that R(k) decreases linearly with k due to the numerical loss of the levels. Typically, about 0.1% of the levels are missing. In order to study the fluctuation properties, one needs to subtract the local mean value, m(R(k)), obtained by averaging over 100 consecutive levels, and then analyze the resulting fluctuation quantity, w(k) = R(k) − m(R(k)). The agreement with the Gaussian distribution is extremely good which is surprising, as the billiard is predominantly regular, with only 28% of the chaotic component. Steiner's theory [11][12][13][14] predicts a Gaussian distribution only for entirely chaotic (ergodic) systems, while the distribution of the fluctuation quantity, w, in integrable systems is expected to be nonuniversal. In the mixed-type case, the distribution of w should be in between, which here is not the case. Moreover, even for the rectangle as an example of the integrable systems, an (almost) Gaussian distribution is found (see below).    The local averaging procedure over 100 consecutive levels is, of course, somewhat arbitrary. Therefore, the results were checked for 50, 200, 400 levels. In all cases, the distribution was found to represent Gaussian extremely well. The standard deviation σ at each level number keeps the same value for all parities, but increases slowly with the number of levels such that for 50, 100, 200, 400 for B = 0.1953, one gets σ ≈ 2.08, 2.63, 2.88, 2.92, respectively. For the B = 0.083 billiard, the corresponding values are: σ ≈ 2.22, 2.92, 3.50, 3.52.
In order to explore an integrable, entirely regular system, an example of the maximally irrational rectangle billiard was studied whose aspect ratio was taken the golden mean g = (1 + √ 5)/2 ≈ 1.61803. The spectrum in this case is known analytically, and is exact: E ln = l 2 /g 2 + n 2 , where l and n are two positive integers. Figure 5 shows the fluctuation quantity w(k), which has no drift, because the spectrum is exact (no lost levels) over a wide range of k. The corresponding distributions at various wavenumber k intervals starting at k 0 are shown in Figure 6. Each histogram comprises about 100,000 objects. They are surprisingly close to a Gaussian distribution, contrary to the theoretical expectation, and thus, the distribution is just close to the case of ergodic chaotic systems. Therefore, the distribution of the number of energy levels (distribution of mode fluctuations) is not a good criterion to distinguish ergodic chaotic systems from the regular integrable systems. In Table 1, the data for the mean, µ, standard deviation, σ, skewness, and kurtosis are given.  Table 1, the values of the mean µ, standard deviation σ, skewness, and kurtosis are given. The fitting (red) curve is the Gaussian distribution with the same µ and σ as those obtained. However, the variance as a function of k is linear for all integrable systems; it is thus universal in the class of integrable systems, see Figure 7. Table 1. Parameters of the distribution of Figure 6 for the wavenumber k-intervals starting at k 0 : the mean µ, standard deviation σ, skewness, and kurtosis. Figure 6 k In Figure 7, the standard deviation of the distributions of w is plotted as a function of k, which clearly accurately follows the prediction by Steiner [11][12][13][14], namely the standard deviation σ rises as √ k, shown also in the log-log plot, and the variance is a linear function of k. Therefore, the dependence of the variance on k is a good signature of chaos, unlike the distribution itself: in ergodic chaotic systems, the standard deviation behaves as log k.

Level Spacing Distribution for the Entire Spectrum
One of the most important statistical measures of the (unfolded) energy spectra is the level spacing distribution, P(S). For integrable systems, one gets Poissonian statistics and P P (S) = exp(−S), while for classically ergodic (fully chaotic) systems the Gaussian orthogonal ensemble (GOE) of random matrix theory applies. The Wigner distribution (Wigner surmise) is 2-dimensional GOE distribution and is a very good approximation for the GOE level spacing distribution (infinite-dimensional), There is a general useful relationship, namely, one using the gap probability, E (S), being the probability of having no level on an arbitrary interval of length S: the level spacing distribution P(S) is, in general, equal to the second derivative of the gap probability, P(S) = d 2 E (S)/dS 2 .
For Poisson statistics: E P (S) = exp(−S), while for the Wigner distribution, one finds: For mixed-type systems, there is typically one dominant chaotic component with the relative density of levels ρ 2 (equal to the relative fraction of the chaotic phase space volume), while the complement is typically a regular component of relative density, ρ 1 = 1 − ρ 2 . If the regular and chaotic levels superimpose statistically independent of each other, then obviously, the gap probability factorizes: and therefore, in this case, the level spacing distribution is given by the Berry-Robnik (BR) formula [24]: The above statements are true provided the Heisenberg time is larger than any classical transport time in the system [8]. (The Heisenberg time is defined as 2πh d(E), where d(E) is the mean energy level density, also the reciprocal mean energy level spacing.) If this is not the case, the chaotic eigenstates can be quantum (or dynamically) localized, which implies localized chaotic Poincaré-Husimi functions in the phase space. The level spacing distribution for such localized chaotic eigenstates becomes (approximately) the known Brody distribution [26,27]: where by the normalization of the total probability and the first moment, one gets: where Γ(x) is the Gamma function. The Brody distribution interpolates the exponential and Wigner distribution as β goes from zero to one. The important feature of the Brody distribution is the fractional level repulsion effect, meaning the power law at small S, P(S) ∝ S β . The corresponding gap probability is: where γ = Γ β+2 β+1 , and Q(a, x) is the incomplete Gamma function: Here, the only parameter is β, the level repulsion exponent in (13), which measures the degree of localization of the chaotic eigenstates: if the localization is maximally strong, the eigenstates practically do not overlap in the phase space (of the Wigner functions), and one finds β = 0 and a Poissonian distribution, while in the case of maximal extendedness (no localization), one finds β = 1 and the GOE statistics of levels applies. Thus, by replacing E W (S) by E B (S), the Berry-Robnik-Brody (BRB) distribution is obtained which generalizes the BR distribution (12) such that the localization effects in chaotic eigenstates are included [28]. Note that in the semiclassical limit,h → 0 or k → ∞, the Heisenberg time becomes arbitrarily large (larger than any classical transport time), the localization effects of chaotic eigenstates disappear, and the BRB distribution becomes the BR distribution. However, the theoretical derivation of the Brody distribution for the localized chaotic states remains an important open problem. Furthermore, while the local behavior at small S, described by the power law P(S) ∝ S β , is certainly correct, the global feature of the Brody distribution is surely approximate, although empirically is well founded.
In Figure 8, the present study, the classical transport time of the billiards is very short; therefore, one expects β ≈ 1, and the level spacing distribution is almost BR one (12). Thus, in the level statistics, one does not detect large localization effects. For the spectral unfolding procedure, the Weyl formula, Equation (5), is used, which at high energies is quite accurate. Figure 8, shows the level spacing distributions P(S) for the two billiards-one of B = 0.1953 (left column), and another one of B = 0.083 (right column)-for the spectral stretches each about 20,000 levels long, starting at k 0 = 2880, for each parity: even-even even-odd, odd-even and odd-odd and for all four parities together (about 80,000 levels), each of them along with the best fitting BRB distribution.
As one can see, in the case of B = 0.1953, Figure 8a,c,e,g,i, the β parameter is very close to one, while the parameter ρ 1 is close to the classical value, ρ 1 = 1 − ρ 2 = 0.7196, being the relative fraction of the volume of the regular part of the phase space (see Section 2). Thus, the system exhibits the Berry-Robnik picture with weak localization effects in the chaotic part of the energy spectrum. P(S) is well fit by the BRB distribution. The picture is very similar for other lower values of k 0 (not shown).
For B = 0.083, one observes a substantial variation of both the β and ρ 1 parameters, the latter one being far from the classical value ρ 1 = 0.8383. This is certainly due to the complexity of the phase portrait shown in Figure 2, where many structural features are quantum mechanically not yet resolved, even at such high lying energies E = k 2 , starting at k 0 = 2880. This will be studied in more detail in the forthcoming paper by us [9], where the Poincaré-Husimi functions in the phase space are analyzed in a similar way as in [4].

Conclusions
In this paper, the statistical properties of the oscillations of the cumulative spectral staircase function (the integrated density of energy levels) around the corresponding mean value were studied, in order to compare the function obtained with the theoretical predictions of Steiner [11][12][13][14] for the fully chaotic and integrable (regular) systems. In billiards, the mean behavior is asymptotically exactly described by the Weyl formula (5). In the case of the integrable maximally irrational rectangle billiard, almost a Gaussian distribution of the fluctuations was surprisingly observed, in contrary to the expectations of Steiner's theory. Nevertheless, the standard deviation as a function of the wavenumber k was found to be proportional to √ k, precisely in agreement with Steiner's prediction for all integrable systems.
In the two mixed-type lemon billiards, B = 0.1953 and 0.083, where regular and chaotic regions coexist in the classical phase space, the Gaussian distribution is found, which is surprising, as the theory predicts a Gaussian distribution only for the fully chaotic (ergodic) systems. Thus, in these two cases, one observes that there is very little difference between the mixed-type systems and the fully chaotic systems. Moreover, this is even more surprising as in the two billiards studied here, the fraction of the chaotic part of the phase space was only 0.28 and 0.16, respectively. This implies that the distribution of the fluctuations of the spectral staircase functions around the mean behavior is not a very significant criterion for quantum chaos. This conclusion is corroborated by the result obtained for the integrable rectangle billiard.
The level spacing distribution of the energy spectra at high lying levels was also explored, starting at k 0 = 2880, and the Berry-Robnik-Brody distribution was found in both lemon billiard cases. In the case B = 0.1953, the results were in agreement with the Berry-Robnik picture [8,24], showing that the localization of the chaotic eigenstates is very weak, the level repulsion parameter β is close to one, and the quantum Berry-Robnik parameter ρ 2 is close to the classical one. On the other hand, in the lemon billiard B = 0.083, one observed relatively strong localization: β is substantially lower than one and varies widely over the four parities. Likewise, the quantum ρ 1 is substantially smaller than the classical value 0.8383 and varies widely over the four parities, which is certainly due to the complexity of the classical phase space, as the fine structure of the classical phase space is not yet resolved by the Poincaré-Husimi functions. This analysis will be the topics of forthcoming paper by us [9], using the approach of [4].