Numerical rank of the basis set

Examination of the matrices $F$ and $G$ shows that they both become singular [188] at almost exactly the point when the basis set size (matrix order) $N$ becomes sufficiently large to represent the solutions at the desired $k$ to sufficient accuracy. This makes (5.14) a singular generalized eigenproblem. Such singularity causes failure of standard generalized matrix eigensolvers (for instance dsygv in LAPACK [5]). This might be seen as a huge setback--however I will show that singularity seems to be an inevitable result of using a basis which obeys the Helmholtz equation (however, the issue of this singularity in the special case of the BIM basis set is unknown). Also I will discuss correct handling of the singular problem.

Why do the matrices become singular? The answer is that physically, wavefunctions $\psi{({\mathbf r})}$ which are arbitrary solutions of the Helmholtz equation at a single wavenumber $k$ are effectively limited to a finite number of degrees of freedom inside a given billiard. This can be seen intuitively by realising that two real plane waves (RPWs) much closer in angular separation than $(k{\mathsf{L}})^{-1}$ (where ${\mathsf{L}}$ is the typical system length) cannot be distinguished when only their values inside the billiard are considered. The overlap of the two RPWs in $\mathcal{D}$ becomes very similar to their individual norms in $\mathcal{D}$, i.e. they tend to the same function in $\mathcal{D}$. (In contrast, when considered over all space, the RPWs are always orthogonal). This would indicate that in $d=2$ a basis size of $N \sim k{\mathsf{L}}$ is sufficient to represent all solutions at $k$. This has been formalised in a more rigorous fashion [64] for convex billiards in $d=2$, as a semiclassical basis size $N_{sc} = kL/\pi$, where $L$ is the perimeter (length of $\Gamma$). This is simply the number of half-wavelengths along the perimeter. (I believe the generalization to higher $d$ would be that $N_{sc}$ is given by the Weyl law for a $d{-}1$ dimensional `volume' ${\mathsf{A}}$).

Not only is this basis size sufficient, but attempts to improve accuracy by increasing $N$ beyond this rapidly fail. Above $N_{sc}$ it becomes exponentially hard to introduce new (non-evanescent) functions into the basis. This can be illustrated by considering an angular-momentum ($l$) basis around some origin inside $\mathcal{D}$. The inside wavefunction is strongly affected by states with $\vert l\vert < kR$ (where $R$ is the maximum radius from the origin). For higher $\vert l\vert$ the billiard falls into the `classically-forbidden' region (the circle of radius $\approx \vert l\vert/k$), so these $l$-states have an effect on $\psi$ which dies exponentially with $\vert l\vert$. For these states, the only effect on $\psi$ is an exponentially-small evanescent wave.

Fig. 5.3 shows the eigenvalues $\lambda^F_n$ of the $F$ (tension) matrix. These $\lambda^F_n$ are simply the tensions produced by eigenvectors ${\mathbf x}^{(N)}$ of unit length in ${\mathbf x}$-space (coefficient space). The figure shows that up to $N_{sc}$ the eigenvalues are $O(1)$, but beyond this they plummet exponentially, finally settling at $\sim 10^{-16}$ which is the rounding error of the double-precision arithmetic used. By analogy with singular values this defines a numerical rank [188,80]: the number $N_r$ of eigenvalues above a certain cut-off $\epsilon$, which is chosen as roughly the working precision. Examining the figure, the $F$ matrix for $N=250$ has a rank of about 182, not much larger than $N_{sc} = 164$. The eigenvectors with $\lambda^F_n < \epsilon$ fall into the numerical null-space, and have a numerically irrelevant (tiny) effect on the value of ${\mathbf x}^{{\mbox{\tiny T}}}F {\mathbf x}$, and therefore on the value of $\psi$ anywhere within $\mathcal{D}$. An eigenvector heading towards the null space is shown in Fig. 5.4b; it has negligible value inside $\mathcal{D}$ and resembles a high angular-momentum state outside. What is surprising is that increasing $N$ further has little or no effect on the rank $N_r$. As Vergini has realised[194], this behavior has profound consequences for the representation of billiard eigenstates by RPWs. Namely, once the rank of the problem is reached (at a particular $k$) it is very hard to improved accuracy simply by `throwing' more real plane wave basis states at the problem (a method usually successful in other basis solution techniques). Most of the basis functions thus `thrown' will simply fall into the numerical null space without contributing to the accuracy of the $\psi$ found.

This self-truncation of the basis set size would appear quite bizarre to a numerical analyst trained on the likes of lattice and FEM approaches. In these latter approaches the basis set size $N$ can be increased without limit, and increased $N$ always results in higher accuracy of solution. However, because the basis sets of this chapter are already solutions of the wave equation, self-truncation occurs when the basis becomes (semiclassically) overcomplete. The exception to this is if suitable evanescent waves are added close to the boundary (see Appendix J). Therefore if highly accuracy is needed (or if the semiclassical limit is to be reached), EPWs must be added, as demonstrated by Vergini [194]. Discussion of the Helmholtz solution function space is continued in Section 6.1.3.