Symmetry enhanced variational quantum spin eigensolver

The variational quantum-classical algorithms are the most promising approach for achieving quantum advantage on near-term quantum simulators. Among these methods, the variational quantum eigensolver has attracted a lot of attention in recent years. While it is very effective for simulating the ground state of many-body systems, its generalization to excited states becomes very resource demanding. Here, we show that this issue can significantly be improved by exploiting the symmetries of the Hamiltonian. The improvement is even more effective for higher energy eigenstates. We introduce two methods for incorporating the symmetries. In the first approach, called hardware symmetry preserving, all the symmetries are included in the design of the circuit. In the second approach, the cost function is updated to include the symmetries. The hardware symmetry preserving approach indeed outperforms the second approach. However, integrating all symmetries in the design of the circuit could be extremely challenging. Therefore, we introduce hybrid symmetry preserving method in which symmetries are divided between the circuit and the classical cost function. This allows to harness the advantage of symmetries while preventing sophisticated circuit design.

The most popular application of variational quantum algorithms is to find the ground state of complex many-body systems. For fulfilling this task, the Variational Quantum Eigensolver (VQE) algorithm has been designed to target the ground state of a manybody system through minimizing the average energy [20,22]. The VQE has been extensively applied to quantum chemistry problems [16,[48][49][50] and experimentally realized on superconducting [16,49,55,56] and ion trap [1,10,57] quantum simulators. Several attempts have been made to enhance the VQE performance, including: minimizing the number of required measurements [58][59][60][61][62][63], improving the initialization [52,64,65], speeding up the classical optimization [66][67][68] and designing better circuits [69][70][71][72][73][74]. Several important phenomena in physics, such as topological phases [75], are described by the knowledge of a few low-energy eigenstates and not just the ground state. Therefore, the generalization of VQE to higher energy eigenstates is of high importance. However, the existing VQE protocols for simulating excited states [76][77][78][79] are very resource demanding which makes their scalability and practicality in doubt. Therefore, improving the performance of VQE in terms of feasibility for simulating large systems is a key milestone to achieve quantum advantage with NISQ devices.
Symmetry is one of the most profound concepts in physics, especially in quantum mechanics [80]. Most physical systems reveal various forms of symmetries for which a precise mathematical description has been developed [81]. Symmetries have been exploited to improve data mitigation [82], quantum machine learning [83], quantum state representation [84] and variational quantum optimization [32] in NISQ devices. In addition, the VQE algorithm can also hugely benefit from the incorporation of symmetries. There are two ways to incorporate symmetries in the VQE algorithms for simulating the ground state: (i) designing the circuit to naturally generate the quantum states with the relevant symmetry [52,[85][86][87][88][89][90][91]; and (ii) adding extra terms to the cost function to penalize the quantum states without the relevant symmetry [22,92]. Two key open problems still exist. Firstly, whether the symmetries can also be exploited for efficiently simulating the excited states. Secondly, which of the above methods, or a hybrid combination of them, are more effective for incorporating symmetries in the VQE algorithm for enhancing its performance.
In this paper, we show that symmetries can indeed significantly improve the VQE for simultaneously simulating several low-energy eigenstates. The improvement becomes even more pronounced for excited states. We introduce two different approaches for incorporating symmetries. First, in hardware symmetry preserving method, we include all the symmetries in the circuit. Second, we add symmetries as proper penalization terms to the cost function. Interestingly, our analysis shows that the first method is more effective with respect to both quantum and classical resources. However, designing a circuit which can integrate all the symmetries can be notoriously difficult. Hence, we introduce hybrid symmetry preserving method in which the two approaches are combined in order to harness the symmetries while keeping the circuit simple. Thanks to significant enhancement in resource efficiency, our proposal paves the way for achieving quantum advantage. In addition, it is very timely and can be implemented on existing quantum simulators.

Ground state VQE
In this section we provide a brief review on the VQE algorithm for preparing the ground state of a given Hamiltonian using a shallow quantum circuit [20]. In the VQE algorithm, a parameterized quantum circuit, represented by unitary operator U ( θ), is used to prepare a quantum state |ψ( θ) = U ( θ)|ψ 0 for a given N qubits. This parameterized quantum circuit is often referred to as the ansatz, with θ = (θ 1 , θ 2 , . . . , θ L ) being some tunable parameters in the circuit and |ψ 0 is the input state. By varying θ one can explore some part of the Hilbert space, which is called reachable set. In very deep circuits, and thus large number of parameters L, one may generate any possible quantum state of N qubits and thus the reachable set will be the entire Hilbert space. However, we would like to keep the circuit as shallow as possible and restrict ourselves to the most relevant part of the Hilbert space. In particular, in VQE algorithms we are interested in a fairly shallow ansatz for which the reachable set contains the ground state of the Hamiltonian of interest H. So far, several choices of the ansatz with different levels of complexity have been proposed [22,31,49,[93][94][95][96]. After choosing the ansatz, one can measure the average energy H = ψ( θ)|H|ψ( θ) through some appropriate measurements on the quantum device. This measured average energy will then be fed into a classical optimizer to be minimized through adaptively updating the parameters θ in the quantum circuit. Eventually, the optimization will be finished by obtaining the optimal parameters θ * . If the reachable set contains the ground state then the output of the circuit |ψ( θ * ) will be the ground state of H.
The conjecture behind the VQE is that a shallow circuit is enough to realize the ground state of the Hamiltonian. The price for keeping the circuit shallow, i.e. saving quantum resources, is to add a classical optimizer which demands extra classical resources. In this paper, we use L-BFGS algorithm as the classical optimizer [97]. If the optimization landscape suffers from Barren plateau [98,99] or the presence of significant number of local minima then the classical optimization may converge very slowly or even never reach the right quantum state. Therefore, in any VQE algorithm it is crucial to quantify both quantum and classical resources. Since single qubit operations are almost perfect, we use the number of two-qubit entangling gates (e.g. controlled-not) in our circuit as a quantification of quantum resources. On the other hand, for classical resources one has to notice that L parameters have to be optimized iteratively. Therefore, a logical definition for Classical Resources (C R ) can be the multiplication of the number of parameters L and the number of optimization iterations n I , namely The choice of ansatz is very crucial in all variational quantum algorithms. Perhaps the most widely used quantum circuit in the literature is the hardware efficient circuit [49] which is schematically shown in Fig. 1(a). In this circuit the single qubit rotations are defined as where σ α is the Pauli operator with α=x, y, z. In this paper, we use Controlled-Not (CNOT) gate as the two-qubit entangling gate in our circuits which is defined as Although, the hardware efficient circuit has been heavily used for solving various problems, it suffers from Barren plateau [98,99] which makes its training extremely difficult, in particular, when the number of layers increases. In addition, hardware efficient circuit does not conserve any symmetry. Therefore, as an alternative, one can use a more complex entangling gate with tunable parameters such as [100] N (θ x , θ y , θ z ) = e i(θxσ 1 where σ j x,y,z are the Pauli matrices acting on qubit j. The circuit of this unitary operation is depicted in Fig. 1(b). For the special case of θ x =θ y =θ z this unitary operator conserves the number of excitations as well as the total spin. By combining these two qubit gates and phase shift gates P (θ) = 1 0 0 e iθ , as shown in Fig. 1(c), one can make a quantum circuit U ( θ) which conserves the number of excitation, In the absence of phase shift gates, see Fig. 1(d), this quantum circuit can also preserve the total spin, namely [U ( θ), S 2 tot ] = 0, where S 2 tot = S 2 x + S 2 y + S 2 z . Note that, the convergence of the VQE would be successful only when the energy gap between the ground state and first excited state is larger than the standard deviation of the Hamiltonian with respect to state of the simulator, namely where ∆E = E 2 − E 1 is the energy gap between the ground state (with energy E 1 ) and the first excited state (with energy E 2 ). The situation becomes even more serious when there are degenerate eigenstates. In these cases, VQE algorithm may converge to an arbitrary superposition of the degenerate eigenvectors. The superposition may vary from one VQE implementation to another depending on the choice of the initial state or the circuit parameters. However, degeneracy among the eigenvectors is usually a consequence of the presence of some symmetries in the system. Indeed, degenerate eigenstates, despite having the same energy, have different eigenvalues with respect to the symmetry operators. Despite some attempts for exploiting the symmetry properties either through updating the cost function [22], or incorporating the symmetries into the design of the circuit [88], a systematic analysis in this issue is still missing. In addition, a few existing proposals [76,77], its generalization for higher energy excited states is very challenging. In the following sections, we try to address these two major issues through exploiting the inherent symmetries in the target Hamiltonian.

Excited state VQE
The most natural generalization of the VQE is to go beyond the ground state and target higher energy eigenstates. This can be very important in some topological systems in which the topologically relevant states are typically not the ground state [75]. In addition, in some physical phenomena, e.g. integer and fractional quantum Hall effects, the physics is fully described by low energy spectrum and not just the ground state. In order to realize higher energy eigenstates via VQE algorithms, two main methods have been developed: (i) penalty approach [77]; (ii) subspace-search VQE (SSVQE) algorithm [76]. In this section we briefly review these two methods and discuss their pros and cons.
In the penalty method, we assume that the first k − 1 lowest energy eigenstates are known through different VQE circuits. Then one can generate the k-th eigenstate of the Hamiltonian by penalizing the k − 1 lowest energy eigenstates in the cost function of VQE [77]. In order to target the k-th eigenstate of the system one can update the desired Hamiltonian where β i are some sufficiently large positive scalar for i-th eigenstate. By minimizing the cost function H , VQE yields the k-th eigenstate. This method provides a general and systematic excited state preparation algorithm. Nevertheless, the measurement of the projector terms requires calculating the overlaps | ψ( θ)|E i | 2 (for i = 1, · · · , k − 1), which are difficult to realize. So far, two methods have been proposed for computing the overlap between two quantum states. The first method requires an extended circuit in which the depth is increased by a factor of ∼ 2 [101]. The second approach utilizes the swap test [102,103] which demands doubling the number of qubits and requires complex many-body controlled gates. Considering the limitations in NISQ devices the penalty approach is unlikely to be beneficial in practice.
The weighted SSVQE provides an alternative method to generate all the k lowest energy eigenstates of a given Hamiltonian H [76]. In this algorithm, one uses a set of {|φ i } k i=1 orthogonal initial states (namely, φ i |φ j = δ ij ) as the input of a single parameterized quantum circuit, described by the unitary operator U ( θ). Since the initial states are orthogonal, the outputs U ( θ)|φ j , generated by the same circuit, are orthogonal too. In the weighted SSVQE algorithm, one has to minimize the cost function [76] where w 1 > w 2 > · · · > w k are real positive numbers. Minimizing the cost function in Eq. (7) produces all the k lowest energy eigenstates such that |E i = U ( θ * )|φ i . The major advantage of the weighted SSVQE is that it provides all the k lowest energy eigenstates in one single optimization procedure, without requiring any overlap of quantum states. However, the algorithm becomes more resource demanding as the number of target eigenstates increases.

Hardware symmetry preserving ground state simulation
In this section, we show how the exploitation of symmetries in the design of the circuit hardware can improve the VQE algorithm for simulating the ground state. We consider a 1-dimension chain of N spin-1/2 particle interacting via Heisenberg Hamiltonian where J = 1 sets the unit of energy and is the vector of Pauli operators at site i. The Heisenberg Hamiltonian in Eq. (8) represents a model of interacting system which cannot be mapped to free fermions and supports several symmetries including the conservation of the total spin, i.e. [H 0 , S 2 tot ] = 0, as well as its components in all directions, namely [H 0 , S α ] = 0. The first symmetry implies that every eigenstate |E k of the system has a specific total spin s which is an integer number for even N or half-integer for odd N , such that E k |S 2 tot |E k = s(s + 1). The second symmetry also guarantees that each eigenstate conserves the z component of the total spin such that E k |S z |E k = s z , with −s ≤ s z ≤ s. In particular, for even N , the ground state |E 1 is a global singlet with both s = 0 and s z = 0. Thus we represent it as |E S1 = |E 1 . The first excited state is a global triplet state with total spin s = 1 and triply degenerate with s z = 0, ±1.
For the sake of simplicity we use |E = |E 4 . The second excited state is another triplet state with s = 1 for which we similarly define |E The fourth eigenstate is a global singlet with s = 0 which we will represent it as We consider the VQE simulation of the ground state of H 0 using two different ansatzes, namely: (i) the hardware efficient circuit, shown in Fig. 1(a), which conserves no symmetry; and (ii) the S zconserving circuit, shown in Fig. 1(c), which conserves the number of excitation but do not preserve the total spin. Later, we will also consider circuits which can realize both S z and S tot symmetries. In both ansatzes, we initialize the circuit in the quantum state |0, 1, 0, · · · , 0, 1 . While in the hardware efficient circuit the choice of the initial state does not matter, in the case of S z -conserving circuit this choice is crucial and should have similar s z as the ground state |E S1 . To see the importance of exploiting the symmetry in the design of the circuit, we perform VQE for a system N = 16 qubits on both ansatzes. In Figs. 2(a)-(b) we plot the average energy and the obtainable fidelity as a function of optimization iteration, respectively. The error bars are computed for ∼ 100 repetitions of random samples. Interestingly, while only 5 layers of excitation conserving circuit is enough for a fast (∼ 300 iteration) convergence to F > 0.95, the hardware efficient circuit with the same depth can only reach F 0.4. In fact, the hardware efficient circuit can only reach fidelities above 0.95 when it contains at least 15 layers. Even for such a circuit, the optimization needs ∼ 2000 iterations to converge. This means that the hardware efficient circuit demands way more classical resources (C R = 870400 for F = 0.95) than the S z -conserving circuit (C R = 77500 for F = 0.95).
To better understand the difference between the performances of the two circuits one can determine the entangling power of the two circuits in terms of the number of CNOTs, as a quantum resource, as well as the number of parameters L, as a quantifier of classical resources (see Eq. (1)). To measure the entangling power, one can compute the average entanglement, quantified through von Neumann entropy, between the two halves of the system. For any quantum state |ψ( θ) at the output of circuit, one can compute the reduced density matrix of the left side of the system by tracing out the N/2 qubits on the right side, namely ρ L ( θ) = Tr R |ψ( θ) ψ( θ)| . The entanglement between the left and the right side of the system is then quantified by The average entangling power of the circuit is then defined as where the integration is performed over the whole pa-rameter space. For the numerical simulation, we approximate S av V by averaging S V ( θ) over 500 random samples of θ. In Fig. 2(c) we plot S av V versus the number of CNOT gates when the number of layers vary from 1 to 120. Note that, for the same number of CNOTs the hardware efficient circuit has many more layers than the S z -conserving circuit. Interestingly, the two ansatzes reach the fidelity F = 0.95 when their number of CNOTs are equal to 225. However, as the figure shows, for this number of CNOTs the hardware efficient have slightly more entangling power than the S z -conserving circuit. Both circuits can reach this fidelity only when their S av V is ∼ 3 times more than the entanglement in the real ground state. This is due to the fact that both circuits require similar number of CNOTs to reach the same fidelity. It is very insightful to consider the average entangling power S av V versus the number of parameters L as an indicator of required classical resources. In Fig. 2(d) we plot S av V as a function of L for the two circuits when the layers vary from 1 to 120. Remarkably, the S z -conserving circuit demands much less parameters for reaching the fidelity F = 0.95 as it demands a circuit with L = 155 (5 layers) in contrast to L = 512 (15 layers) for the hardware efficient circuit. This clearly shows that for the VQE simulation of the ground state while implementing the symmetry in the hardware may not reduce quantum resources, it significantly enhances classical resource efficiency.

Hybrid symmetry preserving excited state simulation
In this section, we focus on the effect of symmetry for the simulation of excited states using the weighted SSVQE. In particular, we consider generating the first 8 energy eigenstates of the Heisenberg Hamiltonian H 0 , namely |E S1,2 , |E . We first use the hardware-efficient ansatz, shown in Fig. 1(a), which supports no symmetry. Three different hardware efficient circuits are trained to target 1 (i.e. |E S1 ), 4 (i.e. |E S1 , |E ) different eigenstates. In Figs. 3(a)-(c) we present the average energy H 0 as a function of circuit layer for targeting 1, 4 and 8 eigenstates, respectively. For each case, the procedure is repeated for ∼ 100 random initial samples over which the results are averaged to be statistically meaningful. As the number of target eigenstates increases, the circuit needs more layers in order to converge for all the corresponding eigenstates. Importantly, as shown in Fig. 3(c), the eigenstates |E 5 to |E 8 do not converge properly even with 25 layers. In addition to the circuit layer, i.e. quantum resources, the classical optimization also gets more demanding as the number of target eigenstates increases. To quantify the required classical resources, in Fig. 3 the number of iterations n I needed for the convergence of the VQE circuit as a function of circuit layer for various number of target eigenstates. In order to see how the required resource varies as the number of target eigenstates increases, in Table. 1 we show the minimum required quantum resources and optimization resources for SSVQE simulation of k eigenstates using hardware efficient ansatz. The number of circuit layers and the number of optimization iterations are chosen such that the energy precision for all the k eigenvectors becomes less than 0.5J, namely

(d) we plot
These results clearly show that targeting more eigenstates makes optimization slower and demands more iterations. Thus, one may wonder whether the exploitation of symmetries can help to make the SSVQE more resource efficient.
Considering one symmetry is not enough to truly distinguish different eigenstates. For instance, |E S1,2 and |E (0) T1,2 all have s z = 0. In order to further discriminate the global singlets from the global triplets, one has to also take the total spin s into account. In this section, we use a hybrid approach for implementing the symmetries. Namely, we use the S z -conserving circuit to target quantum states with specific s z and update the cost function to target quantum states with the right total spin s. By using the S z -conserving circuit in the SSVQE algorithm one can only target the eigenstates with a given s z . For instance, four eigenstates, namely |E S1,2 and |E ) can be targeted on the same S z -conserving circuit with initial states satisfying s z = +1 (s z = −1). This division of eigenstates between different quantum circuits can significantly reduce the circuit complexity in terms of the required number of CNOTs.
In order to even further simplify the SSVQE circuit one can exploit the total spin symmetry as well. Implementing total spin symmetry in the hardware requires extra CNOT gates, which will be discussed later. Here, we take a hybrid approach and include the total spin symmetry in the cost function. In particular, for the preparation of |E S1 and |E S2 via SSVQE, one can adopt the S z -conserving circuit with two orthogonal initial states |φ 1 = |0, 1, 0, · · · , 1 and |φ 2 = |1, 0, 1 · · · , 0 and the following cost function where β = 1000 is a positive constant which is taken to be sufficiently large as in this type of penalty terms the final error in estimation of the cost function is of the order O(1/β) [104]. In Figs. 4(a)-(b), we plot the average energy H 0 and fidelity as a function of optimization iteration, respectively, in a circuit of system size N = 14 with 18 layers. The error bars in the figure are computed through averaging over 100 random initial samples. Remarkably, both eigenstates can be generated with fidelity above 0.95 after only ∼ 500 iterations. This is much better performance than the hardware efficient circuit, see Fig. 3(c), in which the output did not even converge properly to |E S2 for a smaller system of size N = 10 with even 25 layers and ∼ 5000 iterations.
where, again, β is a sufficiently large positive scalar. This type of cost functions are different from the one used for the singlets, see Eq. (10), as the operator in the penalty term requires four point correlation functions. In Ref. [104], it is shown that β for these type of cost function can be chosen to be smaller, lower bounded by the energy gap. In fact, for our triplet cost function in Eq. (11), choosing β = 2 is enough for convergence. In Figs. 5(a)-(b), we plot the average energy H 0 and the obtainable fidelity as a function of optimization iteration, respectively, for a circuit of system size N = 14 with 20 layers. Similar to the previous case, the results are repeated for 100 random samples over which the error bars are estimated.
In the numerical simulation, we set β = 2. Indeed, the SSVQE successfully generates |E T2 , simultaneously, achieving the fidelity F > 0.95. It is worth emphasizing that in comparison with the preparation of |E S1 and |E S2 , targeting the triplet eigenstates is more difficult in terms of required circuit layers. However, thanks to a stronger penalization term in the cost function, the error bars for triplets are smaller than their singlet counterparts, in particular at initial iterations, see Fig. 4. Our procedure shows significant improvement over the results without symmetry, shown in Fig. 3(c), in which |E T2 fails to be generated for a fairly small system size of N = 10 with even 25 layers of hardware-efficient ansatz and ∼ 5000 iterations. Similarly, one can target |E ) using the same circuit with proper choice of initial states with s z = +1 (s z = −1). For the sake of brevity, we do not present the results of these simulations.

Hybrid versus Hardware symmetry preserving
In the previous section, we showed how hybrid symmetry preserving method can enhance the performance of VQE through a combination of S zconserving circuit and updated cost function. In this section, we introduce hardware symmetry preserving approach in which both of the symmetries are included in the design of the circuit and thus the cost function only minimizes the average energy. In Fig. 1(d), we present one layer of a quantum circuit which conserves the S tot . This ansatz is very similar to S z -conserving circuit except phase gates which are removed. As mentioned before, this circuit not only conserves the z component of the spin, namely s z , but also preserves the total spin s as well. Thus, by choosing a proper initial state with a specific total spin s and a given spin component s z one can guarantee the preservation of these symmetries in the output state |ψ( θ) . For the case of global singlets, i.e. s = 0, one simple initial state is |ψ 0 = ⊗ N/2 |ψ − , where |ψ − = (|01 − |10 )/ √ 2 is a two qubit singlet state. To generate this initial state, one has to use at least N/2 extra CNOTs at the beginning of the quantum circuit. In the SSVQE, if one wants to target two global singlets, i.e. |E S1 and |E S2 simultaneously, one has to generate another global singlet initial state which demands extra CNOTs. In fact, by increasing the number of target eigenstates, generating proper initial states which are all orthogonal and satisfy s = 0 become more complex and demand extra CNOT gates. In the case of global triplet with s = 1, the situation is simpler. For instance, by taking the quantum state ⊗ N/2 |ψ − and locally converting one of the two-qubit singlets |ψ − into a two-qubit triplet (with a desire s z ), one can generate N/2 different orthogonal global triplet initial states. For higher total spins (i.e. s > 1), the simple converting of two-qubit singlets into triplets does not create an initial state with a given s. Therefore, S tot -conserving symmetry circuits become more complex, see Ref. [88] for more detailed discussion.
To better understand the impact of symmetry on VQE simulation, we compare the performance of the hybrid symmetry preserving approach, which uses S zconserving circuit, and the hardware symmetry preserving approach, which utilizes S tot -conserving circuits. As an example, we focus on generating the ground and the first excited state of the Heisenberg Hamiltonian, namely |E S1 (with s = 0 and s z = 0) and |E (0) T1 (with s = 1 and s z = 0), with the size of N = 14. In Fig. 6(a), we plot the obtainable fidelity for targeting the ground state |E S1 as a function of circuit layers for the two circuits. As the figure shows, the hardware symmetry preserving approach reaches the fidelity F > 0.95 with less layers. Despite requiring extra CNOTs for preparing its initial global singlet state |ψ 0 = ⊗ N/2 |ψ − , the hardware symmetry preserving method is still more efficient in terms of two-qubit entangling gates, using 117 CNOTs versus 156 CNOTs in hybrid preserving symmetry method. In Fig. 6(b), we plot the classical resources, defined in Eq. (1), as a function of layers which shows significant advantage for the hardware symmetry preserving method. The superiority of the hardware symmetry preserving method becomes more evident when one targets the first excited state |E (0) T1 . In Fig. 6(c), we plot the fidelity as a function of layers for both quantum circuits. The hardware symmetry preserving method shows significant advantage as it reaches the fidelity F > 0.95 only after 4 layers (with 156 CNOTs) in comparison with 7 layers of hybrid symmetry preserving method (with 273 CNOTs). In Fig. 6(d), the corresponding classical resources are compared. Interestingly, the hardware symmetry preserving method not only demands much less classical resources, but also benefits from less fluctuations.
For the sake of completeness, we also compare the performance of hybrid and hardware symmetry preserving methods using SSVQE for targeting |E S1 and |E S2 as well as |E  Table. 2. As the results show, hardware symmetry preserving method outperforms the hybrid symmetry preserving approach.
This analysis shows that the hardware symmetry preserving approach is the most efficient way for exploiting symmetries in SSVQE in terms of both quantum and classical resources. However, the price that one has to pay is the more complex circuit design. In particular, for total spin s > 1 the initialization may indeed need complex quantum circuits [88].

Generality of symmetry method
So far, we have focused on Heisenberg Hamiltonian for which we exploit S z and S tot symmetry for preparing several low-energy eigenstates. Here, we show that symmetry method can be generalized to other Hamiltonians with different symmetries. In particular, we will show that exploiting symmetries can also enhance the VQE performance for simulating free fermionic systems. For example, we consider the Ising Hamiltonian with transverse field where J z is the exchange coupling and h x is the strength of the transverse magnetic field. It is well known that this Hamiltonian has a quantum phase transition at J z /h x = 1. At the critical point where the quantum phase transition takes place, the ground state is highly entangled and has a complex form. Indeed, it has been shown that at the critical point the convergence of VQE requires more quantum and classical resources [105]. In other words, one needs a deeper quantum circuit as well as more iterations to reach the target state. Therefore, we focus at the critical point, as the most complex point in the phase Since the symmetry operator S x is a global action, designing a quantum circuit that implements this symmetry is very difficult. Therefore, hardware symmetry preserving approach becomes infeasible. Consequently, we use updated cost function for generating the first two eigenstates of the Hamiltonian. We use the circuit shown in Fig. 7(a). For starting the state preparation, we initialize the circuit in the quantum state ⊗ N |+ , where |+ = (|0 + |1 )/ √ 2. While for generating |E 1 a simple cost function H I is enough, for preparation of |E 2 one can simply update the cost function as cost = H I + β( S x + 1) 2 (13) where β is a positive scalar to make the two terms of the same order (here, we put β = 1). In Figs. 7(b)-(c), we plot the average energy H I and the average fidelity over 50 random samples as a function of optimization iterations, respectively, for a system of size N = 8 with 4 layers of ansatz. As shown in the plot, as optimization iteration increases, the VQE successfully generate |E 2 with fidelity F > 0.99. For generating more eigenstates one has to combine the symmetries with SSVQE algorithm as we did above for the Heisenberg Hamiltonian.

Conclusion
The VQE has emerged as one of the leading NISQ algorithms with the potential of achieving quantum advantage. So far, it has been successfully applied for simulating the ground state of condensed matter systems and several chemical materials. Nonetheless, generalization of the VQE to higher energy eigenstates is very resource consuming, demanding deep circuits and lengthy classical optimizations. In this paper, we address this problem by exploiting symmetries to enhance both quantum and classical resource efficiencies of the VQE algorithm. The symmetry enhanced resource efficiency becomes even more effective when higher energy eigenstates are targeted. Indeed, some of the excited states cannot be reached without including symmetries in the VQE algorithm. We have considered two ways for incorporating symmetries. In the first approach, which we call it hardware symmetry preserving method, all the symmetries are included in the circuit. In the second method, the symmetries are integrated in the cost function. Our results show that the hardware symmetry preserving method significantly outperforms the penalization of the cost function. However, implementing all symmetries in the design of the circuit may not be practical. Therefore, we have introduced the hybrid symmetry preserving method in which some of the symmetries are included in the circuit and the rest are incorporated in the cost function. This allows to simultaneously improve the resource efficiency and keep the circuit design simple. Our proposal achieves significant resource efficiency and thus paves the way for achieving quantum advantage on NISQ simulators. In addition, it is very timely and can be adopted in existing quantum simulators. Note that our method is applicable only we have prior knowledge about the symmetries of the system. For any given Hamiltonian, in order to investigate the symmetries of the system, one can use a small system as symmetries are not length dependent. In addition, usually the symmetry values of the low energy eigenstates can also be determined for small sizes and directly be used for large systems. In other words, the fact that symmetries are not length dependent can be very beneficial in specifying the symmetries and incorporating them in VQE algorithms.

Data availability
All the source codes for generating the data is available [106]. In addition, the data which have been for generating plots can be provided upon reasonable request from the authors.