Folded Spectrum VQE: A Quantum Computing Method for the Calculation of Molecular Excited States

The recent developments of quantum computing present novel potential pathways for quantum chemistry as the scaling of the computational power of quantum computers could be harnessed to naturally encode and solve electronic structure problems. Theoretically exact quantum algorithms for chemistry have been proposed (e.g., quantum phase estimation), but the limited capabilities of current noisy intermediate-scale quantum devices motivated the development of less demanding hybrid algorithms. In this context, the variational quantum eigensolver (VQE) algorithm was successfully introduced as an effective method to compute the ground-state energies of small molecules. This study investigates the folded spectrum (FS) method as an extension of the VQE algorithm for the computation of molecular excited states. It provides the possibility of directly computing excited states around a selected target energy using the same quantum circuit as for the ground-state calculation. Inspired by the variance-based methods from the quantum Monte Carlo literature, the FS method minimizes the energy variance, thus, in principle, requiring a computationally expensive squared Hamiltonian to be applied. We alleviate this potentially poor scaling by employing a Pauli grouping procedure to identify sets of commuting Pauli strings that can be evaluated simultaneously. This allows for a significant reduction in the computational cost. We applied the FS-VQE method to small molecules (H2, LiH), obtaining all electronic excited states with chemical accuracy on ideal quantum simulators. Furthermore, we explore the application of quantum error mitigation techniques, demonstrating improved energy accuracy on noisy simulators compared with simulations without mitigation.


I. INTRODUCTION
Computing the electronic structure of molecules and materials is crucial for the prediction of chemical or structural properties.Theoretical chemists and physicists have acknowledged the essential challenges that must be addressed, but the exponentially scaling dimensionality of electronic structure problems seems to be insurmountable on classical computing resources.As a result, the theoretical study of large molecules or materials using ab initio methods such as coupled cluster is impractical.Hence, less costly methods involving approximations are generally employed at the cost of a loss in accuracy and predictive power.The emergence of quantum computing presents potential novel pathways for theoretical chemistry, as quantum resources show exponentially scaling computational power that could be harnessed to naturally encode and solve quantum problems.Although exponential speedup may not be achieved, a polynomial acceleration could be groundbreaking for quantum chemistry applications [1].
In this context, the Variational Quantum Eigensolver (VQE) [2] was introduced as an effective algorithm to find the lowest eigenvalue of a quantum observable.In particular, it can compute the ground-state energy of a molecular Hamiltonian.The capability of VQE for electronic ground-state computation of small molecules has been extensively studied [3], but the effective and direct computation of excited states remains elusive.
In this study, we propose a variant of VQE that aims at computing molecular excited states.It uses the Folded Spectrum (FS) method to reorder the Hamiltonian's eigenspectrum, thus allowing for the direct com-putation of highly excited states.Although this method is documented in the literature, its quantum implementation was considered too costly due to the squared number of terms of the measured operator [4,5].Here, we show that a Pauli grouping procedure reduces the required number of measurements, thereby making the cost of the FS method reasonable.The effect of Pauli grouping is particularly significant for second quantized molecular Hamiltonians as a result of their particular structure.Finally, we present FS-VQE results obtained on a noisy quantum simulator, and show the successful use of quantum error mitigation techniques on this algorithm.

II. VARIATIONAL QUANTUM EIGENSOLVER
The Variational Quantum Eigensolver [2] is a hybrid quantum-classical algorithm (see figure 1).Its purpose is to find the lowest eigenvalue of a given quantum operator.It can be applied to quantum chemistry problems to obtain the electronic ground-state of a molecule, by focusing on the molecular Hamiltonian Ĥ.
The algorithm relies on an ansatz to prepare a trial electronic wavefunction on the quantum computer.The ansatz takes the form of a parametric quantum circuit whose parameters, denoted θ, are angles in rotation gates.Details on the ansatz design are given in section II A.
The quantum subroutine prepares a parametric trial wavefunction Ψ(θ) on a qubit register.This quantum state can be assessed by measuring the qubits : from the measurement results, the expectation value of the molecular Hamiltonian ⟨ Ĥ⟩ can be computed on a classical computer (see section IV A).This value corresponds to the average electronic energy of the trial wavefunction Ψ(θ).A classical optimizer is then used to adjust the parameters θ in the ansatz in order to minimise the value of ⟨ Ĥ⟩.By means of the variational principle, the minimal expectation value obtained for a set of parameters θ opt is an upper bound on the Hamiltonian's ground energy.The quantum state prepared with the optimal, final angles θ opt is a representation of the molecule's ground-state electronic wavefunction.

A. Ansatz
In the context of VQE, the ansatz is a parametric quantum circuit aiming to explore the wavefunction search space.The ansatz design can take various forms as different properties are targeted [6].The number of parameters in the circuit is key to the success of the optimization procedure ; a very large number of parameters may lead to intractable optimization.
The so called chemically motivated ansatz class includes ansätze inspired by quantum chemistry methods [6].Their advantage is that the prepared states are by design physically relevant (number of electrons and total spin are conserved).However, they often require a large number of parameters and deep quantum circuits, which limits both the optimization success and their feasibility on NISQ hardware.
Another approach is to design hardware motivated ansätze [6].Such ansätze are constructed to be efficiently implemented on quantum computers.Strong constraints in terms of quantum gates, qubit connectivity, number of two-qubit gates, global circuit depth, etc. are defined in accordance with the capability of the target hardware.These ansätze are computationally advantageous, but they do not offer a guarantee on the physical properties of the prepared trial states, thus limiting the convergence.

Unitary Coupled Cluster Ansatz
Unitary coupled cluster (UCC) is a widely used chemically motivated ansatz for electronic wavefunctions in quantum computing.It is a unitary variant of the wellknown coupled cluster (CC) theory.Like coupled cluster, UCC is based on a reference wavefunction (often Hartree-Fock) and it creates linear combinations of excited determinants using excitation operators T : where T1 is the operator of all single excitations, T2 the operator of all double excitations, etc. â † k and âk are, respectively, the fermionic creation and annihilation operators acting on orbital k.Indices i, j denote occupied orbitals and a, b virtual orbitals.Parameters θ are optimized to obtain the CC wavefunction.
Because the CC operator e T is not unitary, it cannot be directly implemented on a quantum circuit.To create a unitary variant of CC, the cluster operator needs to be modified to become anti-hermitian, as the exponentiation of an anti-hermitian operator is unitary: Therefore, the anti-Hermitian cluster operator T − T † is considered.The UCC ansatz state is created similarly to the CC state : For the UCC state to serve as a quantum ansatz, the operator e T − T † must be expressed in terms of quantum gates.Since excitation operators do not commute, a Trotterization step is required to decompose the exponentiated operator [7].The Trotter decomposition is given by : with Ti the excitation operators defined in eq. ( 2) and (3) and p the Trotter decomposition order.In this work, the order p = 1 was implemented, as it was proven to be an exact and general form of UCC under the condition of an appropriate ordering of excitations [8] (details in section II A 2).This large form of UCC is truncated to a rank k corresponding to the highest excitation considered (for example, UCCSD is k = 2 : only single and double excitations).The resulting smaller ansatz is, therefore, : with parameters θ to optimize.The trotterized and truncated UCC operator can be translated into quantum gates in two steps.The first step is to map the â † and â of excitation operators into Pauli strings (tensor products of Pauli matrices, acting on several qubits) as described in section II B. Then, each exponentiated Pauli string can be translated into a Pauli gadget [9], as shown in figure 2. The ansatz parameters θ are the rotation angles in the gates R z , denoting single qubit rotations around the Z axis.The sequence of Pauli gadgets gathered into one single circuit constitutes the UCC ansatz.

Ordering of excitation operators
The Trotterized form of UCC is a product of noncommuting terms, making the ordering of excitation operators an important hyper-parameter of the ansatz.In Ref 8, authors have proposed a universal ordering of the excitations allowing us to reach any fermionic state .
Their method is employed in this study.Considering a single determinant reference state |Φ 0 ⟩ (chosen as the Hartree-Fock wavefunction here), we iterate through the occupied indices j of |Φ 0 ⟩.For each index, all single excitations involving index j (e θ a j â † a âj ) are added to the ansatz, followed by all double excitations with index j (e θ ab ij â † a â † b âiâj ) and so on for all excitations of higher rank involving orbital j.This procedure is repeated for all occupied indices j of the reference state, eventually adding all excitations of the UCC ansatz to the circuit.

B. Fermion-to-qubit mapping
Second quantized fermionic operators can be mapped to qubit operators, implementable on a quantum circuit.
Different mapping schemes are available, the most common being the Jordan-Wigner (JW) [10] and Bravyi-Kitaev (BK) [11] transformations.In the second quantized formalism, fermionic operators are expressed as sums of creation and annihilation operators.Fermion-to-qubit mapping is a systematic formula to translate creation and annihilation operators into Pauli strings.Any second quantized operator : can be mapped to : with Pi being Pauli strings and o i scalars.
In this work, Jordan-Wigner mapping was employed.

Jordan-Wigner Mapping
In this formalism, each qubit represents a fermionic state (a spin-orbital for molecules), with the qubit state |0⟩ corresponding to an unoccupied state, and |1⟩ to an occupied state.The creation and annihilation operators are mapped using the transformation in equation (10) for a N -qubit register corresponding to N electronic spin-orbitals.
with gates applied to qubit j.The σ + and σ − gates act as qubit creation and annihilation operators, while the Z gates are required to conserve the anti-commutation relations : C. VQE procedure The procedure implemented in this work is summarized in figure 3.
•Select a basis for the system •Build the trial wavefunction on the quantum register, by applying the UCC ansatz to the HF state.
• Classically compute : •The classical optimizer updates

Repeat until convergence
•Build the final wavefunction with the angles measure the qubits and compute : •The final energy is ...

Repeat (shots)
for ground state for excited states using the integrals and and the counts.

III. EXCITED STATES
The standard VQE algorithm applied to molecular systems allows one to compute the ground-state of the electronic wavefunction.It is not primarily designed to compute excited states, as it relies on the minimisation of the average energy.Several approaches have been proposed to reach excited states with quantum algorithms, including Quantum Subspace Expansion (QSE) [12], Variational Quantum Deflation (VQD) [13] similar to Orthogonally Constrained VQE (OC-VQE) [14,15], witnessing eigenstates (WAVES) [4] or Quantum Equation of Motion [16,17].
The Folded Spectrum (FS) method has also been reported in the literature [6], but the presence of a quadratic term in Ĥ is regarded prohibitive and it is expected to scale as O(N 8 ) relative to system size [14].To the best of our knowledge, no extensive study of this method has been reported.

A. Folded Spectrum method
The principle of the FS method is to minimise the expectation value of the FS operator ( Ĥ −ω) 2 instead of the Hamiltonian Ĥ, with ω an arbitrary target energy.
This method is also known as state-specific variance minimisation in the Quantum Monte Carlo (QMC) literature, where it has been actively employed and studied for many years [18][19][20].
Let |Ψ⟩ be an eigenstate of the Hamiltonian Ĥ.It satisfies the time-independent Schrödinger equation : The linearity of the Schrödinger equation allows one to write equation (13) for all |Ψ⟩ eigenstates of Ĥ, E the associated eigenvalues, and ω an arbitrary scalar.
The FS operator ( Ĥ − ω) 2 and the Hamiltonian Ĥ share the same eigenstates but with a reordering in the eigenvalues (corresponding to a fold around ω) [21].The lowest lying eigenstate of the folded operator is the one with an energy E i closest to ω (see figure 4).
By minimising the expectation value of the FS operator ( Ĥ − ω) 2 , one can find an eigenstate of Ĥ such that (E i − ω) 2 is minimal, and thus obtain an excited state of the Hamiltonian, close to the target energy ω.In practice, we perform expectation value minimization with a VQE procedure, and the obtained wavefunction is an approximation of the true eigenstate given by the ansatz.The cost function of interest is : All excited states of the Hamiltonian may be obtained by modifying the parameter ω over a wide enough range of energies.
One major limitation of the FS method is that it requires a squared Hamiltonian, containing a large number of fermionic terms and thus of Pauli strings compared to the Hamiltonian itself.However, by using a Pauli reduction and grouping procedure as described in section IV B, the number of required measurements can be considerably reduced.

spectrum
Folded spectrum Colored lines represent eigenstates along a vertically ascending axis.The eigenspectrum of Ĥ (left) is folded around ω in the spectrum of ( Ĥ − ω) 2 (right), causing its eigenvalues to reorder.The lowest eigenvalue of the folded spectrum, circled in red, is the excited eigenstate of Ĥ originally closest to ω.

A. Computing expectation values
The expectation value of quantum operators represented by Pauli strings (as in equation ( 9)) can be decomposed as shown in equation ( 15), with Pi Pauli strings and o i scalar coefficients.
The expectation value of Ô can be computed by classically summing the expectation values of each Pauli string ⟨ Pi ⟩ weighted by the pre-computed coefficients o i .In molecular Hamiltonians the coefficients o i are typically denoted by h i and are formed through linear combinations of one-body and two-body integrals.
Once the trial state |Ψ⟩ is prepared on the quantum register, we measure the qubits and repeat the state preparation and measurement procedure several times (shots).Ultimately, we obtain some counts that form estimates for |Ψ⟩ populations.Qubit measurements are usually performed in the computational basis denoted {Φ i } 2 n i=1 , corresponding to the values 0 or 1 for each of the n qubits : {Φ i } = {|0...00⟩ , |0...01⟩ , ..., |1...11⟩} in binary order.In this basis, the spectral decomposition of |Ψ⟩ is : with Measurable quantities are the populations for each basis vector of the computational basis : Note that the counts provide estimates of the populations, due to finite sampling.The final precision ϵ in the results is directly correlated with the number of shots taken s as ϵ ∼ 1 √ s (see Appendix 1).Results can be made arbitrarily close to the theoretical value by increasing the number of shots, but this can lead to considerable computing time.To reduce the number of quantum measurements, a Pauli grouping routine can be used, as discussed in the next section.
The derivation of expectation values from the counts results is explained below for diagonal and non-diagonal operators.

Diagonal Operators
The expectation value of diagonal Pauli operators in the computational basis (i.e.tensor products of I and Z) can be directly computed from the counts.Such operators can be decomposed in the computational basis as a sum of projectors : with λ i corresponding to the eigenvalues of the operator, namely +1 or −1 for products of I and Z Pauli operators.
The expectation value of Pdiag is therefore : Consequently : is directly accessible from the quantum measurement, by classically summing the count results weighted by the eigenvalues of the Pauli operator.Note that all diagonal Pauli operators can be evaluated from the same counts measurement, as the λ i coefficients are treated classically.

Non-diagonal operators
Given a non-diagonal operator P in the qubit basis, it is always possible to find an appropriate basis change to diagonalize it.X and Y Pauli matrices are nondiagonal in the computational basis, but they can be diagonalized using the following basis changes: Therefore, for a Pauli string P , it is possible to find a rotated Pauli string P diagonal in the qubit basis, using these basis changes.In general, one can write : with The expectation value of P is then : with Finally, To evaluate the expectation value of a non-diagonal Pauli string P , it is therefore necessary to apply a postrotation gate R to the quantum state |Ψ⟩, which creates a state Ψ in a basis where the Pauli string is diagonal.In practice, the additional post-rotation operator is built with single-qubit gates added after the ansatz.A Hadamard is applied to the qubits where the Pauli operator is X, and a R x π 2 is applied to the qubits where it is Y .
The expectation value of a non-diagonal P is then computed similarly to a diagonal Pauli string, using the rotated state Ψ (equation ( 28)).

B. Pauli strings reduction and grouping
The number of Pauli strings in the Hamiltonian scales polynomially with the system size, and naively the FS operator can contain up to the square of this number.Evaluating each term one by one can lead to a very large number of measurements, which lowers the potential advantage of the quantum algorithm.

Pauli reduction
When computing the FS operator ( Ĥ −ω) 2 , the number of Pauli strings primarily obtained is approximately the square of the number of terms in Ĥ.It is possible to simplify and reduce this sum by using the commutation and anti-commutation relations between Pauli matrices [22,23].In Ref 24, authors studied a collection of systems of increasing sizes and concluded that the actual number of Pauli strings in Ĥ2 after Pauli reduction has an effective scaling below O(N 6 ) instead of the expected O(N 8 ) with N the number of spin-orbitals.This result can be extended to our work : the number of terms in the FS operator has a much more favorable scaling with respect to the system size due to Pauli reduction.More formal analyses are required to consolidate this result and assess the feasibility of the FS method for larger systems.

Pauli grouping
To further reduce the number of quantum evaluation required, one can partition the operators into groups of simultaneously diagonalizable Pauli strings [25,26].All the Pauli strings in the same group can have their expectation values determined with a single quantum evaluation, by adding a classical post-processing step.In the formalism of section IV A 2, it means that all Pauli strings in the group share the same post-rotation R in equation (24).Formally speaking, a group of operators is simultaneously diagonalizable if and only if the operators commute [27].This reduces the problem to identifying groups of commuting Pauli strings in the qubit Hamiltonian or FS operator.In particular, we want to find a partitioning with a minimal number of groups, leading to a minimal number of quantum evaluations.
Two distinct definitions of commutation can be considered to partition the Pauli strings: qubit-wise commutativity (QWC) or general commutativity (GC) [28].The former defines that two Pauli strings commute if the Pauli matrices commute at each index.For instance the group {IX, XX, XI} is QWC since all Pauli matrices for qubit 1 commute, and equally for qubit 2. General commutativity is fulfilled if the two Pauli strings commute, regardless of the single-qubit case.The group {XX, Y Y , ZZ} is GC although none of the pairs is QWC.The general rule is that each pair must fail to commute at an even number of indices.QWC is, in fact, a special case of GC where the strings fail to commute at 0 indices.Figure 5 shows an example of QWC and GC partitioning for the electronic Hamiltonian of H 2 .Finding the optimal Pauli partitioning (in QWC or GC) is equivalent to a graph partitioning problem known as the minimum clique cover problem [29], and it is NP hard [30].Efficient heuristic algorithms to find a good Pauli partitioning are therefore essential to tend toward scalability for the FS method [31,32].
In this work, we use QWC Pauli partitioning as it is more straightforward to implement.We present some results using Jordan-Wigner qubit mapping in figure 6, and table II in appendix 2 reports some examples for both J-W and B-K transformations.The number of quantum evaluations is systematically decreased by grouping the Pauli strings.As expected, more evaluations are required for the FS operator compared to the Hamiltonian for the same system.Additional results on Hamiltonian grouping for other transformations and systems can be found in Refs 31,33,34.As shown in Ref 28, GC partitioning is more efficient than QWC, and it would lead to even fewer quantum measurements for the FS method.

Measurement cost
Pauli grouping for molecular Hamiltonians shows interesting features that may be extended to the FS operator.Electronic Hamiltonians under certain fermion-toqubit transformations (such as JW mapping described in section II B) have the convenient property of including a large number of diagonal Pauli strings (with only I and Z operators).In fact, the product of σ + and σ − operators defined in equation ( 10) can be decomposed as shown in equation (29).
Therefore, all terms in the Hamiltonian involving a creation and an annihilation operator for the same spinorbital will be mapped to a diagonal Pauli string in the qubit basis.These terms are the one-body terms with p = q (h pp ) and the two body terms with p = r and q = s h prpr or p = s and q = r (h prrp ).For a Hamiltonian describing n electrons in N spin-orbitals, there are N diagonal one-body operators and N 2 diagonal twobody operators.All these diagonal Pauli strings can therefore be grouped together and evaluated simultaneously.
Additionally, asymptotically dominant terms in the molecular Hamiltonian are two-electrons operators of the form h pqrs â † p â † q âr âs with p ̸ = q ̸ = r ̸ = s [28].The number of such terms in a molecular Hamiltonian describing N spin-orbitals scales as N 4 ∼ O(N 4 ).Under Jordan-Wigner mapping, two of these terms with disjoint indices, namely â † p â † q âr âs and â † i â † j âk âl with {p, q, r, s} ∩ {i, j, k, l} = ∅ involve disjoint qubits for X and Y gates, and authors of Ref 28 demonstrated that they commute.They showed that using the Baranyai's graph coloring theorem [35], it is possible to partition these Similarly to the Hamiltonian, the FS operator can be partitioned into commuting groups to reduce the number of measurements.Here, we present empirical data on the effect of Pauli grouping for the FS operator.Future studies could aim to establish analytical results on the scaling of the number of evaluations needed for the FS method.

C. Reference state
When computing the molecular ground-state, the Hartree-Fock (HF) determinant can often be used as a reference because of the significant overlap between the HF state and the FCI electronic ground-state.This is not generally true for excited states.In this case an excited single determinant or a superposition of two or more determinants having overlap with the target wavefunction can be used as a reference.
In this study, we have selected relevant references for all electronic excited states by exciting the ground HF determinant with single and double excitations and symmetrizing the spin function when necessary.This procedure can be generalised to larger systems, but we expect that more sophisticated reference states may be required for molecules with strong multi-reference character.This question remains an essential challenge for the scalability of the FS method.The reference states used in this paper are specified in the results section.

V. NOISE ROBUSTNESS AND ERROR MITIGATION
The presented algorithm is based on the Variational Quantum Eigensolver, which is designed to be amenable on near-term quantum devices.To assess the feasibility of our method on noisy devices, we evaluated its noise robustness by including noise models in our simulations and using zero noise extrapolation (ZNE) and state preparation and measurement (SPAM) error mitigation techniques [36].

A. Noise Model
We designed a tunable and realistic noise model based on the models provided for IBM quantum devices [37].This model includes: • gate errors consisting of a depolarizing channel characterized by one and two qubit gates error rates p 1 and p 2 respectively, followed by thermal relaxation and dephasing processes, driven by T 1 and T 2 characteristic times applied for the gate lengths t gate1 and t gate2 for one and two qubits gate.
• readout error characterized by an error probability p SPAM for each qubit.
This model can be considered realistic for quantum computers, provided that the parameters are adjusted to the calibration data of the real device.However, it excludes noise sources such as state leakage or cross-talk which are more complex to model.This approximation is widely adopted for its simplicity, but may not be suitable for some devices where leakage or cross-talk effects are not negligible [38].
Parameters T 1 , T 2 , t gate1 , t gate2 p 1 , p 2 and p SPAM are provided for IBM machines.Some typical values for the present devices are given in Table I [39].Table I: Values of the noise parameters for different values of the scale factor λ. λ = 1 corresponds to the order of magnitude of the experimental values in IBM machines [39], and λ = 0 is the ideal value in a noiseless device.In practice, ideal values were fixed at T 1∞ = 2000µs and T 2∞ = 1000µs.
To vary the noise level, we scaled the parameters according to a scale factor λ between their experimental value (λ = 1) and their ideal value (λ=0).The scaling is performed exponentially for T 1 and T 2 and linearly for all other parameters.Ideal values for T 1 and T 2 were fixed at two orders of magnitude above the total length of the circuit, that is, T 1∞ = 2000µs and T 2∞ = 1000µs here.

B. State Preparation and Measurement mitigation
The SPAM technique aims to mitigate errors introduced during the state preparation and measurement stages [40,41].A confusion matrix A is measured by running small calibration circuits on the quantum processor.A represents the noise channel corresponding to the probability that the state preparation or measurement outcome will be incorrect for each qubit.This matrix is then inverted, and the inverse channel A −1 is classically applied to the next experiments as a postprocessing step, thus obtaining quasi-probabilities with mitigated SPAM error.The nearest probability distribution is then selected as the mitigated measurement result.This method assumes that the SPAM noise remains constant over multiple experiments close in time on the same device.The calibration circuits should be run regularly to make this assumption reasonable.A more detailed description of the SPAM error channel can be found in Ref 38.In our implementation, the confusion matrices are directly extracted from the noise model's parameters p SPAM for each qubit.

C. Zero Noise Extrapolation
Zero noise extrapolation (ZNE) is a technique that aims to mitigate the effect of noise when evaluating expectation values on QPUs [42].In ZNE, the hardware noise level is represented by a parameter γ, with γ = 1 corresponding to the actual noise level of the quantum computer, γ > 1 being a noisier hardware and vice versa.
The principle of ZNE is to intentionally increase the noise level (γ = 3, 5, 7...), and to evaluate the same expectation value E γ for the different values of γ.The points obtained are plotted on a E γ vs γ curve and fitted with an analytic model.The model provides an extrapolated value at γ = 0 that is retained as the mitigated result E 0 , which should be an approximate of the ideal value E.
In our implementation, ZNE was employed to mitigate the summed expectation values of each Pauli groups separately.For H 2 , the Folded Spectrum operator contains 24 Pauli strings that can be partitioned using general commutativity into two commuting groups G 1 and G 2 : For this system, each evaluation of ⟨( Ĥ −ω) 2 ⟩ requires measuring two different circuits, one for each group.A ZNE procedure can be used to mitigate the noise on the summed expectation value of each commuting group.
To increase the noise level γ, we used a method named unitary folding that consists of replacing a unitary gate sequence G by a folded version (for example, GG † G).The logical operations of G and its folded versions are the same as GG † = 1, but the effective number of gates in GG † G is three times greater, corresponding to γ = 3.Further folding is performed to implement γ = 5, 7, 9.... On noisy devices, folding the circuits effectively corresponds to scaling the gate noise.Thus, this method allows us to artificially implement different values of γ by increasing the circuit depth.This method has the advantage of being very general since it can work for any unitary circuit.
A quadratic model in equation ( 31) was selected to fit and extrapolate the E γ (γ) curve.
ZNE relies on the assumption that folding the circuits corresponds to increasing the noise level by γ overall, which implies that gate noise is the main source of noise.This assumption is not always valid, and in particular, SPAM error is not scaled with circuit folding, which can make the extrapolation process erroneous.To alleviate this limitation, we employ ZNE mitigation in conjunction with SPAM mitigation, such that in theory SPAM errors are negligible in the ZNE fitted data.
An example of this procedure is given in figure 7. (bottom).Noise scaling is performed by unitary folding with folding factors 1, 3, 5 and 7, and extrapolation is performed using a quadratic model.

VI. COMPUTATIONAL DETAILS A. State tracking
Both the energy and the wavefunction show continuity along the potential energy surface (PES) of the same electronic state.Therefore, the final energy and final angles for one molecular geometry are good starting points for another close molecular geometry along the PES.We take advantage of this property by setting ω to the previous energy computed on the PES, as well as setting the initial θ parameters to the angles found in the previous calculation.In other words, with k index representing the evolution along the PES, sets the initial trial wavefunction for point k + 1 as the final wavefunction of point k, and sets the target energy for the next point to the previous energy computed.
This state tracking method allows one to reduce the optimization time.To ensure continuity of the wavefunction along the PES, it is necessary to have continuous molecular orbital (MO) coefficients as the geometry changes.The MO coefficients are pre-computed for each geometry using the PySCF package [43].Phase jumps are possible in the RHF computation, as the MO phase can freely change between independent calculations.These phase jumps do not affect the energy, but they do break the continuity in wavefunction, lowering the effectiveness of state tracking.To avoid phase jumps in the MO coefficients, we compute at each step between close geometries k and k + 1 : with C k the MO coefficient matrix at geometry k and S k the overlap matrix at geometry k.The P matrix has a diagonal with ≈ ±1 elements.Negative signs indicate a phase jump between the two geometries.In this case, we rectify the phase of the corresponding MOs in C k+1 matrix, and use the rectified MOs in the calculation.This ensures continuity in the ansatz parameters and allows us to facilitate convergence and reduce computation time by employing the state tracking method.
We observed that state and energy tracking, beyond reducing computation time, can also help the optimization convergence for points where it initially fails.When noticing non-converged points on the PES, one can start from a previous converged point and state track towards the desired geometry with smaller geometry steps.This technique usually allows one to obtain better convergence.However, it requires a large number of calculations, since the step size needs to be small enough to allow good continuity.

Preventing jumps between electronic states
To some extent, state tracking helps prevent jumps between close electronic states as the wavefunction tends to be continuously evolved along the PES, which is particularly useful when degenerate or quasidegenerate states are present.However, we observed that jumps still occurred in our computations when the ω parameter was closer to another electronic state for a particular geometry (which is particularly frequent in the presence of large energy gradients or when two electronic states are very close in energy).This behavior is expected for the FS method, but it can be undesirable when trying to follow a particular electronic state on the PES.To prevent jumps, a continuity constraint term can be added to the optimized cost function, ensuring continuity of the parameters along the PES : (35) with η a scaling factor that can be adjusted to balance the relative weights of the two terms in the cost function.This regularization method was used for some isolated points in the PEC, and η was fixed to 0.1 in our implementation, based on empirical trials to find a value that influences optimization toward a region close to the previous calculated point, while still maintaining the right landscape and minima.

B. Classical optimization
The optimization of the variational parameters θ is performed using a classical optimizer.The dimensionality of the ansatz and the presence of noise in the cost function make the optimization difficult.
In addition to the finite sampling error that is inherent in quantum computing, NISQ hardware is characterised by the presence of noise within the quantum circuit, making robustness an important feature of quantum algorithms to be applied to current hardware.
Noise-tolerant optimizers are therefore the most appropriate.We use the SPSA (Simultaneous Perturbation Stochastic Approximation) optimizer for QASM (Quantum Assembly Language) simulations.In addition to being noise tolerant, SPSA also has a constant number of 2 evaluations per iteration that does not scale with the number of parameters [44].SPSA uses a stochastic procedure to update the parameters: at each iteration, the perturbation ∆θ is a randomly generated vector that has a component in every dimension of the optimization problem.The cost function is evaluated at (θ + ∆θ) and at (θ − ∆θ), and the numerical gradient for each parameter is approximated only from these two measurements.Because the perturbation vector is randomly generated, additional shifts due to noise in the cost function have a minor impact on the optimization process.The noise is, in a sense, absorbed by the stochasticity of the optimizer.

Scaling of parameters
In the UCC ansatz, the variational parameters are angles in rotation gates, ranging from −π to π.To assist the classical optimizer in finding the optimal angles, we scale the parameters to range from −cπ to cπ with c a predetermined constant.This simple procedure helps prevent optimization failures caused by optimization steps being too small to obtain non-zero numerical gradients, especially for QASM simulations.

Shots scheduler
Each evaluation of the cost function F θ is performed through a number s of state preparation and measurement procedures, where s is named the number of shots.The more shots s are taken, the more precise the mea-sured estimate of the cost function.The precision follows ϵ ∼ 1 √ s (see appendix 1).When measurement accuracy is not crucial, it is advantageous to use a smaller number of shots, as a large s requires a significant amount of computational time.For this reason, we use an increasing number of shots during the optimization, leading to uncertain measurements far from the optimum where the cost gradients are large and increasing precision when approaching the optimal parameters.For all computations, we used an inverse exponential scheduler of the form: with k > 0 that brings the number of shots from s min =1000 to s max =10000 with an exponential trend as iterations are performed.The final measurement after optimization convergence and post processing of the parameters (see section VI C) is performed with 30000 shots, allowing one to attain a better estimate of the final wavefunction and energy.

Quadratic fitting
The characteristics of the UCC ansatz search space can be harnessed to improve the parameters θ opt found by the optimizer.Let us consider the energy space in the UCCSD formalism, having dim(θ) dimensions.In this representation, the parameters θ ideal corresponding to eigenstates of the electronic Hamiltonian are located at minima or maxima in the energy space for the relevant excitations.More precisely, the eigenstates are located on vertices of parabolas in the energy space.
It is possible to take advantage of the particular location of the eigenstates to refine the solutions found by the optimizer.After optimization, one can probe the energy space around each parameter.If the solution found by the optimizer is close to an eigenstate, the energy space around each θ should be a parabola.It is therefore possible to sample a few points around the optimized solution, fit a quadratic equation, and choose the vertex of the fitted parabola as a refined solution.We employ this method as a post-processing step to improve the cost function.The refinement is usually very low (as the optimizer already locates the vertex of parabolas well enough), but in some cases a few tenths of a percent can be gained on the cost function.

Rounding of parameters
When the electronic wavefunction is built using the UCC ansatz, it is common for some excitations to be irrelevant because the wavefunction usually does not contain all possible Slater determinants in the molecule search space.As a result, some parameters of the quantum ansatz have an optimal value of zero, and therefore do not participate in the circuit.This can be harnessed to reduce the depth of the ansatz circuit as described in Ref 45.Similarly, some parameters can have an ideal value of π or π 2 when one determinant is completely excited to another, or two determinants have the same contribution to the wavefunction, respectively.This is particularly common when considering systems with internal symmetries.Here, we take advantage of this feature to improve the accuracy of the optimizer solution, by including a rounding post-processing routine.After the optimization (and quadratic fitting, see VI C 1) of the ansatz parameters, we detect close to zero (or close to a fraction of π) parameters, and evaluate the cost function when rounding those angles to zero (or to the corresponding fraction of π).If the cost function is improved by rounding them, the adjusted parameters are maintained.This procedure usually only improves the result very slightly, but the refinement can go up to a few hundredths of a percent in the cost function.

VII. RESULTS
The FS-VQE method was applied on two small molecules H 2 and LiH.The current capability of NISQ hardware (in terms of quantum volume and gate fidelity) is too limited for FS-VQE to obtain reasonable results on real quantum devices.Thus, we restricted our computations to small molecules and small active spaces that are tractable on simulators.
A. Excited states of H2 H 2 was described with the STO-3G basis including the 1s orbital for each atom, resulting in 4 spin-orbitals for the system and 4 qubits after Jordan-Wigner transformation.
H 2 in STO-3G basis is described by two spatial orbitals σ g and σ u with up and down spin functions.The reference states we used for the 3 excited states (T 1 , S 1 and S 2 ) are : The UCCSD circuit for this system of 4 spin-orbitals and 2 electrons is composed of 3 excitation operators that can be implemented in a compiled quantum circuit of depth 71 (with 44 CNOT gates).Computations were performed on Qiskit's QASM simulator acting like an ideal noiseless quantum computer, including finite sampling error.A shot scheduler between 1000 and 10000 shots was used during optimization (see section VI B 2),  and the final measurement was performed with 30000 shots.Figure 8(a) shows the results of the FS-VQE algorithm for excited states H 2 compared to the exact FCI states in the same basis in solid black lines.The FCI energies were obtained by numerically diagonalizing the electronic Hamiltonian matrix to obtain its eigenvalues.The ground-state results were obtained with standard VQE.FS-VQE allows recovering the complete potential energy curves for the 3 excited states of H 2 , at chemical accuracy.The absolute error is shown in the subplot of figure 8(a).

B. Excited states of LiH
LiH is a 4 electron system that can be described with 6 spin-orbitals in a minimal basis (considering s orbitals only for both atoms).Its excited energies were computed with FS-VQE using 6 qubits in Jordan-Wigner mapping.In this configuration, the UCCSD gate includes 8 excitations, resulting in a circuit of depth of 311 with 212 CNOT gates.The computations were per-formed on Qiskit's Statevector simulator allowing the measurement of the qubits' exact state, thereby avoiding finite sampling error.Note that the results LiH may differ greatly from the experimental data because the basis set only includes s orbitals, which is a poor approximation for the lithium atom.This minimal description allows us to put the FS-VQE to the test but does not aim for physically accurate results.
Here, LiH is described by three spatial orbitals σ 1 , σ 2 and σ 3 , each with up and down spin functions.The references we used for each excited state are as follows:  The mitigated FS-VQE algorithm was used to compute the highest excited state S 2 of the H 2 molecule in STO-3G basis, at a fixed bond length of 0.74Å, while adjusting the level of noise with the scaling factor λ as described in section V B. 20000 shots were taken for each circuit evaluation.Figure 9 compares the results of noisy FS-VQE simulations without mitigation with the corresponding results obtained with combined SPAM and ZNE mitigation methods.
The accuracy of the excited S 2 state energy computation is improved by employing mitigation techniques.When using mitigation, the energy computation reaches chemical accuracy compared to FCI for λ < 0.4, while only λ = 0 is below 1 kcal/mol without mitigation.The noise parameters for λ = 0.4 are given in Table I.
This may indicate that error mitigation techniques will be useful tools throughout the early fault-tolerant era of QPUs to extract more accurate data from noisy quantum devices.Our results suggest that chemical accuracy for small systems using our algorithm and the described mitigation methods could be reached for quantum devices with improved performances of about an order of magnitude compared to those of the present machines.In the long term, we expect noise mitigation methods to be one of the tools in the error correction arsenal for the FT era, but logical qubit-based error correction techniques will be inevitably required for large-scale quantum computing.
These experiments are a proof of concept that the FS-VQE algorithm can be combined with mitigation techniques to deal with noise in quantum computations.More detailed and extensive analysis of the use of mitigation techniques in FS-VQE is left for further studies.

A. Accuracy
In principle, the FS-VQE method with the UCCSD ansatz allows one to recover Coupled Cluster accuracy, provided that the optimization converges.For small systems such as H 2 or LiH with frozen core, CCSD is complete, and in theory we can recover FCI energies in the selected basis.However, for larger systems with more orbitals, single and double excitations are not generally sufficient to reach FCI accuracy and CC accuracy is expected when implementing UCCSD ansatz.Larger simulations or experiments would be needed for confirmation, as our work is restricted to very small systems.
In general, our simulations achieve good accuracy on simulated noiseless quantum computer, and all points in figure 8 have an error within the range of chemical accuracy (1 kcal / mol) compared to FCI, as evidenced by the subplots.However, a wide range of errors (between 1.5 mHa and 10 −16 mHa) and different patterns of error curves are observed.These differences can be attributed to the optimization process: the optimization terminates when a threshold value (10 −9 in our implementation) is reached in the cost function gradient, leading to various stages of convergence between different runs.
When noise models are included in the simulations, the addition of mitigation techniques is required to reach chemical accuracy.A simple implementation of SPAM mitigation and zero-noise extrapolation is sufficient to significantly improve the results of FS-VQE, as evidenced in figure 9.This result is promising for the next early fault-tolerant quantum era with lower error rates, where error mitigation is expected to play a major role and where FS-VQE could produce useful mitigated results.

B. Scaling and cost
One asset of the FS-VQE method is that the same ansatz circuit can be used for ground-state and excited state calculations.Excited states only require the evaluation of additional Pauli operators, meaning additional state preparation and measurement (SPAM) procedures using the same hardware requirements and the same quantum circuit structures (with possibly varying postrotation gates, representing minor changes).
The FS-VQE method has the disadvantage of involving the Hamiltonian square, making the number of Pauli string expectation values to evaluate larger than in standard VQE.After a Pauli grouping procedure (described in section IV B) the resulting number of required measurements is greatly reduced, as shown in figure 6.
Many questions remain open about the feasibility of the FS-VQE algorithm (and variational quantum algorithms in general) for large systems, as the classical part of the hybrid algorithm potentially retains intractable stages for large systems.Among these may be mentioned the number of measurements needed, the classical storage of measurement results, or the precomputation of the qubit Hamiltonian and of the FS operator.In particular, the Pauli reduction of the FS operator is a difficult classical task that would need further investigation to become scalable.The Pauli grouping procedure is also a crucial challenge, as it was shown to be a NP-hard problem [28], while the feasibility of the quantum FS method is highly dependent on it.In addition, optimizing variational parameters becomes more and more challenging as the system size increases, especially due to Barren plateaux.A careful design of the ansatz can overcome some of these challenges [46].
The scaling analysis of the proposed algorithm can be divided in terms of the number of circuit evaluations required, the circuit depth, and finally the number of shots to reach a target accuracy.The number of circuit evaluations corresponds to the number of Pauli groups, which scale as O(N 6 ) (discussed in section IV B) with respect to system size.The circuit depth depends on the ansatz.In this work, we implement the UCCSD ansatz that has a O(N 4 ) scaling in circuit depth, corresponding to the number of double excitations in the cluster operator.However, several ansätze have been proposed in the literature [47,48] that show a more favorable scaling, such as O(N ), while maintaining the same accuracy as UCCSD.Finally, the number of shots needed to achieve a certain accuracy ϵ scales as O( 1 ϵ 2 ) (see Appendix 1) for each Pauli string.Several techniques have been developed to address the measurement problem in VQAs, including Pauli grouping, measurement weighting, or shadow tomography [3].
The Folded Spectrum method is a general minimisation procedure that can be implemented within algorithms other than VQE to find excited states.Pauli operators are the building blocks of gate-based quantum computing, and we expect Pauli grouping procedures to remain relevant beyond variational algorithms, in which case Spectrum Folding could be one advantageous method, beyond variational algorithms, to compute molecular excited states on quantum computers.
Employing mitigation techniques seems to be central to obtaining meaningful results on noisy devices, but it comes at the cost of running deeper circuits, with more shots.Our implementation of ZNE requires four times more shots compared to the non-mitigated algorithm, for circuits of depth multiplied by γ=1,3,5,7.This overhead is non-negligible and needs to be addressed in the future.

IX. CONCLUSION AND PROSPECTS
In this work, we demonstrate that Folded Spectrum method is a successful approach for computing excited states using the VQE algorithm.The concept of evaluating the FS operator instead of the Hamiltonian to reach excited states is a well-known technique in QMC, and it could also be extended to algorithms other than VQE in quantum computing.Moreover, the Folded Spectrum procedure is agnostic to the quantum ansatz and fermion-to-qubit mapping scheme, and future improvements at any stage of the VQE algorithm can directly benefit this method.
Folded Spectrum allows one to directly compute any excited state around a target energy, which is a considerable asset compared to other methods where excited states are computed sequentially.This advantage is especially important for studying larger systems that have an increasing number of electronic states.It can be particularly useful for the computation of highly excited electronic energies and of great interest for the study of photo-chemical processes and light-matter interaction.
A major challenge to enable the scaling of variational methods is the preparation of good reference states for larger systems, which is particularly challenging for multi-reference states.As explained in Ref 1, a severe limitation to the scalability of ground-state VQE is the exponentially vanishing overlap of local reference states with the FCI ground-state.Our method also faces this issue, and further research is needed to allow for a systematic and scalable determination of reference states.However, we believe that finding excited states could be a less difficult task if the approach is to compute any excited state around a target energy rather than to search for a specific state (like it is the case for ground-state computations).In this scenario, our algorithm could in principle find a local minimum in the Folded Spectrum landscape that would correspond to an excited eigenstate.The increasing density of states in larger systems would benefit this approach and the reference state would thus play a less crucial role.
The main limitation of the FS method is the need to evaluate a squared Hamiltonian.In the quantum computing formalism of fermion-to-qubit mapping, this disadvantage can be alleviated by partitioning Pauli operators into commuting groups that can be evaluated simultaneously.The resulting number of evaluations needed to compute the FS operator expectation value is substantially reduced by the Pauli grouping.Despite this improvement, the number of shots required by the FS-VQE method for large systems is still prohibitive on quantum hardware, and further progress is needed to make the method scalable for practical applications.This is particularly central when dealing with noisy quantum processors, as error mitigation or error correction techniques are key to obtaining meaningful results, but come at the cost of additional quantum resources both in number of shots and number of qubits.
X. ACKNOWLEDGEMENT

Measurement precision
Measurement precision is directly related to the number of shots taken s.Let us consider the evaluation of ⟨ Ô⟩ with spectral decomposition : (1.1) The measurement of ⟨Ψ| Ô |Ψ⟩ relies on many repetitions of preparation of and measurement to evaluate the populations |⟨Φ i |Ψ⟩| 2 .Let V=⟨ Ô⟩ be the target value of the process.
⟨ Ô⟩ can be interpreted as the expected value of a random variable X having possible outcomes {o i } 2 N i=1 , with probabilities given by the Born rule : (1.2) Each shot of the experiment is a measure of X.By taking s shots we obtain a set of results X 1 ,...,X s .Thus, In this formalism, Chebyshev's inequality states that : with ϵ the precision of the result and σ 2 the variance of X, making σ 2 s−1 the variance estimate of the sample by means of the central limit theorem.
σ 2 can be bounded by a constant [49], so it can be deduced that in a worst-case scenario : Consequently, for a number of shots s, the precision on the expectation value ⟨Ψ| Ô |Ψ⟩ is of the order of 1 √ s .The extraction of classical information from a quantum system is therefore limited by a finite number of shots.This result is a direct consequence of the probabilistic nature of quantum mechanics.

Figure 3 :
Figure 3: Summary of the pre-processing steps and VQE algorithm to obtain a molecule's energy.

Figure 4 :
Figure 4: Illustration of Folded Spectrum method.Colored lines represent eigenstates along a vertically ascending axis.The eigenspectrum of Ĥ (left) is folded around ω in the spectrum of ( Ĥ − ω) 2 (right), causing its eigenvalues to reorder.The lowest eigenvalue of the folded spectrum, circled in red, is the excited eigenstate of Ĥ originally closest to ω.

Figure 5 :
Figure 5: Example of Pauli grouping for the 15 Pauli strings in H 2 Hamiltonian.Colors show the QWC partitioning that reduces the number of evaluations to 5, and red circles show the GC partitioning that only require 2 evaluations.

N 4 ∼ 3 ∼
O(N 4 ) terms into N −1 O(N 3 ) groups such that the operators within each set have disjoint indices, and therefore commute.In other words, instead of measuring each of the O(N 4 ) terms individually, one can perform O(N 3 ) measurements only to compute the expectation value of the asymptotically dominant Pauli strings in the molecular Hamiltonian.

Figure 7 :
Figure 7: Example of zero noise extrapolation curves coupled with SPAM mitigation on the expectation value of two Pauli commuting groups G 2 (top) and G 1(bottom).Noise scaling is performed by unitary folding with folding factors 1, 3, 5 and 7, and extrapolation is performed using a quadratic model.
(a) Potential Energy Curves of H2 (b) Potential Energy Curves of LiH

Figure 8 :
Figure 8: Results of FS-VQE calculations for H2 (a) and LiH (b) excited electronic states.Ground-states were obtained with a standard VQE.Colored markers show the (FS)-VQE results, and solid black lines are the FCI energies obtained by numerical diagonalisation of the Hamiltonian.The lower subplots represent the absolute error from FCI.The error plots were rescaled for clarity, and all non-visible points are below 10 −3 mHa for H2 and below 10 −8 mHa for LiH.(a) : Potential energy curves of H2 in STO-3G basis obtained on QASM simulator (ideal quantum computer with finite sampling shots error) using 30000 shots for the final evaluation.The error bars of finite sampling are represented but smaller than the markers.(b) : Potential energy curves of LiH in a minimal s-orbitals-only basis.LiH computations were performed on a statevector simulator, ignoring finite sampling error.

Figure 8 (
Figure 8(b) shows the results of FS-VQE for the potential energy curves of LiH.The solid black lines are the FCI states, obtained by diagonalization of the Hamiltonian.Absolute errors compared to FCI are presented in the subplot.

Figure 9 :
Figure 9: Comparison of the final energy accuracy obtained with non-mitigated and mitigated FS-VQE on the S 2 state of H 2 at 0.74Å.Horizontal axis is the simulated level of noise.