Quantizing neutrino billiards: an expanded boundary integral method

With the pioneering fabrication of graphene the field of relativistic quantum chaos emerged. We will focus on the spectral properties of massless spin-1/2 particles confined in a bounded two-dimensional region, named neutrino billiards by Berry and Mondragon in 1987. A commonly used method for the determination of the eigenvalues is based on a boundary integral equation originating from Green’s theorem. Yet, in the quantization one might face problems similar to those occurring for non-relativistic quantum billiards. Especially in cases where the eigenvalue spectrum contains near degeneracies the identification of complete sequences of eigenvalues might be extremely elaborate, if not unfeasible. We propose an expanded boundary integral method, which yields complete eigenvalue sequences with a considerably lower numerical effort than the standard one. Actually, it corresponds to an extension of the method introduced in Veble et al (2007 New J. Phys. 9 15) to relativistic quantum billiards. To demonstrate its validity and its superior efficiency compared to the standard method, we apply both methods to a circular billiard of which the eigenvalues are known analytically and exhibit near degeneracies. Finally, we employ it for the quantization of a neutrino billiard with a hole, of which the spectrum contains many close lying levels and exhibits unusual fluctuation properties.


Introduction
The search for manifestations of classical chaos within the field of quantum chaos in properties of the associated quantum system has been a focus of research during the last 5 decades [1][2][3]. One of the concepts of quantum chaos, namely the understanding of the features of the classical dynamics in terms of the spectral properties of the corresponding quantum system, like nuclei, atoms, molecules, quantum wires, and quantum dots or other complex systems [4][5][6][7][8][9][10], has been elaborated extensively. According to the 'Bohigas-Giannoni-Schmit' conjecture [11][12][13][14] for chaotic systems the spectral fluctuation properties of a generic chaotic system coincide with those of random matrices from the Gaussian orthogonal ensemble (GOE), if time-reversal invariance is conserved, and from the Gaussian unitary ensemble (GUE), if it is violated. For integrable systems they have been shown to coincide with those of random numbers from a Poisson process [15]. Billiards are well suited for the study of the dependence of the features of a quantum system on the degree of chaoticity of its classical dynamics since it depends only on their shapes. Furthermore, the corresponding non-relativistic quantum mechanical system, referred to as quantum billiard in the sequel, can be investigated experimentally using flat, cylindrical microwave resonators of corresponding shapes. In such systems signatures of classical chaos were observed in the fluctuation properties of their resonance frequency spectra, in the statistical properties of their electric field distributions, in the fluctuation properties of the scattering matrix elements of chaotic scattering processes, also in systems, where time-reversal invariance is broken [3,[16][17][18][19].
Most of the works on quantum chaos focus on systems described by the non-relativistic Schrödinger equation. Only recently a new field of interest emerged, namely relativistic quantum chaos [20,21]. These studies were triggered by the pioneering fabrication of graphene [22]. In the low-energy region it exhibits extraordinary properties attributed to the conically shaped band structure near the Fermi energy. Consequently, in that energy region the electronic properties of graphene are described by a Dirac equation so that graphene features relativistic phenomena [23,24]. Accordingly, the question arose whether the conclusions drawn on properties featured by generic non-relativistic quantum systems with a classically chaotic counterpart also apply to those of graphene dots, frequently referred to as graphene billiards; see [20,21,25] for an overview on numerical and experimental studies. With this focus on the spectral properties of graphene billiards in the relativistic energy region interest in those of a massless spin-1/2 paricle confined in a two-dimensional domain -named neutrino billiard in the seminal work [26] of 1987, where neutrinos still were considered to be massless -re-emerged. Confinement of the particle to the billiard domain is achieved by ensuring a vanishing outward probability current and realized by imposing appropriate boundary conditions. The Dirichlet boundary condition, commonly used for non-relativistic quantum billiards, would comply with this requirement, yet it is inadequate, because it destroys the self-adjointness of the associated Dirac Hamiltonian. One possible boundary condition was proposed in [26] implying the need of the extension of the, by now well developed, numerical methods employed in quantum billiards to neutrino billiards. There already exist numerous works on the study of neutrino billiards with shapes of classically integrable billiards or resulting from a conformal map of the unit circle. The first ones were actually presented in [26]. It was found that for circular neutrino billiards the spectral properties follow Poisson statistics whereas for one with the shape of an Africa billiard with a fully chaotic classical dynamics they do not coincide with those of random matrices from the GOE as would be expected, but instead exhibit GUE statistics, i.e. they behave like chaotic systems with violated time-reversal invariance. In [26] the eigenvalues were determined by numerically solving the boundary integral equation resulting from Green's theorem for massless relativistic spin-1/2 particles. It has the enormous advantage compared to the corresponding equation for quantum billiards, that its integrand does not exhibit singularities. Yet, its numerical evaluation may be accompanied by problems similar to those encountered in quantum billiards, e.g. for too closely-lying eigenvalues, inner corners and holes in the billiard domain like the Monza billiard [27]. The purpose of the present article is to introduce a procedure for the quantization of such billiards which incorporates an expanded boundary integral method of the type proposed in [27].
In section 2 we will introduce the Dirac equation of a massless spin-1/2 particle in the most commonly used coordinate systems, in section 3 the boundary conditions defining a neutrino billiard are specified and its salient features are reviewed, in section 4 the boundary integral is given and compared to that of quantum billiards with mixed boundary conditions, in section 5 we introduce the boundary integral method, in section 6 the expanded boundary integral method is derived, tested in circular neutrino billiards exhibiting nearly degenerate eigenvalues and finally applied to the Monza billiard, which comprises a hole.
2. The Dirac equation for a free, massless spin-1/2 particle in general coordinates For a massless spin-1/2 particle moving freely in the (x, y) plane the Dirac Hamiltonian is given as and , x y x y y y = -- Its solutions may be written explicitly in terms of plane waves where f 0 is the angle between k and the x-axis.
We consider a massless spin-1/2 particle which is confined to a bounded two-dimensional region, referred to as neutrino billiard in the original work [26]. Generally, the boundary of billiards with polygonal shape is parametrized in cartesian coordinates, whereas curved shapes resulting from a conformal mapping of the circle are prevalently defined in terms of polar coordinates 0r1 and 0f2π by a polynomial in z in the complex plane with c l real or complex coefficients, r<1 parametrizing the interior and r = 1 the boundary. The coordinate transformation from cartesian coordinates (x, y) to w z re i = f ( )yields for the gradient in the complex plane where the star denotes complex conjugation and w z w z z The simplest example is a neutrino billiard with circular shape for which polar coordinates are most suitable, or in the complex-plane representation w(z)=z and w′(z)=1. Solutions of equation (7) are again given in terms of plane waves, which with k R kR cos 0 q q = -· ( ) can be expanded in polar coordinates as [28] where J m (x) denotes the mth Bessel function of the first kind. An ansatz of the form equation (5)  . In these coordinates z f r e i = f ( ) with f (r=r 0 ) defining the boundary of the billiard for some value r=r 0 and Another commonly used parametrization for convex billiards is with the gradient given as and the boundary again defined by some value r=r 0 . We are eventually interested in solutions of the Dirac equation in bounded domains which contain the origin and in view of equation (8) Inserting this ansatz into equation (7) and using The same set of solutions is obtained when choosing the parametrization (9) or (10).

Neutrino billiards
The expansions equation (18) provide the solutions of the Dirac equation for a free spin-1/2 particle in an unbounded system. To obtain those for a massless spin-1/2 particle which is confined to a bounded twodimensional domain of which the boundary is given at r=1 by w z e , 0, 2 in the sequel referred to as w(f), we proceed as described in [26]. The appropriate boundary condition is obtained from the requirement that the Hamiltonian of a closed system should preserve self-adjointness implying that it is Hermitian, i.e. the eigenvalues should be real. This condition ensures the impenetrability of the walls, that is, that there is no outward flux is determined from the derivative dr/df or in complex-plane . Furthermore, we would like to emphasize, that this choice of boundary conditions is not the only one guaranteeing self-adjointness of the Dirac Hamiltonian and zero outgoing current [30]. Further information on boundary conditions for the confinement of relativistic particles to a bounded region, called 'MIT bag model', can be found on pages 108-129 of [31].
3.1. Implications of geometric symmetries of a neutrino billiard for its eigenstates At this point a remark concerning the impact of symmetry properties of the billiard shape on those of the wave functions is necessary. In a quantum billiard a mirror symmetry of its shape implies that wave functions are either symmetric or antisymmetric with respect to the symmetry axis, that is, follow Neumann, respectively Dirichlet boundary conditions along this axis. Consequently numerical efforts and the computation time required for the determination of its eigenstates can be drastically reduced by separating the boundary integral equation according to these boundary conditions into two independent ones. This, however, is not the case for a neutrino billiard of which the shape has a mirror symmetry, e.g. with respect to the x-or y-axis, that is, is invariant under reflections y→−y at the x-axis, respectively, x→−x at the y axis.
, and x y x y x y x y , , , , respectively, in accordance with the boundary condition equation (20), implying real coefficients a l in equation (18).
If, on the other hand, the billiard shape has an n-fold rotational symmetry, then the Hamiltonian can be brought to a block diagonal structure in its eigenrepresentation with each of the n blocks corresponding to an invariance with respect to a rotation by m which is fulfilled for the eigenfunctions equation (27) if it is for those of H NB itself. Note that, if the boundary shape has mirror symmetries with respect to two perpendicular axes, it has a twofold symmetry so that for each of its eigenfunctions the component ψ 1 is symmetric and ψ 2 is antisymmetric with respect to a rotation by π, yielding for the solutions equation We conclude this section with a remark on time-reversal invariance, i.e. on the invariance of H NB under T K i y s =ˆˆwith K denoting complex conjugation [2,32]. The transformed eigenfunction does not fulfill the boundary condition, that is, the Dirac Hamiltonian of neutrino billiards is not time-reversal invariant. Consequently, the spectral fluctuation properties of a typical neutrino billiard with the shape of a chaotic billiard are described by the GUE, if it does not exhibit any geometric symmetries.

Green's theorem and a boundary integral for neutrino billiards
Only in a few cases, namely for certain shapes corresponding to classically integrable billiards, like a circle neutrino billiard [26], or for billiards of which the boundary results from a conformal transformation of the circle [33], solutions can be found by imposing the boundary condition equation (20) on an ansatz of the form equation (18), e.g. in terms of the eigenfunctions of the circular neutrino billiard [33,34] similar to [35]. Generally a plane-wave decomposition [36] or the scaling method [37] based on the ansatz equation (18) are suitable only for convex billiards. The boundary integral method on the other hand originates from the Green theorem, which provides an exact integral equation for a wave function in the interior of the billiard in terms of that on the boundary, and is applicable to arbitrary shapes. Furthermore, like the scaling method, the boundary integral approach has the huge advantage that the eigenvalue problem is reduced from a two-dimensional differential equation to a one-dimensional boundary integral. It was generalized to neutrino billiards in [26]. The starting point for the derivation of the boundary integral are the Dirac equation for y † , and the corresponding matrix equation for the free-space Green function where  denotes the 2×2 unit matrix and the Green function G r r , Here, we used that in two-dimensional space r r rr which is the scalar free-space Green function with H x  figure 1 for the definition of these quantities) for the free-space Green function Multiplying equation (33) from the left with r y ( ) † and equation (32) from the right with G r r , 0 ¢ ( ), subtracting the latter from the former and performing the r integration over the billiard domain Ω leads to The first equality results from the application of Gauss's theorem to transform the two-dimensional integral over Ω into a one-dimensional line integral over s along the domain boundary ¶W. Note that, if r¢ is chosen at a corner of the boundary, the factor 1 2 needs to be replaced by 2 in q p with θ in denoting the inner angle of the corner [38,39]. Using the boundary condition equation (20) we obtain for each wave function component a separate boundary integral equation, ). Furthermore, χ=1 for r¢ Î W ¶W ⧹ and χ=1/2 for r¢ Î ¶W and the arc length s and the parameter f are related through s w In order to attain a non-singular integrand, we apply in equation (38) the boundary condition equation (20) to (37) and obtain The wave functions inside the billiard are calculated by inserting the solutions of equations (40) or (42) into (37) and (38). Note, that f and f′ are not necessarily chosen as in equation (5), but generally denote the parameter which defines the boundary of the billiard. In appendix A we compare the boundary integral equation (37) to the corresponding one for a quantum billiard with Robin boundary conditions in order to illustrate their commonalities.
Generally, the boundary integral equation equation (40) needs to be solved numerically in order to determine the eigenvalues and eigenfunctions of a neutrino billiard. Only for a few exceptions, the eigenstates can be found either directly from the plane-wave expansion equation (18) or by solving the boundary integral equation analytically. Yet, when solving the boundary integral equation numerically one might face the problem of missing levels due to the presence of nearly degenerate eigenvalues. Such situations demand an expanded boundary integral method of the same type as the one for quantum billiards [27]. Thus, even though we can obtain exact analytical results either directly from an expansion in terms of plane waves or from the boundary integral equation, the computation of all eigenvalues of the circular neutrino billiard by numerically solving the latter is a non-trivial task. Therefore we chose it as an example to illustrate the efficiency of this method in neutrino billiards. In appendix B we outline in detail the analytic procedure for solving the boundary integral equation for the circle billiard. For systems with nearly degenerate eigenvalues and also for billiards with holes [27] or inner corners [40], we propose an expanded boundary integral method which is based on [27]. First, we briefly review in section 5 the procedure for solving the boundary integral equation (41) numerically and then outline in section 6 the expanded boundary integral method.

Boundary integral method for neutrino billiards
In order to solve the boundary integral equations (40) with (41) numerically, the boundary parameter f is dicretized and the integral is approximated by a sum, turning the boundary integral equation into a matrix equation . 44 Here, the matrix elements Q k , ; j i f f ( )are obtained from equation (41) by replacing f′ by f j and f by f i . One possible choice is to partition the boundary parameter f into equal-size pieces i yielding a Riemann sum. Another one is to use increments generated by the Gauss-Legendre algorithm. For a given N this discretization generally yields the integral and thus the eigenvalues with a considerably higher accuracy than the former one. Note that, if f does not coincide with the arc length s, partioning of 0, 2 f p Î [ ) into pieces Δf i corresponds to a discretization of the billiard perimeter into segments s w , and thus to a higher density of supporting points in regions where the boundary curve is stronger bent; see figure 2. Because the integrand in equation (40) does not contain any singularities, the numerical approximation equation (44) approaches the boundary integral with increasing number of partitionings N. However, the larger N, the higher will the numerical effort and computation time be. For a given k-range an optimal value for N is obtained in terms of multiples c λ of the number of de Broglie wave lengths k 2 l = p fitting into the perimeter  of yields the matrix equation which has solutions at discrete values of k n corresponding to the eigenvalues of the neutrino billiard. We would like to stress that, in distinction to, e.g. the boundary integral problem for quantum billiards with Robin boundary conditions (see equation (A.6)) where one has to deal with the singularity exhibited by the H k 0 r ( ) Hankel function at ρ=0, this is not the case for neutrino billiards because there the corresponding matrix elements Q k , ; i i f f ( )vanishes as outlined in section 3. Even though this singularity is logarithmic in ρ and thus lifted when performing the integral over f, it affects the accuracy of the eigenvalues.
The eigenvalues k n correspond to singular values of k ( ), that is, to zeros of the determinant of , k det 0 . Since equation (47) corresponds to an approximation of the actual boundary integral equation, k ( ) will in general not be exactly singular, i.e. its determinant will generally be non-vanishing. Thus, in order to find the eigenvalues, one may either use the singular value decomposition method for complex matrices and then sort the eigenvalues by size of their absolute values and plot the smallest one or else k det  | [ ( )]|versus k and seek for local minima. In order to save computation time this can be done by successively refining the discretization. However, this procedure might not work if two or more eigenvalues are too close to each other, as it may occur, e.g. for eigenvalue spectra exhibiting Poisson statistics, or if a dip is too narrow. Note, that the problem of spurious eigenvalues possibly occurring when solving the Dirichlet boundary problem for a quantum billiard using a double layer equation with no singularities of the integrand instead of the original boundary integral equation [38], does not arise for neutrino billiards containing no holes since there singularities are absent. We will use the Monza billiard [27] to illustrate the efficiency of the expanded boundary integral method in neutrino billiards of which the shape comprises holes.

Expanded boundary integral method for neutrino billiards
In order to identify missing levels we use an expanded boundary integral method which corresponds to an extension of the one proposed for quantum billiards in [27] to neutrino billiards. The basic idea of this method is similar to that of Vergini and Saraceno [37]. They expanded the boundary norm, which vanishes at the eigenvalues of the corresponding quantum billiard, in a Taylor series around a reference value k n 0 ( ) to obtain a generalized eigenvalue problem of which the solutions provide the eigenvalues in a narrow interval around k n 0 ( ) . In our case the boundary norm does not vanish at an eigenvalue k n of the neutrino billiard, yet k ( ) in equation (47) has a singular value [27]. Accordingly, we choose a reference value k n  To find appropriate initial values for k n 0 ( ) and for the k-interval, we first solve equation (47) for a larger k range and then identify regions of missing levels by, e.g. analyzing the fluctuating part of the integrated level density [38]. In order to determine the values k n of the missing levels we solve a generalized eigenvalue problem similar to the one proposed in [27]. There, a perturbation parameter ò was introduced and then k n , k ( ) and u k ( ) were expanded in terms of powers of it where we suppressed the argument k n 0 ( ) in the  (·) and u (·) -terms Inserting these expansions into equation ( In a first approximation we solve this generalized eigenvalue problem up to o(ò 2 ), The o(ò 2 ) contribution to the eigenvector u is obtained with this result and equations (54) from (58) as This procedure yields the corrections k n but not ò itself. Yet, for the consistency of this approach it is sufficient to know these terms to ensure that the required accuracy is achieved which is essentially determined by their values and the size of the next-leading order contribution k n 3 3  ( ) . For their determination we need to differentiate the matrix elements of k ( ), that is, according to equation (46) the matrix Q defined in equation (41) Note that also the derivatives of Q k , ; j i f f ( )do not contain singularities, because each differentiation of the Hankel functions with respect to k is accompanied by an additional factor ρ ij which vanishes df µ| |with δf→ 0, and thus cancels the singularities of the associated Hankel function. Furthermore, the wave functions were determined by inserting the eigenvalues (61) into equation (47) and proceeding as in the standard method.
Missing along the k-axis and using the above described procedure for each reference value k n 0 ( ) . All the resulting sequences need to be merged into one. Generally this is possible by simply cutting the individual sequences at a value k k c u = ( ) where they start to overlap. In cases where an eigenvalue obviously is obtained twice, the one with the higher accuracy, that is, producing a value of k det n  [ ( )]closest to zero, is chosen. Yet, due to the numerical errors, it can happen that either a genuine eigenvalue is thereby removed or doubled [27].
In order to avoid that k c if there is some eigenvalue closer than 3 to k c u ( ) . In such a case we define a new value for the border between the sequences, k 3 is satisfied.

Test of the expanded boundary integral method in the circular neutrino billiard
For a test the expanded boundary integral method and comparison with the standard one we determined the first 1000 eigenvalues with each method independently and compared them with the exact eigenvalues of the circular neutrino billiard obtained from equation (B.7) which indeed exhibits many near degeneracies. Here we chose 0.001  = for the expanded, and c λ =8 for the standard boundary integral method. This yielded in both cases the eigenvalues with a similar accuracy, yet for the latter one several eigenvalues could not be detected because they were too close to a neighboring one. One example of two nearly degenerate eigenvalues is given in table 1. Figure 4 shows the real part, the imaginary part and the phase of the components ψ 1 and ψ 2 of the associated pair of wave functions. The distance of the eigenvalues equals k 5 10 6 D »´-| | corresponding to 1/ with A=π denoting the area. Accordingly, the determination of all eigenvalues with the standard method of finding local minima of the determinant k det  | ( )| requires in the vicinity of degeneracies a discretization c λ 25. Thus, the numerical efforts related with the identification of regions where eigenvalues are missing and then finding them is large in the standard method whereas we did not encounter the problem of missing levels in the expanded one, that is, it is more reliable and consideraly less timeconsuming in such situations. We made similar experiences with rectangular, elliptic and triangular neutrino billiards which also exhibit nearly degenerate eigenvalues.

Application of the expanded boundary integral method to the Monza neutrino billiard
In the above mentioned examples the eigenvalues can be obtained directly based on the plane-wave expansions equation (18). In the present section we will consider a case where this is no longer possible, namely the Monza billiard [27] which for the following reasons provides a stringent test of their efficiencies. Firstly, since the Monza billiard contains a hole, one has to face with spurious eigenvalues. Secondly, the eigenvalue spectrum comprises closely-lying eigenvalues. This is attributed to the unidirectionality of its classical dynamics, implying that the Monza billiard exhibits the particular property that the motion of a particle launched into it with a certain rotational direction will follow this direction forever [42][43][44]. Consequently, the classical phase space is divided into two disjoint regions. They are separated by a family of marginally stable bouncing-ball orbits which are reflected back and forth perpendicularly at the inner and outer boundary parts and thus do not perform rotational motion, that is, they form the invariant separatrix between clockwise and counterclockwise motion. The separation of the phase space, in turn, gives rise to an extraordinary structure in the eigenvalue spectrum of the non-relativistic quantum and the neutrino Monza billiard, which can be split into doublets and singlets. The singlets correspond to the bouncing-ball modes which are of measure zero in classical phase space [27,[43][44][45]. The appearance of doublets has its origin in the separation of classical phase space into two parts. Yet, for the quantum case the eigenvalues are not degenerate, as would be expected if switching between clockwise and counterclockwise motion were strictly forbidden also in this limit. They are, actually, split by a distance which is small compared to the mean spacing of the doublets [27]. This splitting is attributed to dynamical tunneling [27,43,44] between the two regions of phase space through the barrier region associated with the bouncing-ball orbits. Note, that the splittings do not decrease exponentially with increasing eigenvalue number, as might be expected in the semiclassical limit, but rather algebraically. Accordingly, almost all eigenvalues are nearlydegenerate in the quantum Monza billiard. It was shown in [27] that the classical dynamics in each invariant half of phase space is fully chaotic and ergodic and the spectral properties of the doublets of the Monza quantum billiard were found to exhibit GUE rather than GOE behavior [27]. This reminds on the spectral structure  present in quantum billiards with threefold symmetries [46]. Yet, in contrast to the latter the Monza billiard does not exhibit any geometric symmetries. These exceptional features of the classical and quantum system immediately brought us to the question to what extent they are also present in the Monza neutrino billiard. We chose for the Monza billiard the same shape as in [27]. It is composed of four parts with straight inner and outer boundary segments at the same distance r d , where two parts have length a, respectively b, and lie opposite to each other with respect to the hole. They are joined by three ring shaped parts of width r d defined by parameters α, q, and r. We chose for the parameters the same values as in [27], q=1/12, a=1/2, b=1/3, r=1/3, α=1 and computed its levels up to k;65 using the standard and the expanded boundary integral method. The discretization parameter was chosen to be c λ =12 and the parameter ˜setting the spacing between the reference value k n u ( ) between individual runs equation (72) equaled 0.001  = . The resulting eigenvalue sequences contain for both methods ;100 spurious eigenvalues. This is attributed to the fact that, due to the hole, the boundary integral equation involves a combination of boundary integrals along the inner and outer boundary which might result in a zero of A k det ( ) within the numerical error, even though the equation is not exactly fulfilled along the inner boundary like in the example shown in figure 5. Checking whether the wave function components vanish inside the hole does not always work, because of the unavoidable numerical error. Therefore, in order to test whether the eigenvalues are genuine we checked the orthogonality of the corresponding eigenfunctions n y with respect to those corresponding to the upper neighboring eigenvalues, C m n n n m , , 1, 2, , nm n m max y y = = + + ¼ + | ⟨ ⟩ | with m max ; 10. Again, because of the numerical accuracy we cannot expect that C nm vanishes exactly for n m ¹ in the case of orthogonality. We found out that by setting the limit to C nm  0.051 for orthogonality most of the spurious eigenvalues could be eliminated. Finally, we used the Weyl formula equivalent for neutrino billiards in order to identify non-genuine eigenvalues or regions of missing ones. It turned out that the fluctuating part of the integrated spectral density, N k N k N k n n n fluc smooth = -( ) ( ) ( ), n=1, L, 2170, which is depicted in figure 6 exhibits fast and clearly visible slow oscillations which hamper this search and are due to the bouncing-ball orbits mentioned above. Furthermore, the right panel shows the length spectrum, i.e. the absolute value of the Fourier transform of the fluctuating part of the spectral density which exhibits peaks at the lengths of the periodic orbits. Here, below l;10 the bouncing-ball orbits clearly dominate.
Their contribution can be derived explicitly [45] based on the trace formula for a rectangular neutrino billiard with r d denoting the distance between the inner and outer boundary, bbo  the total length of the boundary where bouncing-ball orbits occur. It, and also the length spectrum deduced from it are shown as red dashed lines in figure 6, thus illustrating that it indeed describes the slow oscillations exhibited by the fluctuating part of the integrated spectral density. The summand of the trace formula contains an extra phase (−1) m as compared to the one applicable to the bouncing-ball orbits in the corresponding quantum billiard which is due to the differing boundary conditions. After subtracting this contribution from N fluc (k), missing and spurious levels are clearly visible as jumps in it as demonstrated in figure 7.
We did not encounter any missing levels in the expanded boundary integral method whereas we had to increase the numerical effort, associated with locating the regions where an eigenvalue was missing and their identification which required a discretization of c λ 25, considerably in order to achieve a complete sequence of 2170 eigenvalues below k≈65 in the standard one.
Like in the stadium billiard, bouncing-ball orbits are of measure zero in the phase space of the Monza billiard. Yet, they have a clearly visible effect on the spectral properties of the Monza neutrino billiard, as illustrated in figure 8.
Before analyzing the fluctuation properties in the level sequence of a neutrino billiard the eigenvalues need to be unfolded such that the spectral density is uniform. This is generally achieved by replacing the k n by the smooth part of the integrated spectral density which is given by [26] N k n Ak smooth 4 with A denoting the area of the Monza billiard. Shown are the distribution of the spacings between adjacent eigenvalues P(s), the corresponding cumulative distribution I(s), the variance Σ 2 (L) of the number of eigenvalues in an interval L and the Dyson-Mehta statistic Δ 3 (L) [47], which gives the local average least-square deviation of the integrated spectral density of the unfolded eigenvalues from a straight line over an interval of length L. The curves are compared to the corresponding ones for Poisson statistics (dashed lines), GOE (solid lines), GUE (dashed-dotted lines) and the case of random matrices consisting of two diagonal blocks of the same dimension containing GUE matrices (red lines). The left panels show the spectral properties before eliminating the bouncing-ball orbits, the right panels after their extraction. Like in the stadium billiard this is done by replacing the eigenvalues k n by k N k N k n n n smooth bbo = + ( ) ( )instead of just the smooth part. Note that unfolding refers to the slowly varying part of the integrated spectral density which in the present example includes the slow oscillations caused by the bouncing-ball orbits. Then the spectral properties agree well with the 2GUE curves. This implies, that in distinction to the Monza quantum billiard, the two doublet partners corresponding to clockwise and counterclockwise motion in the Monza billiard, are well separated, i.e. are well split, indicating that dynamical  tunneling is much stronger in the neutrino than in the quantum Monza billiard. Yet, deviations are observed even after extraction especially in P(s) for small spacings. These may be attributed to remnants of contributions from the singlets and a weak coupling between the doublet partners which may be accounted for in a random matrix ensemble of the form used in [44] in which the two GUE blocks are weakly coupled. We conclude this section with two examples for the wave functions, one showing a bouncing-ball mode, the other one three chaotic ones, which correspond to three successive eigenvalues, see figures 9 and 10. They were computed by inserting the eigenvalues, which were determined with the expanded boundary integral method, into the original boundary integral equations (37) and (38). In order to check their accuracy, i.e. to verify that the boundary condition is fulfilled, we first computed the boundary functions by choosing r¢ Î ¶W which indeed was the case. Actually, the absolute value of the corresponding smallest eigenvalue of equation (47) was less than 10 −4 . Then, we computed the wave functions inside the billiard by choosing r¢ Î W ¶W ⧹ . The distance between the eigenvalues k 399 and k 400 is by a factor of 30 smaller than that to the adjacent ones, so they form a doublet, whereas k 398 is framed by this doublet and another one, thus it is a singlet state and indeed looks like a slightly distorted bouncing-ball mode.

Conclusions
We introduced an expanded boundary integral method which is especially suited when dealing with neutrino billiards of which the eigenvalue spectrum exhibits near degeneracies. In order to test the efficiency of the standard and expanded boundary integral method we applied both to the circle neutrino billiard and also to neutrino billiards, of which the eigenvalues can be determined directly from the plane wave expansion equation (18). It turned out that the numerical effort associated with the computation of complete sequences of eigenvalues is considerably more time-consuming when using the standard boundary integral method, because  the discretization of the billiard boundary needs to be notedly refined and several search iterations might be needed before a missing eigenvalue is detected. Finally we applied the expanded boundary integral method to the Monza neutrino billiard which is interesting per se because it has a hole and above all belongs to the family of unidirectional billiards. We found out, that like in the corresponding quantum billiard its eigenvalue spectrum can be separated into singlets and doublets of close lying ones. Yet, in distinction to the latter, the spectral properties coincide with those of a superposition of two independent GUEs, thus indicating that dynamical tunneling is stronger in the neutrino billiard than in the corresponding non-relativistic quantum billiard.
The Schrödinger equation of a quantum billiard with domain Ω is given by that of a free particle,  , k 399 =26.505 and k 400 =26.507. The orthogonality between the state 398 and 399, respectively 400 equaled 0.002, that between the latter two equaled 0.006 5.