A Simple GPU-Accelerated Solver for the Schrödinger Operator with Applications to Ground States and Hamiltonian Simulation

We extend the tensor-product direct solver from the Laplacian to the Schrödinger operator $-\mathrm{\Delta }+V$. When the potential $V_1$ is separable, the operator $-\mathrm{\Delta }+V_1$ is inverted or exponentiated at cost $O(N^{1+1\mathrm{/}d})$ in $d$ dimensions via per-axis eigendecomposition. On a single NVIDIA A100 GPU, this costs less than one second for ${10}^9$ degrees of freedom in 3D. For non-separable potentials $V=V_1+V_2$, the same solver provides a preconditioner $(-\mathrm{\Delta }+V_1{)}^{-1}$ for the preconditioned conjugate gradient (PCG) method and a propagator for operator-splitting time integrators. For bounded $V_2$, we prove that the preconditioned operator has a bounded condition number and a clustered spectrum with at most finitely many outlier eigenvalues, independently of the mesh size, and also independently of the domain size when $V_1$ is a confining potential. This explains the mesh- and domain-independent PCG iteration counts observed in practice. We apply this method to ground state computation via inverse iteration for linear problems and via the $a_u$ gradient flow for Gross–Pitaevskii energy in 3D, and also Hamiltonian simulation via the approximated qHOP and Magnus-2 splitting methods from 3D to 9D on a single NVIDIA GH200 GPU.