Quantum simulation of excited states from parallel contracted quantum eigensolvers

Computing excited-state properties of molecules and solids is considered one of the most important near-term applications of quantum computers. While many of the current excited-state quantum algorithms differ in circuit architecture, specific exploitation of quantum advantage, or result quality, one common feature is their rooting in the Schrödinger equation. However, through contracting (or projecting) the eigenvalue equation, more efficient strategies can be designed for near-term quantum devices. Here we demonstrate that when combined with the Rayleigh–Ritz variational principle for mixed quantum states, the ground-state contracted quantum eigensolver (CQE) can be generalized to compute any number of quantum eigenstates simultaneously. We introduce two excited-state (anti-Hermitian) CQEs that perform the excited-state calculation while inheriting many of the remarkable features of the original ground-state version of the algorithm, such as its scalability. To showcase our approach, we study several model and chemical Hamiltonians and investigate the performance of different implementations.


I. INTRODUCTION
Calculating physical properties of excited-state processes of quantum many-body systems is one of the most promising applications of near-term quantum computing [1][2][3].Quantum devices are well suited to deal with many of the distinctive features of excited states such as their strong multiconfigurational character or the presence of conical intersections [4,5].So far, several quantum algorithms have been developed to approximate eigenstates of many-body Hamiltonians, including quantum phase estimation (QPE) [6,7] and the variational quantum eigensolver (VQE) [8,9].VQE has also inspired several related approaches for excited states: The two dominant variants rely on either targeting specific states through adding nonorthogonal penalties to the Hamiltonian [10][11][12][13][14] or by building subspaces while ensuring orthogonality of the lowest-lying eigenstates [15,16].Yet, QPE requires circuit depths beyond what is currently achievable, and VQE relies on high-dimensional classical optimization, which has computational costs that scale rapidly with the system size [17].
Quantum algorithms like QPE and VQE are designed to solve the Schrödinger equation (SE).However, more efficient quantum simulations can be performed if, instead of the standard SE, its contraction (or projection) is solved directly on a quantum computer [18].When solving the corresponding contracted Schrödinger equation (CSE) the prepared wave function ansatz only requires two-body terms, regardless of the number of electrons or orbitals, ensuring the scalability of the algorithm [19].While initially designed to explore ground states of molecular systems [19], quantum eigensolvers based on the CSE have been recently extended to excited states by using the variance of the energy as the cost function [20] or by deflating the CSE to ensure the orthogonality of the eigenstates [21].However, these methods compute the eigenstates individually and therefore the circuit must be run for each desired excited state.
The goal of this work is to demonstrate that when combined with the Rayleigh-Ritz variational principle for mixed quantum states, the CSE can be straightforwardly generalized for the simultaneous (or parallel) calculation of a bundle of lowest eigenstates.Our main result is a novel excited-state quantum algorithm that employs the main features of the ground-state contracted quantum eigensolver (CQE), thus retaining its favorable scaling.
Here we focus on the anti-Hermitian portion of the CSE which has been shown to render accurate approximations for ground-state calculations [22], but our results can be generalized to include its Hermitian part.In the same way, we focus on fermionic systems but our derivations equally hold for bosons.
The remainder of this paper is organized as follows: For completeness, we first introduce both the CSE and the Rayleigh-Ritz variational principle for ensembles, on which our algorithm is based.Next, we generalize the basic equations of the ground-state CQE to excited states and discuss the resulting quantum algorithm.We then present our contracted quantum eigensolvers, discuss different methods of implementing them on a quantum computer, and perform several numerical experiments.The paper ends with some conclusions and a discussion about potential future research directions.

II. THEORY
After we review the CSE and the the Rayleigh-Ritz variational principle for mixed states in sections II A and II B, we derive an anti-Hermitian CSE (ACSE) for arXiv:2311.05058v1[quant-ph] 8 Nov 2023 mixed states in section II C and a quantum algorithm based on this mixed-state ACSE in section II D, which can solve for multiple excited states simultaneously.

A. Contracted Schrödinger equation
The SE of an electronic system governed by a Hamiltonian Ĥ reads: The two-body operator Γpq st ≡ f † p f † q ft fs , where f † p / fp are fermionic creation/annihilation operators, followed by the vector ⟨ψ ν |, can be applied on the left of the SE in Eq. ( 1) to obtain the CSE: Both the CSE in Eq. ( 2) and the SE in Eq. ( 1) have an equivalent set of pure-state solutions [23][24][25]: while the SE clearly implies the CSE, the opposite direction is provable by showing that (2) implies the eigenstate condition of zero variance (i.e., ⟨ψ ν |( Ĥ − E ν ) 2 |ψ ν ⟩ = 0) which in turn implies the SE.Notice that Eq. ( 2) can be written as the sum of two terms (a commutator and anti-commutator) [26,27]: It is well-known that solving only the anti-Hermitian portion of this equation, i.e., gives accurate results both for ground bosonic [28] and ground and excited electronic [20,21] states.Moreover, since the Eq. ( 4) can be interpreted as the residual of a certain cost function, this anti-Hermitian CSE (ACSE) immediately suggests the type of ansatz that can be used to guess the form of the eigenstate |ψ ν ⟩ (see below).

B. Variational principle for ensembles
The Rayleigh-Ritz variational principle is a powerful tool routinely used to study eigenstates of quantum many-body systems [29].Its generalization to mixed quantum states establishes an upper bound for the weighted ensemble energy of the K lowest eigenstates of a Hamiltonian, Ĥ [30]: where ρ(w) = K−1 ν=0 w ν |ϕ ν ⟩⟨ϕ ν | is a density matrix with a positive, decreasingly ordered spectrum, conveniently defined as w = (w 0 , w 1 , . ..) with w ν ≥ w ν+1 ≥ 0. The vectors {|ϕ ν ⟩} can be any set of K orthogonal states.Here E ν ≤ E ν+1 are the exact eigenenergies of the system, arranged in increasing order.The ensemble variational principle in Eq. ( 5) offers a unified approach to variational methods in quantum mechanics: the problem of the ground state is, in fact, just a particular case, corresponding to w = (1, 0, 0, . ..).This variational approach to quantum excitations is currently playing a pivotal role in the extension of ground-state functional theories [31][32][33][34] and hybrid quantum-classical methods [15,16,35] to excited states.
We note in passing that Eq. ( 5) can be written in a state-specific form by employing the purified state [36]: The states |a ν ⟩ are auxiliary orthonormal (ancilla) states added to perform the purification.The only condition is their orthornormality, ⟨a ν |a µ ⟩ = δ νµ .Then, the lower bound of the energy expectation value of the ensemble energy can be written as ⟨ρ(w . .) and I being the identity matrix acting on the auxiliary space (we will skip the writing of I when the notation is obvious).

C. The ACSE for excited states
The generalization of the CSE to ensembles of eigenstates is straightforward.Indeed, since Eq. ( 2) is valid for all the eigenstates of the Hamiltonian Ĥ, one can use it to write a weighted sum for the first K eigenstates: From this equation, the corresponding ACSE for an ensemble of K eigenstates follows: This result suggests a variational implementation of the ACSE for excited states.Consider first a variational ansatz for a set of K orthogonal wave functions, iteratively constructed from unitary two-body exponential transformations: where is an anti-Hermitian two-electron operator and η is a real positive number (whose role will be clear later).The ensemble energy at the (n + 1)th iteration is the weighted sum of the energy expectation value of these states: minimize ⟨Σ n (θ)| Ĥ|Σ n (θ)⟩ with respect to θ, 11: n ← n + 1. 14: end while Thus, at each iteration, the ensemble energy through order η is As in the case of the ground-state calculation [18], the gradient of the ensemble energy can be computed with respect to each A (n) pq,st : where r ν ⟩.This shows that the residual of the energy is the weighted expectation value of the commutators [ Ĥ, Γpq st ].The residual goes to zero when the ensemble is composed of eigenstates, which means that the ACSE in Eq. ( 7) is fulfilled.Hence, an algorithm to find the optimal operator Â using gradient descent should perform the following update of the parameters at each step: which implies that η is the learning rate of the algorithm.Interestingly, the purification introduced in Eq. ( 6) can be used to write a more compact expression for the residual of the ensemble ACSE in Eq. ( 10), namely: If, in addition, one chooses the auxiliary states as a replica of the physical ones (i.e., |a ν ⟩ = |ϕ ν ⟩), then the state can be written as the following unitary transformation of the system's vacuum [36]: |ρ(w)⟩ = V (w)|0⟩, where V (w) = U D(w), U is a unitary acting on the physical space and D(w) is a squeezed operator acting on the duplicate Hilbert space.As a result, the total residual can be written as a vacuum expectation value: ⟨0|[ Γpq st (w), Ĥ(w)]|0⟩, where the notation Â(w) = V † (w) ÂV (w) is used.
One possible way to implement the ACSE in a quan-Algorithm 2 Weighted random CQE 1: Given K > 0, w = (w 0 , .., w K−1 ), ν w ν = 1, δ > 0, 2: choose 0 < η < 1, and N > 0 number of shots, 3: choose K initial states {|ϕ end for end for 15: prepare |ϕ n ← n + 1. 18: end while tum device is to choose w ν as fixed quantities and, for the (n + 1)th iteration, allocate a certain number of shots N ν to measure the contribution of r ν;pq,st to the total residual in Eq. ( 10).Yet it is known that the most efficient way of deterministic assigning shots among the measurements consists of allocating N ν proportionally to w ν [37,38].But since the weights are not integers, this assignment results in a "hard floor" on N total = ν N ν ≥ 1/w K (recall that w K is the minimum of the weights) [39].This is the minimal number of shots needed for an unbiased estimate of the residuals of the ensemble ν w ν r (n) ν;pq,st .Unfortunately, for large K one can expect quite small w K and therefore very large numbers of shots for each unbiased estimate.Random sampling can efficiently perform unbiased estimations of the residuals of the ensemble energy in Eq. ( 7) while using a cheap number of shots.In the next section, based on this sampling, we will present two quantum algorithms.

D. CQE for excited states
To introduce our algorithms, let us start first by choosing a set of weights w, which for convenience we normalize to 1: ν w ν = 1.Next, we choose K initial orthogonal states |ϕ ν ⟩ that can be the K lowest meanfield (Hartree-Fock) wave functions.Weighted random sampling, where the probability of measuring r ν is proportional to w ν , can be used as an efficient unbiased estimator of the residuals of the ensemble energy in Eq. (7).A promising alternative to implementing this anti-Hermitian CQE that does not require a random number generator consists of preparing and measuring the purification presented in Eq. (6).For this parallelized CQE the initial state in Eq. ( 6) can be prepared by applying a suitable linear combination of unitaries [40] to the original Hartree-Fock state |ρ 0 ⟩ = |ϕ HF ⟩ ⊗ |0, ..., 0⟩, with an ancilla term that uses only log 2 K qubits.At each iteration, the states |Λ ± n ⟩ = exp(±iη Ĥ)|ρ n (w)⟩ are prepared and the entries of the matrix A (n) are measured from the equation Importantly, the residual in Eq. ( 10) is exactly zero for any set of eigenstates, not necessarily the lowest ones, so for any combination of eigenstates the optimization will stop at this point.Hence, to guarantee that the lowest set is found, we further prepare the state |Σ n (w)⟩ = exp(θ Â(n) )|ρ n (w)⟩ and minimize the ensemble energy with respect to the value of θ.Besides circumventing local minima, this will also guarantee a faster convergence.As described in Algorithm 1, the process is repeated until a desired convergence is reached.We also sketch the weighted random CQE in Algorithm 2. This algorithm follows similar lines as Algorithm 1 except that the purification (or the parallelization) is replaced by assigning m ν number of shots per state |ϕ ν ⟩ randomly from a multinomial distribution: m ∼ Multinomial(N total , w).Because each of the excited states is treated separately, the algorithm is amenable to distributed parallel programming in which each state is prepared and measured on a separate quantum processor with the results only collected for the classical parts of the optimization.This weighted random sampling algorithm is equivalent to measuring the expectation value of the pure state |ρ(w)⟩, in the sense that the variance of any observable computed by both methods does coincide.As a result, the number of shots needed to achieve a certain measurement error of the residuals is the same for both algorithms.Yet, while the results are certainly the same, the implementation clearly differs in the requirement of computational resources.An advantage, however, of the purification lies in the fact that quantum symmetries can easily be added to the cost function to improve convergence [41,42].

III. RESULTS
We now present the results of both Algorithms when applied to model and molecular Hamiltonians and discuss their advantages and disadvantages.The first system we investigate with the ensemble ACSE is the generic Mqubit Hamiltonian: FIG. 1: Evolution of the projection of the states |ϕ where σ n denotes the Pauli matrix.The initial state is denoted as |ρ 0 (w)⟩ = i∈{0,1} M √ w i |i⟩ p ⊗ |i⟩ a , where p/a denotes the physical/ancilla qubits and i = (i 1 , ..., i M ).The evolution into the exact eigenstates for a random Hamiltonian of the form in Eq. ( 12) for systems sizes M = 2, 3 is presented in Fig. 1.We chose the learning rate η = 0.3 and weights w = (M 2 , M 2 − 1, ..., 1) and then w → w/ i w i .For M = 2, the ground state is reached in 8 iterations, while the exact eigenstate calculation is reached in 20.The highest energy states, having the lowest weights in the cost function, converge the slowest, and, due to orthogonality limiting the degrees of freedom, converge simultaneously.A similar pattern can be seen for another random Hamiltonian for the case M = 3 but due to the larger dimension of the Hilbert space, more iterations are needed for convergence.We investigate also two molecular examples: a noisy backend simulation of H 2 and a noiseless state-vector simulation of H 4 .All calculations were performed using the minimal Slater-type orbital (STO-3G) basis set.The noisy backend is the FakeLagosV2 by IBMQ.
Calculation of H 2 is performed in the spin-symmetry sector S z = 0. Based on the symmetry of the problem, we construct the Hamiltonian in a compressed form with two qubits.Two additional ancillary qubits are used to create the purified ensemble of all four eigenstates, resulting in four qubits in total.The detailed circuit preparation has been reported in previous work [43].For the singlepoint calculation in Fig. 4, performed with the paralleled CQE, the ensemble energy converges to a minimum in only three iterations.Remarkably, we achieve an error of less than 30 mHartree for each state without any error mitigation techniques.We also present the dissociation curve of H 2 in Fig. 4. Energies computed from parallel CQE are in excellent agreement with the full CI results with an average mean unsigned error of 26 mHartree.
It is also worth discussing the role weight values w i play in the rate of convergence.For instance, if all of the weights are equal, only an eigen-subspace can be found, and the individual eigenstates would have to be resolved with classical diagonalization.Giving different values for the weights allows us to perform the entire calculation on a quantum device, resulting in a faster convergence.Indeed, we find that the optimal convergence for H 2 (presented in Fig. 4) is achieved with the weights w ∼ (9, 9, 1, 1), before normalization.To explain our choice, let's observe that, due to system's point-group symmetry, the Hamiltonian matrix is block diagonal with two 2 × 2 sub-matrices on the diagonal.Therefore, since the minimization runs independently within each subblock, we opted for two identical pairs of weights.This indicates that the optimal choice of weights is highly dependent on the molecular symmetries.Linear H 4 is a widely used benchmark system for Energy (Hartree) FIG. 5: The computed and exact lowest eight eigenenergies of the linear equidistant H 4 as a function of the H-H distance.
As the molecule dissociates, the energy levels become highly degenerate due to the non-interacting hydrogen atoms and the system exhibits significant static correlation [45].We take the equidistant form of linear H 4 and use the Jordan-Wigner transformation to map the Hamiltonian from four spatial orbitals to eight qubits.Both algorithm 1 and 2 successfully find the ground and excited states.Yet in the first case, as we are tackling eight states simultaneously, one requires at least three ancillary qubits to prepare all initial states in the expanded Hilbert space.Alternatively, preparing different initial states separately and sampling them using a multinomial distribution (as in our Algorithm 2) becomes particularly valuable with limited qubit resources or when the ancillabased preparations are hard to perform.For the calculation of H 4 shown in Fig. 5, we have used the weight vector (8, 7, ..., 1), before normalization.At a long bond distance, we seed the eight initial guesses with the eight single Slater determinants with the lowest energies.Afterward, each state in the calculation is seeded with the two most important Slater determinants of the corresponding state found in the previous calculation.While the potential energy curves are highly degenerate towards dissociation, as the bond begins the form, the energy curves separate.As shown in Fig. 5, for the dissociation curve on a noiseless simulator, our algorithms give almost exact results (i.e., an error of around 10 −4 Hartree).Most calculations converged in less than 200 iterations.We recall that this convergence speed does depend on the weight being assigned to each element, the initial guess, as well as the optimization method, suggesting opportunities for further exploration and improvement.

IV. CONCLUSIONS
In this paper, we have combined the contracted quantum eigensolver (CQE), originally developed for the calculation of molecular ground states, and the Rayleigh-Ritz variational principle for ensemble states into an excited-state CQE.Quite remarkably, our scheme allows us to compute simultaneously an arbitrary number of lowest eigenstates while preserving the favorable scaling and ease of implementation of the ground-state CQE.Unlike approaches based on the unitary coupled cluster and related ansätze, that give an approximation to the cost function, our algorithm provides a natural choice for the unitary operator through the measured residual.In our experiments with molecular and model systems, we tackle multiple states simultaneously with excellent accuracy in both the weakly and strongly correlated regimes.The ability to optimize near-degenerate states by assigning different weights allows us to study both neardegeneracy and conical intersections, which can be used for nonadiabatic chemistry.Another interesting question for the future is how to use our algorithms for excited states and spectroscopy when additional bosonic degrees are present.
Code availability.-Allcodes to reproduce, examine, and improve our proposed analysis are freely available online [46].

3 FIG. 4 :
FIG. 4: (a) Obtained energies during the optimization for single point calculation of H 2 (bond distance of 0.7 Å).The exact solutions for each state and the ensemble are indicated by black dashed lines.(b) Energies along the dissociation curve computed from 0.5 to 5 Å.The exact results are shown as black lines and our single-point calculations are shown as dots.