Fast computation of high frequency Dirichlet eigenmodes via the spectral flow of the interior Neumann-to-Dirichlet map

Abstract: We present a new algorithm for numerical computation of large eigenvalues and associated eigenfunctions of the Dirichlet Laplacian in a smooth, star-shaped domain in $\mathbb{R}^d$, $d\ge 2$. Conventional boundary-based methods require a root-search in eigenfrequency $k$, hence take $O(N^3)$ effort per eigenpair found, using dense linear algebra, where $N=O(k^{d-1})$ is the number of unknowns required to discretize the boundary. Our method is O(N) faster, achieved by linearizing with respect to $k$ the spectrum of a weighted interior Neumann-to-Dirichlet (NtD) operator for the Helmholtz equation. Approximations $\hat{k}_j$ to the square-roots $k_j$ of all O(N) eigenvalues lying in $[k - \epsilon, k]$, where $\epsilon=O(1)$, are found with $O(N^3)$ effort. We prove an error estimate $$ |\hat k_j - k_j| \leq C \Big(\frac{\epsilon^2}{k} + \epsilon^3 \Big), $$ with $C$ independent of $k$. We present a higher-order variant with eigenvalue error scaling empirically as $O(\epsilon^5)$ and eigenfunction error as $O(\epsilon^3)$, the former improving upon the 'scaling method' of Vergini--Saraceno. For planar domains ($d=2$), with an assumption of absence of spectral concentration, we also prove rigorous error bounds that are close to those numerically observed. For $d=2$ we compute robustly the spectrum of the NtD operator via potential theory, Nystr\"{o}m discretization, and the Cayley transform. At high frequencies (400 wavelengths across), with eigenfrequency relative error $10^{-10}$, we show that the method is $10^3$ times faster than standard ones based upon a root-search.