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 -Δ+ V. When the potential V₁ is separable, the operator -Δ+ V₁ is inverted or exponentiated at cost O\(N^1+1/d\) in d dimensions via per-axis eigendecomposition. On a single NVIDIA A100 GPU, this costs less than one second for 10⁹ degrees of freedom in 3D. For non-separable potentials V = V₁ + V₂, the same solver provides a preconditioner \(-Δ+ V₁\)^-1 for the preconditioned conjugate gradient (PCG) method and a propagator for operator-splitting time integrators. For bounded V₂, 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₁ 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.