Multi-state quantum simulations via model-space quantum imaginary time evolution

The study addresses the challenge of reliably preparing not only ground states but also excited states on noisy intermediate-scale quantum (NISQ) devices. While variational quantum algorithms, especially VQE, have been widely used for ground and excited states, imaginary time evolution (ITE) based approaches such as QITE have shown promise for ground states without costly parameter optimization. However, standard QITE and QLanczos can struggle to retain excited-state components due to exponential decay with imaginary time and numerical issues from linear dependence in Krylov spaces, particularly for higher excited states separated by large gaps. Prior approaches like folded-spectrum QITE (FSQITE) target excited states but require prior knowledge of target energies and handling H², increasing complexity. The paper proposes model-space QITE (MSQITE) to evolve an orthogonal subspace of states, maintaining orthogonality and enabling simultaneous multi-state simulation, with an extension of QLanczos (MS-QLanczos) to accelerate convergence. It also tackles spin contamination, which is especially problematic for excited states, by introducing a spin-shifted propagator.

Existing ITE-based quantum algorithms include VITE (using McLachlan's variational principle with a fixed ansatz) as an optimizer for variational and transcorrelated methods; probabilistic ITE (PITE), which probabilistically implements non-unitary evolution via an ancilla but suffers from exponential decay in success probability and requires SVD; and QITE, which approximates short imaginary-time non-unitary evolution with a unitary determined by linear equations and has been extended and demonstrated experimentally. For excited states, QLanczos constructs a Krylov subspace from QITE-propagated states, but excited-state amplitudes decay exponentially with imaginary time and are often lost to numerical noise from linear dependence. Folded-spectrum approaches (FSQITE) target energies by evolving with exp(−β²(H−ω)²), which can recover excited states and be accelerated by QLanczos but require estimating ω and handling H². In variational excited-state methods, numerous state-specific and state-averaged VQE variants exist; however, recent reports indicate state-averaged schemes using a single unitary for multiple orthonormal states can incur large errors for excited states, foreshadowing limitations for analogous state-averaged MSQITE formulations.

MSQITE starts from an orthogonal model space of n_states zeroth-order configurations {|Φ_l⟩} spanning targeted ground/excited states and evolves the entire subspace under imaginary time while maintaining orthonormality. A short-step evolution at step e with step size Δβ is approximated by a unitary generated by a Hermitian operator A (expanded in chosen Pauli strings) and followed by Löwdin symmetric orthonormalization to minimize state changes and preserve orthogonality. - Unitary approximation and orthonormalization: For each state l at step e, the exact short evolution e^{-Δβ(H−E_l)}|Φ_l^(e)⟩ is approximated by e^{-iΔβ A_l}|Φ_l^(e)⟩ with an energy shift E_l = ⟨Φ_l^(e)|H|Φ_l^(e)⟩. The transformation matrix d is chosen via Löwdin symmetric orthonormalization of the approximately evolved subspace to ensure orthonormality and minimal state change, with d → I as Δβ → 0. - Determining A: A is expanded as A=Σ_j a_j σ_j over selected Pauli strings (e.g., those corresponding to a chosen fermionic excitation ansatz). Two variants are proposed: 1) State-specific MSQITE: a separate A_l is determined for each state by minimizing F^l(a_l)=|| e^{-Δβ(H−E_l)}|Φ_l⟩ − e^{-iΔβ A_l}|Φ_l⟩ ||^2 to second order in Δβ, leading to linear equations M_l a_l + b_l = 0 with M_l and b_l obtained from expectation values on |Φ_l^(e)⟩. An additional term in b_l (proportional to elements of d) enforces orthogonality stability. 2) State-averaged MSQITE: a single A is shared by all states, determined from a pooled linear system M' a + b = 0. While circuit-simpler, this has limited representability for simple operator sets. - Quantum circuits: For state-specific MSQITE, off-diagonal Hamiltonian matrix elements H_IJ among model-space states require Hadamard tests with controlled e^{-iθA}. The circuit is arranged so that neighboring controlled Pauli rotations share basis-change and parity-aggregation subcircuits, minimizing overhead. The off-diagonal circuits incur ~1.2–1.3× the CNOT count of diagonal energy circuits. - MS-QLanczos: Extends QLanczos to a model-space Krylov basis composed of normalized model-space states at multiple imaginary-time steps { |Φ_l^(e)⟩ }. Effective overlap S and Hamiltonian H matrices are constructed using accumulated transformation matrices D^(l) built from the stepwise orthonormalization matrices d, along with a global reference energy shift E_0 to stabilize norms. The generalized eigenvalue problem H c = ε S c yields approximate eigenstates/energies. Stabilization heuristics (excluding near-linearly dependent vectors) are used. - Spin contamination control: A spin-shifted propagator penalizes higher-spin components: evolve with e^{−β (H + λ (S^2 − s(s+1)))} (equivalently shifting by a positive λ so components with s' > s decay faster). Preparing initial states with m_s = s ensures components with s' < s are absent. A moderate λ suppresses spin contamination while controlling Trotter error. - Simulation details: Implemented in QUKET with OpenFermion, PySCF, Qulacs; STO-6G basis and HF orbitals with Jordan–Wigner mapping; Δβ typically 0.1 a.u. (QITE/MSQITE) and 0.05 a.u. (FSQITE). UCCGSD ansatz used for A (except simplified circuits in noise tests). Linear systems solved with SVD truncation threshold 1e−7 of max singular value. For noise simulations (Qiskit), two-qubit depolarizing error rates p2 ∈ {1e−2, 5e−3, 1e−3}, one-qubit p1 = 0.1 p2; zero-noise extrapolation (ZNE) with exponential fit applied to CNOTs. H₄ Hamiltonian tapered to 4 qubits with a reduced ansatz; Krylov stabilization includes backward selection and limiting basis size.

- MSQITE enables simultaneous evolution of multiple states in an orthonormal model space, maintaining orthogonality via Löwdin symmetric orthonormalization and allowing post-evolution diagonalization in the developed model space to obtain energies and states. - Compared to FSQITE and standard QLanczos for excited states, MSQITE (especially with MS-QLanczos) achieves faster and more stable convergence without requiring prior target energy estimates or handling H². - BeH₂ (equilibrium, STO-6G): With a 3-configuration model space (HF and two π_u excitations), MSQITE converges rapidly to X¹Σ_g^+, ¹Δ_g, and ²Σ_g^+; convergence is state-dependent (¹Δ_g around β ≈ 5 a.u.; ²Σ_g^+ several a.u. later) due to strong correlation. FSQITE/FS-QLanczos converge slower and cannot target the ground state when centered on excited-state energies (e.g., ω = E_2 = 15.2263 Hartree for ²Σ_g^+ in folded spectrum). - H₄ square (1 Å): Strongly correlated two-determinant ground state. QITE requires >10 a.u. imaginary time to reach 1 mHartree accuracy; QLanczos ~4 a.u. MSQITE with a 2-determinant model space reaches near-exact energies within a few a.u.; MS-QLanczos is even faster. - N₂ (stretched 1.6 Å): Without spin control, MSQITE excited states collapse into lower-energy higher-spin states (triplets/quintets) and exhibit spin contamination (e.g., 〈S²〉 ≈ 2 and energies trapped near the mean of exact triplet and singlet energies, −108.378 a.u. ≈ (−108.293193 −108.463729)/2). Introducing a spin-shift (λ = 0.5) maintains 〈S²〉 ≈ 0 and drives all targeted singlet states to convergence, eliminating spin contamination. - State-specific vs state-averaged MSQITE (noiseless): State-specific MSQITE with UCCGSD achieves high accuracy across systems, while state-averaged MSQITE deteriorates as system and model space size grow. Illustrative errors at convergence (Hartree): • H₂ X¹Σ_g^+: state-specific <1e−8; state-averaged 2e−7. • H₂ a³Σ_g^+: state-specific <1e−8; state-averaged 2e−7. • BeH₂ X¹Σ_g^+: state-specific <1e−8; state-averaged 2e−4. • BeH₂ ²Σ_g^+: state-specific 8e−8; state-averaged 4e−4. • BeH₂ ¹Δ_g: state-specific <1e−8; state-averaged 2e−4. • N₂ X¹Σ_g^+: state-specific 6e−5; state-averaged 5e−3. • N₂ ³Σ_g^−: state-specific 8e−5; state-averaged 5e−2. The state-averaged scheme may even stall by symmetry (⟨b'⟩ = 0) with limited operator sets; high-rank excitations (e.g., UCCGSDT/Q) are required for completeness in small-electron systems. - Circuit complexity: Off-diagonal MSQITE circuits (Hadamard tests with controlled unitaries) require only ~1.2–1.3× CNOTs versus diagonal energy circuits, reflecting efficient gate sharing. - Noise studies (H₄, depolarizing noise): Circuit depth grows linearly with β, leading to rapid error accumulation; with p2 = 1e−2, energies exceed HF already by β ≈ 0.2 a.u. Off-diagonal elements tend to zero under noise. Zero-noise extrapolation markedly improves results at small β, and mitigated MSQITE shows stability comparable to mitigated QITE. Errors in diagonal and off-diagonal effective Hamiltonian elements after ZNE are of similar magnitude, so overall MSQITE energy accuracy tracks QITE under mitigation. - Overall, MSQITE outperforms standard QITE/QLanczos and FSQITE/FS-QLanczos for excited-state preparation, and with a spin-shift yields spin-pure targets.

MSQITE addresses the loss of excited-state components during imaginary time evolution by evolving an entire orthonormal model space and preserving its orthogonality at each step. This ensures targeted excited-state character is retained and developed, enabling accurate extraction of excited states via a generalized eigenvalue problem in the propagated model space. The Löwdin orthonormalization minimizes state mixing between steps, upholding the validity of the short-time unitary approximation. Extending Lanczos to the model-space (MS-QLanczos) further accelerates convergence by leveraging a richer Krylov basis across states and time steps while controlling linear dependence. Spin contamination, a key impediment in excited-state simulations, is effectively mitigated by a spin-shift in the propagator that penalizes higher-spin components, projecting dynamics onto the desired spin sector without explicit projectors. This removes spin collapse in challenging cases like stretched N₂ and yields clean convergence to singlet states. Empirically, state-specific MSQITE using compact operator sets (e.g., UCCGSD) achieves high accuracy for ground and excited states across test systems, whereas state-averaged MSQITE often lacks sufficient representability unless high-rank excitations are included, limiting its practicality. In noisy settings, although circuit depth growth challenges both QITE and MSQITE, efficient circuit design and error mitigation (ZNE) enable MSQITE to preserve its advantages with noise impacts comparable to QITE. These results demonstrate that MSQITE provides a robust route to multi-state quantum simulations on near-term devices, particularly for strongly correlated problems and excited-state targets.

The work introduces MSQITE, a model-space formulation of quantum imaginary time evolution that maintains orthogonality of multiple states and enables simultaneous, stable preparation of ground and excited states. A Löwdin-based orthonormalization ensures minimal state change at each step, and MS-QLanczos accelerates convergence by exploiting a multi-state Krylov space. A spin-shifted propagator effectively eliminates spin contamination, enabling efficient targeting of specific spin sectors. State-specific MSQITE with modest operator sets (UCCGSD) delivers high accuracy, whereas state-averaged MSQITE requires substantially more complex operators to achieve comparable performance. Noise simulations show that, with zero-noise extrapolation, MSQITE retains advantages over QITE and exhibits similar noise sensitivity. Future directions include applying MSQITE beyond quantum chemistry (e.g., nuclear physics), integrating MSQITE with adaptive and variational algorithms to reduce circuit depth, improving stabilization strategies for MS-QLanczos, and developing symmetry-preserving, low-depth operator sets and error mitigation tailored to multi-state evolutions.

- Circuit depth grows with imaginary time, limiting practicality on NISQ hardware due to accumulated noise; mitigation is required for usable accuracy at modest β. - MS-QLanczos can suffer from linear dependence in enlarged multi-state Krylov spaces, requiring truncation/stabilization heuristics that may omit useful vectors. - State-averaged MSQITE exhibits limited representability with compact operator sets (e.g., up to doubles); achieving high accuracy often demands higher-rank excitations, increasing gate counts and measurements. - Performance depends on the quality and completeness of the initial model space; missing important configurations can slow convergence or require larger β. - Trotterization and short-time unitary approximation errors can accumulate; choice of Δβ and ansatz affects accuracy and stability. - Spin-shift requires selecting a penalty parameter λ; overly large λ can exacerbate Trotter errors, whereas too small λ may not fully suppress contamination. - Off-diagonal element estimation introduces additional circuits (though modestly larger than diagonal), adding to measurement overhead.