Explicit Two-Sided Eigenvalue Bounds for Schrödinger Operators with Singular Potentials via Finite Element Method

Technical Summary: Explicit Two-Sided Eigenvalue Bounds for Schrödinger Operators with Singular Potentials via Finite Element Method

Problem Statement
The paper addresses the computation of explicit, computable two-sided bounds for the eigenvalues $\lambda_k$ of the Schrödinger operator $H = -\Delta + V$ on unbounded domains $\mathbb{R}^N$ ($N=2, 3$). The primary challenge lies in handling singular attractive potentials of Coulomb type, specifically $V(x) = -Z/|x|$ (modeling the hydrogen atom) and two-centered potentials $V(x) = -Z_1/|x-a_1| - Z_2/|x-a_2|$ (modeling the $H_2^+$ molecular ion).

Existing numerical methods typically provide upper bounds via the conforming Galerkin (Rayleigh–Ritz) method. Obtaining rigorous lower bounds is significantly more difficult, particularly for unbounded domains and singular potentials where the potential is unbounded below. Classical lower-bound techniques (e.g., Lehmann–Maehly–Goerisch, Kato–Temple) often require non-conforming methods or specific conditions that fail when the potential $V$ is negative and singular. Furthermore, standard finite element approaches for lower bounds often rely on shifting the potential by its $L^\infty$ norm of the negative part ($\gamma_h = \|c_h^-\|_{L^\infty}$). For Coulomb singularities, $\gamma_h$ diverges as the mesh is refined ($O(h^{-1})$), causing the lower bound to fail to converge.

Methodology
The proposed algorithm combines domain truncation with an extension of the Composite Enriched Crouzeix–Raviart (CECR) finite element method. The approach is structured in three main components:

1. Domain Truncation: The unbounded domain $\mathbb{R}^N$ is truncated to a bounded domain $D(R)$.

   • Upper Bounds: Obtained via Dirichlet boundary conditions ($u=0$ on $\partial D(R)$) using standard conforming $P_1$ Galerkin elements.
   • Lower Bounds: Obtained via Neumann boundary conditions ($\partial u/\partial n = 0$). A key lemma establishes that if the exterior potential minimum $\sigma_{ext}(D(R))$ exceeds the $k$-th eigenvalue, the Neumann eigenvalue on $D(R)$ provides a valid lower bound for the full-space eigenvalue. This condition is verified a posteriori using the computed upper bound.

2. CECR Framework for Sign-Changing Potentials:

   • The method employs a pair-space framework where functions are lifted to pairs $(u, \Pi_{0,h}u)$, separating the kinetic energy (gradient) from the reaction energy (potential).
   • To handle the singularity without divergent shifts, the authors introduce an elementwise reaction constant $\Gamma_h := \max_K c_{h|K}^- \cdot h_K^2 / \pi^2$. Unlike the global $L^\infty$ shift $\gamma_h$, $\Gamma_h$ vanishes ($O(h^2)$ or $O(h)$ depending on grading) as the mesh is refined, even for Coulomb singularities.
   • A fixed shift $\sigma > 0$ (independent of the mesh size $h$) is applied to the potential to ensure the coercivity of the shifted bilinear form. This avoids the need for mesh-dependent shifts that diverge.

3. Perturbation Analysis:

   • The continuous potential $V$ is approximated by its piecewise constant element average $c_h$.
   • A perturbation lemma bounds the error introduced by replacing $V$ with $c_h$. The perturbation ratio $\epsilon_h$ is shown to converge to zero ($O(h^2)$ on graded meshes) under specific mesh grading assumptions near the singularities.

Key Contributions
The paper claims the following specific contributions:

• First Explicit Convergent Lower Bounds: It presents the first numerical algorithm providing explicit, convergent two-sided eigenvalue bounds for Schrödinger operators with unbounded singular potentials on unbounded domains.
• Pair-Space Abstract Theorem: It develops a pair-space framework that lifts the CECR analysis to a two-component space, allowing for a self-contained max-min proof of the lower bound theorem for sign-indefinite potentials.
• Elementwise Reaction Constant ($\Gamma_h$): It introduces $\Gamma_h$ to replace the divergent $L^\infty$ shift. This allows the CECR gap to close at a rate of $O(h^2)$ on graded meshes, bypassing the need for uncomputable form-norm constants.
• Fixed-Shift Two-Sided Bound: It assembles a fully computable enclosure $L_k \le \lambda_k \le U_k$ that does not require $L^\infty$ control of the negative part of the potential.
• Application to $H_2^+$: It specializes the theory to the two-centered Coulomb potential, providing explicit formulas for admissible shifts and perturbation parameters in both 2D and 3D.

Numerical Results
The algorithm is validated on four test cases:

1. 2D Hydrogen Atom ($V = -1/|x|$): The method successfully brackets the exact ground state $\lambda_1 = -1$. The certified lower bound interval narrows significantly with mesh refinement, with the optimal-shift bounds achieving widths of $0.002$ at the finest level.
2. 2D $H_2^+$ Ion: Results for the two-centered potential show convergence, with the certified interval for the ground state narrowing from $[-4.43, -1.317]$ to $[-2.36, -1.317]$ as the mesh refines.
3. 3D Hydrogen Atom: On a ball domain, the method brackets the exact $\lambda_1 = -0.25$. The certified Neumann enclosure width converges at $O(h^2)$.
4. 3D $H_2^+$ Ion: For the molecular ion on a box domain, the method provides a certified interval $[-1.045, -0.544]$ (width $0.501$) for the ground state, which contains the reference value $\approx -0.5513$.

In all cases, the CECR lower bound gap $(\mu_{k,h}^\sigma - \sigma) - L_{k}^{\mu,\sigma}$ is observed to close at the predicted $O(h^2)$ rate on graded meshes. The paper notes that while the 3D gaps are wider than in 2D due to slower decay of the elementwise constant on coarse tetrahedral meshes, the convergence behavior is consistent with theory.

Significance and Claims
The paper asserts that this work constitutes the first finite element framework capable of providing explicit and convergent lower bounds for eigenvalues associated with unbounded singular potentials. The significance lies in overcoming the "coercivity gap" caused by attractive singular potentials, which previously caused standard lower-bound techniques to diverge.

The authors emphasize that the method is explicit: all constants (interpolation constants, perturbation bounds, shift parameters) are derived in closed form from the mesh and potential data. While the current implementation uses double-precision arithmetic, the authors explicitly state that conversion to fully verified (IEEE-754-safe, interval-arithmetic) enclosures is a separate verification step left for future work. The paper also notes that while the current bounds are projection-based, they serve as a rigorous starting point ($\rho$) for sharper post-processing techniques like Lehmann–Goerisch bounds, which could further tighten the intervals.

The work does not claim to solve the multi-electron problem directly but identifies the extension to magnetic Hamiltonians and multi-electron systems with antisymmetry constraints as natural future directions.