Abstract
Determining quantum excited states is crucial across physics and chemistry but presents significant challenges for variational methods, primarily due to the need to enforce orthogonality to lower-energy states, often requiring state-specific optimization, penalty terms, or specialized ansatz constructions. We introduce a novel variational principle that overcomes these limitations, enabling the \textit{simultaneous} determination of multiple low-energy excited states. The principle is based on minimizing the trace of the inverse overlap matrix multiplied by the Hamiltonian matrix, Tr(S⁻¹H), constructed from a set of \textit{non-orthogonal} variational states {|ψ_i}. Here, H_ij = ψ_i | H | ψ_j and S_ij = ψ_i | ψ_j are the elements of the Hamiltonian and overlap matrices, respectively. This approach variationally optimizes the entire low-energy subspace spanned by {|ψ_i} without explicit orthogonality constraints or penalty functions. We demonstrate the power and generality of this method across diverse physical systems and variational ansatzes: calculating the low-energy spectrum of 1D Heisenberg spin chains using matrix product states, finding vibrational spectrum of Morse potential using quantics tensor trains for real-space wavefunctions, and determining excited states for 2D fermionic Hubbard model with variational quantum circuits. In all applications, the method accurately and simultaneously obtains multiple lowest-lying energy levels and their corresponding states, showcasing its potential as a unified and flexible framework for calculating excited states on both classical and quantum computational platforms.