Sequence of penalties method to study excited states using VQE

We propose an extension of the Variational Quantum Eigensolver (VQE) that leads to more accurate energy estimations and can be used to study excited states. The method is based on the introduction of a sequence of increasing penalties in the cost function. This approach does not require circuit modifications and thus can be applied with no additional depth cost. Through numerical simulations, we show that we are able to produce variational states with desired physical properties, such as total spin and charge. We assess its performance both on classical simulators and on currently available quantum devices, calculating the potential energy curves of small molecular systems in different physical configurations. Finally, we compare our method to the original VQE and to another extension, obtaining a better agreement with exact simulations for both energy and targeted physical quantities.


I. INTRODUCTION
In recent years, the development of quantum computing hardware made strides from the first working prototypes to devices with more than a hundred qubits [1]. Companies such as IBM and Google have revealed plans to build thousands or even a million qubits devices by the end of the decade. However, it has to be highlighted that the number of qubits is not the only metric to take into account in a quantum computing platform and the current amount of operations that can be performed is strongly limited by hardware noise and decoherence. While quantum error correction schemes have been proposed and even tested on hardware [2][3][4][5][6], their large-scale application remains still far in the future.
In order to provide meaningful results in the near-term, hybrid algorithms that combine quantum devices with classical computers have been proposed [7][8][9][10][11][12][13][14][15][16][17]. In these schemes a quantum computer is only in charge of a subroutine, while the classical computer governs the whole algorithm. This greatly reduces the number of quantum operations required and those methods have been demonstrated to be naturally robust with respect to certain hardware errors [18].
One such scheme is the Variational Quantum Eigensolver (VQE) [19,20], where a parameterized quantum circuit is used to approximate the ground state of an interacting quantum system. In this case a quantum computer is used to measure the energy and its derivatives, then a classical optimizer will tune the variational parameters according to those. VQE has been widely studied and extended to improve the final result precision or to lower the hardware requirements [21][22][23][24].
More recently, extensions of VQE to the study of excited states have been proposed. As an example, in the quantum subspace expansion method [25], excitation op- * r.carobene@campus.unimib.it erators are applied to the variational ground-state, while the quantum imaginary time evolution algorithm [10] can be used to construct the Lanczos subspace [26]. Alternatively, one can directly compute excitation energies using the quantum equation of motion method [11], or even use algorithms of quantum machine learning [27]. Finally, another category of VQE variants is based on the idea of modifying the cost function to guide the optimization to the desired state [28][29][30][31][32][33].
In this paper, we expand on the last approach and propose a method based on the introduction of a sequence of increasing penalties in the cost function that leads to better approximated eigenstates and can be used to study excited states. We chose to name our algorithm Sequence of Penalties VQE (SPVQE). Based on a cost function modification, the proposed method does not need more resources than VQE.
The structure of this paper is as follows: in Section II we review the VQE algorithm, its constrained modification (CQVE) and present the Sequence of Penalties VQE method, while in Section III we apply it to the study of excited and ionized state of small chemical compounds, assessing the performances both on classical simulators and on real hardware. Finally, Section IV concludes the paper with some considerations and outlooks on the proposed method.

A. VQE: Variational Quantum Eigensolver
Here we briefly review the Variational Quantum Eigensolver (VQE), for a more detailed explanation we suggest to refer to [18,19]. Consider a physical system described by the Hamiltonian H. The aim of the VQE is to prepare a variational ground state approximation of this system on a quantum computer. First, the Hamiltonian is mapped into a Pauli operator [34][35][36] Figure 1. Sketch of the cost function space of H + 3 projected on one and two dimensions. A comparison of the minima before and after applying a penalty constrain on the total spin shows that the desired state becomes the global minimum after the transformation.
tion values on quantum hardware. Then, the qubits are prepared in the state |ψ is a combination of parameterized quantum operations (gates) depending on the set of parameters θ, and θ 0 is the initial choice of parameters.
We use the quantum device to measure the expectation value of H, defined as E(θ) = ψ(θ)|Ĥ |ψ(θ) , and its derivatives with respect to the parameters θ [37][38][39][40]. Finally, we feed these quantities to a classical optimizer that determines new values of θ in order to decrease E(θ). When this iterated procedure converges, we obtain a set of parameters θ that defines the desired state. Now suppose that, given some operatorÂ, we want to prepare the eigenstate corresponding to its eigenvalue a. Given that these configurations may not correspond to the ground state of the system, we can not rely solely on lowering the energy.
More generally, we may be interested in studying an excited state in a region in which quantum states with different physical properties have similar energies. In all these cases, VQE may fail to correctly approximate the desired state.
In perspective of using VQE to study molecules with complex energy spectra, we must develop a robust method that preserves every physical property we may desire to fix.

B. CVQE: Constrained Variational Quantum Eigensolver
Some of the proposed variations of VQE already include in the optimization process information about additional operators [41][42][43]. In particular, we will briefly recall the Constrained Variational Quantum Eigensolver [29].
While in the standard VQE approach we minimize the energy as stated in Section II A, the CVQE method introduces a redefined cost function F (θ) by adding a penalty multiplier for each of the operators we may desire to fix: hereÂ i with i ∈ {1, . . . , m} is a set of operators, µ i are their corresponding multipliers and a i the eigenvalues that we want to fix. Chosen µ i sufficiently large, the desired state will correspond to a global minimum of the cost function. In Figure 1 we show as an example the variational landscape computed for the H + 3 molecule, obtained by scanning the Hamiltonian and the spin operator on two random parameters. We see that the desired state is transformed from a local to a global minimum through the application of the constrain. We highlight that, among all the possible states satisfying the given constraints, the global minima will be always the lowest energy state. For this reason, we envision the combination of our method with other algorithms, such as the Quantum Subspace Expansion [25] or the quantum Equation of Motion [11], to provide even more accurate estimations of excited state properties of physical systems.
In a practical usage of CVQE, however, the combination of an unknown good starting point θ 0 and restricted local knowledge of the loss landscape often leads to convergence into local minima. Moreover, a careful fine tuning of µ i is often required since multipliers too large yield to narrow global minima, while multipliers too small do not highlight sufficiently the global minima with respect to local ones.

C. SPVQE: Sequence of Penalties Variational Quantum Eigensolver
We propose a simple, yet effective, variation of CVQE: the Sequence of Penalty VQE (SPVQE).
We want to have a penalty large enough to exclude any local minima, while avoiding changing the cost function space too rapidly. Therefore, we propose to repeat the constrained optimization increasing penalty multipliers by steps [44].
For clarity of exposition, we will constrain a single operator. The method can be extended to multiple constraints modifying the penalty term.
First, we choose the maximum value of the penalty multiplier (µ max ) and the number of steps N s . If the gap between the desired state energy E and the ground state one (E GS ) is approximately known, we can choose the value of µ max so that it satisfies the inequality [45]: It is possible to obtain an estimation of the gap by using classical approximated methods. Even if this classical estimation does not reach chemical accuracy, it provides an idea of the magnitude of the gap. In general, N s and µ max are hyperparameters of the SPVQE method, therefore there is not a rigorous framework to exactly set the parameters and we have to rely on heuristics. Equation (2) provides a lower limit for the maximum multiplier, but we do not have an upper limit for the single step multiplier and therefore there is no guarantee that the choice made will leada to the feasible region. In a theoretical perspective, µ max can be chosen arbitrarily large since N s compensates for the narrowing of the global minimum, hence the problem is reduced to the choice of N s .
In practice, we can first set N s by considering how much time we can afford on the quantum platform, given that computational time grows linearly with N s . In our examples a number of steps of ≈ 10 was always sufficient, as we will show in Figure 4.
We run N s instances of VQE constrained with increasing penalties. At each iteration k ∈ {1, . . . , N s }, we compute the penalty using µ = µmax * k Ns , using the optimal set of parameters θ k−1 of the previous step as starting point. Finally, the best loss computed over all the iterations is returned as the result.
The complete method is schematically presented in Algorithm 1.

D. Numerical simulations
To demonstrate a practical SPVQE application, we considered molecular systems where the constrained operators were the particle number or the total spin operator:Â We chose a hardware efficient ansatz composed by a layer of single qubit rotations R y , each one with an independent parameter, followed by an entangling layer of CNOTs. This structure is repeated D times and, at the end of the circuit, another layer of R y rotations is applied. As an example, for 4 qubits the circuit reads In the following experiments we set D = 3, unless where explicitly stated. To limit the number of qubits required, we considered a minimal Slater-type orbital basis set constructed with 6 primitive Gaussian orbitals (STO-6G) [20,23]. The qubit mapping of the Hamiltonian is done with parity mapping [35] with the two qubit reduction. For noiseless simulations we used the SciPy Conjugate Gradient (CG) optimizer [46]; whereas for noisy simulations and hardware calculations we used the Nakanishi-Fujii-Todo (NFT) optimizer [47] that proved to give the best results. Some final measurements were performed after the VQE procedure to improve results in noisy simulations and hardware calculations. These are presented in B. The simulations and quantum hardware computations that can be found in Section III and in the appendices are performed using IBM's open source library Qiskit [48].

III. RESULTS
In this Section we will use the SPVQE algorithm to find ground and excited states of molecular systems and compare the results to VQE and CVQE calculations. In particular, we will consider the Born-Oppenheimer approximation of molecular Hamiltonians, that in second quantization have the form whereâ † i andâ i are fermionic creation and annihilation operators corresponding each to a different spin-orbital and the coefficients h ij and g ijkl are coefficient accounting for electron kinetic energy, interactions electronsnuclei and electron-electron repulsion. E 0 is the shift obtained by fixing the atomic positions and calculating their repulsion energy.
Details on the systems studied can be found in A, while complete results including nuclear repulsion energy, frozen core energy and parameters obtained at every step are available on [49].
A. Using SPVQE to study excited states Calculations were performed on a classical computer simulating a noiseless quantum environment. The circuit describing H 2 and its excited states needs two qubits, for a total of 8 parameters when a D = 3 circuit is considered.
The dissociation profiles computed with SPVQE are presented in Figure 2 and compared to the numerically exact solutions obtained with the classical Full Configuration Interaction (FCI) method [50].
In all simulations presented in Figure 2 SPVQE is able to target the physically correct state among the various ones that are present.
Then, we repeated these simulations using CVQE and VQE and show a comparison of the different methods in Figure 3. VQE points correspond to the best value obtained among 1000 simulations with random starting parameters. CVQE simulations were repeated 100 times and SPVQE calculations 2 times.
The standard VQE algorithm struggles to find a good approximation for every configuration. In particular, it it needs a higher penalty (in particular at low bond lengths) and this misleads the algorithm that fails to find the correct minimum. Lastly, SPVQE reaches the correct state at every configuration, hitting machine precision.

B. Dependence on the number of steps
In this section we analyze the performance of SPVQE as a function of the number of sequence steps N s . In Figure 4, we show the results for three different molecules: Na − , H 2 O and H + 3 . While H 2 O is not an excited state, the addition of a constraint is still needed because the molecule has an excited level with a very similar energy, making the convergence of standard VQE to a state with the right physical properties difficult [23,29].
For H 2 O we varied the length of one of the H − O bond from 1 Å to 3 Å, while for H + 3 we varied the length of all H − H bonds from 1 Å to 3 Å. These configurations are difficult to handle for the standard VQE, which is not able to reach a correct approximation.
For Na − and H 2 O the simulations were carried in the active space, freezing core orbitals to diminish the number of needed qubits. This resulted in an ansatz with 24 parameters on 6 qubits for Na − ,considering 2 electrons in the 4 most external molecular orbitals, and one with 16 parameters on 4 qubits for H 2 O, considering just 4 electrons in 3 orbitals. H + 3 simulations are performed without frozen core approximation, requiring 4 qubits and 16 parameters. Thus, every marker on the graphs represents the average result obtained considering a SPVQE iteration for each ionic configuration.
The plots show that even with few iterations of the SPVQE method, the total spin is correctly constrained to the desired value, while the energy decreases accordingly. We note that the first marker of each graph, namely the single-step optimization,is equal to a CVQE calculation.

C. Robustness against starting point choice
Finding the ground state of a quantum system is a hard problem even on quantum computers [51]. For this reason, VQE and its modifications are not guaranteed to reach the global minimum of the energy (or the modified cost function) in polynomial time in every scenario. The  . Results for Na − , H2O and H + 3 with µmax = 1 constraint on the total spin operator. Every point corresponds to a full SPVQE calculation with a different number of steps. Increasing the number of sequence steps improves the accuracy of our variational approximation, both for the energy and the total spin, which is the constrained operator.
convergence of the algorithm depends on different factors. The choice of starting parameters is one of them and can hinder the convergence of the algorithm. Despite this importance, often we do not have the possibility to make an educated guess on the starting set of parameters; thus a certain robustness with respect to the choice of the starting parameters is desired.
In this section we analyze the robustness of SPVQE and compare it to CVQE. To this aim, we performed multiple simulations of the H + 3 molecule for bond lengths varying from 0.4 Å to 2.5 Å. Results of these calculations are shown in Figure 5. Every marker is obtained considering multiple simulations with random starting parameters and represents the mean error of the computed energies, numbers of particle, spins. The coloured areas of the violin-plot mirror the underlying data distribution.   Figure 5. Energy, total spin and particle number error comparison between the variational approximation of the ground state of H + 3 at different bond length obtained with CVQE and SPVQE. The error is calculated as |Ox-VQE − OFCI|, where O is the observable considered. Every marker corresponds to the mean error of 100 different simulations with random starting parameters. Colored areas mirror the underlying data distribution. SPVQE proves to be more robust in every analyzed situation, with lower errors and smaller deviations from the mean value. This is due to the fact that SPVQE is generating itself a sequence of improving educated guesses.
To show that this effect is not due to some inherent property of the H + 3 molecule, in C we report the calculations performed on different chemical systems.

D. Hardware experiments
Finally, we tested SPVQE on real quantum hardware, using the IBM Quantum processor ibmq_quito.
We analyzed the trihydrogen cation (H + 3 ) for bond lengths varying from 0.3 Å to 2.5 Å. A standard VQE approach works properly for bond lengths smaller than ≈ 1.5 Å: after that point, the energy of H − 3 gets lower than the desired configuration, making the algorithm converge to the wrong minimum. Thus, the optimization process fails to preserve the particles number (it finds a state where N = 4) and also the spin (it computes S 2 = 1). Therefore we chose to constrain S 2 = 0.
As we stated before, current quantum hardware has limited performances and requires extreme care to reduce errors impact. Since ibmq_quito has a medium-low quantum volume [52] of 16, we reduced the depth of the ansatz described in Section II D choosing D = 2. Once the circuit is transpiled using ibmq_quito basis gates {CX, ID, RZ, SX, X}, its final depth is 19. For every bond length configuration we performed 4 independent calculations with 1024 shots each. Starting parameters have been chosen randomly.
Hardware experiments were prepared simulating the SPVQE algorithm in a noisy environment. The obtained results, reported in D, confirmed us that SPVQE could mantain the desired advantages even in the presence of noise.
Hardware results are showed in Figure 6.
We can see that the profile computed using SPVQE correctly approximate the exact calculated one. VQE fails to correctly approximate the state of H + 3 at high bond distance, when H − 3 energy level became lower than the H + 3 one. Error bars correspond to the standard deviation obtained among the 4 different results and thus represent the statistical error.
VQE shows the biggest statistical uncertainty at at 2 Å. For this configuration, the VQE wave functions always present the wrong spin. Therefore, this point can be interpreted as corresponding to a molecular configuration where it is difficult to find the global minimum of the Hamiltonian.
Looking at the lower panel of Figure 6, we see that SPVQE is able to constrain S 2 to the desired value even on hardware, while VQE starts to miscalculate it at 1. fact that VQE just minimizes the energy, thus finding a state not characterized by a precise spin, but rather by a mixed spin state.

IV. DISCUSSION
In this paper we introduced a simple, yet effective, modification of the VQE algorithm to enable excited and ionized state simulations. The proposed method, Sequence of Penalties VQE, was tested for small molecules both on classical simulators and on real quantum hardware and proved to work correctly in the studied cases. Moreover, the algorithm showed an increased robustness against the choice of the starting point. For these reasons, we envision to use SPVQE in combination with other methods to study excited states of physical systems with higher accuracy on quantum computers.
The choice of the parameters, number of iterations and maximum penalty is done heuristically and can be further optimized in future studies. Then, another outlook is the integration with error mitigation techniques [53][54][55][56], in order to improve the accuracy of the obtained results and scale the method to bigger physical systems.
In the end, we think that SPVQE should be taken in consideration for calculations regarding excited states, considering the straightforwardness of its implementation, its reliability and the accuracy improvements that it brought. To reduce computational burden we restricted the active space with the Qiskit ActiveSpaceTransformer. All the computed values, including frozen core energies and nuclear repulsion energies, are available on GitHub [49].
In Table I we summarize some of the most important simulation details for every studied molecule. Here, bold symbols correspond to computations on hardware. To improve both CVQE and SPVQE results in noisy and real quantum environments, we compute the expectation value of the Hamiltonian one last time after the algorithm has converged. In fact, even if we have the value of the cost function F and the penalty P correspondent to the best result found, the introduction of a penalty leads to a larger uncertainty. We can show that with error propagation theory: Measuring a last time the expectation value of the Hamiltonian, using the best computed parameters, give us a measure of the best energy with the error associated minimized, because we have: Appendix C: Robustness of SPVQE for different molecular systems While in Section III C we presented the result for H + 3 , we now show the same analysis for different molecules. This analysis shows that the robustness of SPVQE is not limited to the H + 3 system. We simulated, in a noiseless environment, different molecules applying SPVQE and CVQE multiple times with different and random starting parameters.
In particular, we studied the H 2 O molecule for different bond lengths (we varied the position of one hydrogen atom in respect to H − O. To simulate a molecule with 10 electrons, we froze the core orbitals, leaving just 4 particles to be placed in 3 molecular orbitals.
We highlight that results both of CVQE and SPVQE were obtained setting the same number of iteration for the optimizer.
We also repeated the same analysis for the N 2 molecule where we froze the core orbitals considering only 6 electrons in 3 molecular orbitals.
Results are shown in Figures 7 and 8.

Appendix D: Noise simulations of SPVQE
We simulated H + 3 dissociation profile considering a noise model imported from the ibm_perth IBM Quantum processor. Trihydrogen cation has 6 spin orbitals on the lowest shell: we therefore need a minimum of four qubits and a 16 parameters ansatz. For every method, to help the optimization, we used as starting ansatz parameters the optimal ones calculated for the last computation of that same method.
The results are presented in Figure 9.   We can see that SPVQE outperforms both the standard VQE and the CVQE approach. The VQE fails to converge to the right physical state for bond lengths greater than 1.4Å.
Compared to CVQE, the SPVQE approach converges to the desired value with higher precision both for total spin (the constrained operator) and energy.