A spectral-element method for computing many eigenvalues and their asymptotics of the Schrödinger operator with Robin boundary condition

Technical Summary: A Spectral-Element Method for Computing Many Eigenvalues and Their Asymptotics of the Schrödinger Operator with Robin Boundary Condition

Problem Statement
The paper addresses the numerical computation of the eigenvalues and eigenfunctions of the Schrödinger operator $L_{SO} := -\Delta + V(x)$ subject to Robin boundary conditions on both simple and complex geometries. The specific eigenvalue problem seeks $\lambda \in \mathbb{R}$ and $u \not\equiv 0$ such that:
$$(-\Delta + V(x))u(x) = \lambda u(x), \quad x \in \Omega$$
$$\frac{\partial u}{\partial n} + \sigma u = 0, \quad x \in \partial\Omega$$
where $\Omega \subset \mathbb{R}^d$ is a bounded domain, $V(x)$ is a real-valued potential, and $\sigma \in \mathbb{R}$ is the Robin parameter.

The primary challenge identified is the need to compute thousands of high-index eigenvalues with high accuracy to study the asymptotic behavior of the Robin-to-Neumann (RtN) gaps, defined as $g_n(\sigma) = \lambda^\sigma_n - \lambda^0_n$. Standard low-order finite element or finite difference methods typically resolve only a limited portion of the spectrum for a fixed number of degrees of freedom (DOFs), making them unsuitable for investigating high-frequency asymptotics. While spectral methods offer high resolution, they are generally restricted to simple geometries (e.g., rectangles, disks).

Methodology
The authors propose a unified framework utilizing spectral methods for simple geometries and spectral-element methods (SEM) for complex/irregular geometries.

1. Variational Formulation: The problem is cast in a weak form using a conforming finite-dimensional subspace $X_N \subset H^1(\Omega)$. This leads to a generalized matrix eigenvalue problem $AU = \lambda_N MU$, where $A$ and $M$ are stiffness and mass matrices, respectively.
2. Spectral Methods (Simple Geometries): For domains like intervals, rectangles, disks, cubes, and balls, the authors employ global spectral basis functions. In 1D, Legendre polynomials are used to construct the basis. In higher dimensions, tensor-product bases are applied.
3. Spectral-Element Methods (Complex Geometries): For irregular domains, the domain $\Omega$ is decomposed into non-overlapping curvilinear elements. Within each element, high-order polynomial approximations (using Gauss-Lobatto-Legendre nodes and Lagrange bases) are employed. The method utilizes coordinate transformations from a reference square (or cube) to physical elements, incorporating Jacobian factors to handle geometric mappings.
4. Resolution Strategy: The methods are designed to achieve "spectral-type accuracy," where the number of eigenvalues computed with reasonable accuracy is proportional to the number of DOFs. This allows for the efficient computation of thousands of eigenvalues on a standard High-Performance Computing (HPC) node with 128 GB of memory.

Key Results
The authors validate their methods and apply them to study the RtN gaps and their cumulative averages $G_n(\sigma) = \frac{1}{n}\sum_{m=1}^n g_m(\sigma)$.

• Accuracy and Resolution:
  • 1D: With $N$ DOFs, the spectral method computes approximately $0.6N$ eigenvalues with relative errors below $10^{-8}$.
  • 2D: For simple geometries (spectral method), approximately $0.4N$ eigenvalues are resolved with high accuracy. For complex geometries (spectral-element method), approximately $0.25N$ eigenvalues are resolved.
  • Scale: Using $N \approx 25,000$ DOFs, the authors successfully computed over 5,000 reliable eigenvalues with relative errors $< 10^{-8}$ for various 2D complex domains (e.g., L-shaped domains, mushrooms, quarter stadiums).
• Asymptotic Analysis:
  • The computed eigenvalues confirm existing theoretical results for the Laplacian ($V \equiv 0$) in 2D, specifically that $\lim_{n \to \infty} G_n(\sigma) = \frac{2|\partial\Omega|}{|\Omega|}\sigma$.
  • The study extends these findings to the Schrödinger operator with non-trivial potentials $V(x)$.
  • Numerical results across 1D, 2D, and 3D domains with various potentials and boundary parameters $\sigma$ provide extensive numerical evidence supporting the derived asymptotic behaviors.

Proposed Conjecture
Based on extensive numerical evidence and observations, the authors speculate and formulate a unified conjecture regarding the RtN gaps for the Schrödinger operator with Robin boundary conditions on a connected bounded domain $\Omega \subset \mathbb{R}^d$ with piecewise smooth boundary:

1. Cumulative Averages: For any $\sigma \in \mathbb{R}$,
   $$\lim_{n \to \infty} G_n(\sigma) = \frac{2|\partial\Omega|}{|\Omega|}\sigma$$
   and for large $n$,
   $$G_n(\sigma) = \frac{2|\partial\Omega|}{|\Omega|}\sigma + O(n^{-1/d}).$$
2. 1D Specifics: When $d=1$ and $\Omega=(a,b)$, the individual gaps satisfy:
   $$\lim_{n \to \infty} g_n(\sigma) = \frac{4}{b-a}\sigma$$
   and
   $$\lambda^\sigma_n = \lambda^0_n + \frac{4}{b-a}\sigma + O(n^{-2}).$$

Significance and Claims
The paper claims to fill a gap in the literature where efficient, accurate numerical methods for computing many eigenvalues of the Schrödinger operator with Robin boundary conditions were previously limited. By combining the geometric flexibility of finite elements with the high-order accuracy of spectral methods, the proposed approach enables the resolution of thousands of eigenvalues, a task unattainable with standard low-order discretizations.

The primary contribution is not merely the algorithmic development but the systematic numerical investigation of the RtN gaps. This large-scale computation provides strong empirical support for a unified conjecture on the asymptotic behavior of these gaps, linking algebraic spectral quantities (cumulative averages) directly to geometric properties of the domain (area and boundary length) and the Robin parameter, even in the presence of non-trivial potentials. The work confirms theoretical results for the Laplacian and extends the understanding of these asymptotics to the more general Schrödinger operator.