Accelerated variational algorithms for digital quantum simulation of many-body ground states

Digital quantum simulators can provide universal quantum computation outperforming all existing classical computers. One of they key applications of quantum simulators is to emulate the ground state of many-body systems, as it is of great interest in various fields from condensed matter physics to material science. Traditionally, in an analog sense, adiabatic evolution has been proposed to slowly evolve a simple Hamiltonian, initialized in its ground state, to the Hamiltonian of interest such that the final state becomes the desired ground state. Recently, variational methods have also been proposed and realized in quantum simulators for emulating the ground state of many-body systems. Here, we first provide a quantitative comparison between the adiabatic and variational methods with respect to required quantum resources on digital quantum simulators, namely the depth of the circuit and the number of two-qubit quantum gates. Our results show that the variational methods are less demanding with respect to these resources. However, they need to be hybridized with a classical optimization which can converge slowly. Therefore, as the second result of the paper, we provide two different approaches for speeding the convergence of the classical optimizer by taking a good initial guess for the parameters of the variational circuit. These approaches show significant improvement in the optimization procedure and thus allow for more efficient simulation of the variational methods on fairly shallow circuits.


I. INTRODUCTION
Simulating strongly correlated many-body systems at and out of equilibrium is one of the key tasks in condensed matter physics. On classical computers, simulating a general quantum many-body system is inherently intractable, due to the exponential growth of the Hilbert space with respect to the system size. The only genuine approach for understanding a quantum system is to emulate its behavior on a quantum simulator, which is a quantum device with more controllability and versatility [1][2][3][4]. Quantum simulators are now rapidly emerging in various physical systems, including optical lattices [5][6][7][8][9], Rydberg atoms [10], ion traps [11][12][13][14][15], photonic circuits [16,17], quantum dot arrays [18], dopants in silicon [19,20] and superconducting devices [21][22][23][24][25][26][27][28][29][30][31][32]. There are two different approaches for implementing quantum simulation, namely analog and digital. In the analog approach, the particles interact via some Hamiltonian, possibly time-dependent, and the quantum simulation is the continuous evolution of this system [7-9, 25, 26, 30]. Although, reasonably strong noise can be tolerated [33][34][35], the analog quantum simulation is not universal, in the sense that only limited types of unitary operations can be implemented. In contrast, the digital quantum simulation is implemented through quantum gates acting on either one or two qubits in the system [27][28][29]. The main advantage of digital quantum simulation is the capability of performing universal computation, i.e. implementing * Electronic address: chufan.lyu@std.uestc.edu.cn † Electronic address: vmontenegro@uestc.edu.cn ‡ Electronic address: abolfazl.bayat@uestc.edu.cn arbitrary unitary operations [36]. The major drawback is the high fragility of such simulators to noise and their demand for complex error correction schemes [37]. Many important phenomena in condensed matter physics, material science, and chemistry are explained from the ground state of a certain Hamiltonian. This includes electronic structures of matter [38], molecular formations [39], magnetization [40] and quantum phase transitions [41]. There has been a significant effort to simulate the ground state of many-body systems on classical computers. Semi-classical approaches, such as density functional theory [42], are extensively used for characterizing the electronic structures but fail to explain all quantum effects. Quantum Monte Carlo [43] has been developed to go beyond such mean field approaches but they also suffer from sign problem. Density matrix renormalization group based methods [44], are perhaps the most successful approach for characterizing the ground state features of many-body systems, although they are mostly limited to 1-dimensional systems. Machine learning [45][46][47][48][49][50] has also been exploited to capture the ground state features, but it is hard to scale them up for large systems. The limitations of these approaches truly reveals the limited power of classical simulations for unraveling quantum features.
The emergence of quantum simulators has opened new vistas for simulating many-body systems. To date, there are three different methods for simulating the ground state of many-body systems on quantum simulators: (i) imaginary time evolution [51,52]; (ii) adiabatic evolution [35,53]; and (iii) Variational Quantum Eigensolver (VQE) [54]. The first two approaches are entirely performed on quantum simulators, and thus their success heavily relies on the quality of the quantum hardware. In particular, for large systems, these approaches either demand deep quantum circuits or long coherence times. However, the current Noisy Intermediate Scale Quantum (NISQ) simulators are prone to various types of error in their initialization, readout, and gate operations, while their finite coherence time limits both the system size and the depth of the circuit. To overcome these obstacles, the VQE method, a hybrid combination of a quantum circuit and a classical optimizer, has been theoretically investigated [55][56][57][58][59][60][61][62] and experimentally implemented on various NISQ simulators, including photonic chips [54], ion traps [13][14][15], nuclear magnetic resonance systems [63], and superconducting quantum devices [22,31,32]. In fact, the use of the classical optimizer allows one to simplify the quantum hardware into a shallow circuit. In other words, the resources needed for simulating the ground state is divided between quantum hardware and classical optimizer, enabling the imperfect NISQ simulators to emulate the ground state of many-body systems. Recently, VQE has been generalized to simulate excited states [64,65], non-equilibrium dynamics [66], Gibbs ensemble [67], and approximate solutions for combinatorial optimization problems [68] in many-body systems. Based on these, it is highly desirable to have a quantitative analysis for the required resources, either quantum or classical, in different approaches for simulating the ground state of many-body systems.
In this paper, we use a digital quantum simulator to prepare the ground state of a Heisenberg Hamiltonian. We first compare the adiabatic evolution and the VQE with respect to their demand of quantum resources, namely the depth of the circuit or equivalently the number of two-qubit entangling gates. Our results quantitatively show that the VQE is far more efficient in terms of quantum resources than the adiabatic scheme. We then focus on the classical optimization part of the VQE to speed up its convergence and thus minimize the total time needed for the performance of the VQE. We specifically develop two different strategies for improving the initial guess for the circuit parameters to start the optimization at a closer point to the global minimum. Our analysis shows that these strategies can significantly speed up the convergence of the optimization, in particular, for larger systems.

II. MODEL
Our goal is to prepare the ground state of a manybody system on a Digital Quantum Simulator (DQS). For simplicity, and without loss of generality, we consider a 1-dimensional chain of N spin-1/2 particles interacting via Heisenberg Hamiltonian where J > 0 is the exchange coupling and σ i = (σ i x , σ i y , σ i z ) is the vector of Pauli operators at site i. The Heisenberg model is one of the key models in both condensed matter and quantum information technologies, and thus, has been extensively explored in both ground state [69] and non-equilibrium dynamics [70][71][72][73]. The Heisenberg Hamiltonian has SU(2) symmetry and commutes with the total spin in the α = x, y, z direction (i.e. [H, S α ] = 0 where S α = i σ i α ). In addition, the Hamiltonian also commutes with the total spin operator S 2 tot = S 2 x + S 2 y + S 2 z . This implies that every eigenstate of the system has a specific total spin S tot with a fixed S z . In particular, for even N , assumed in this paper, the ground state |GS is a global singlet with both S tot = 0 and S z = 0. Furthermore, the Heisenberg Hamiltonian supports mirror-inversion symmetry as [H, M] = 0, where M is the mirror-inverting operator with M|i 1 , i 2 , · · · , i N = |i N , · · · , i 2 , i 1 . This implies that the eigenstates of the system also follow the same symmetry. It is worth emphasizing that our methodology is easily generalized to other Hamiltonians, such as XXZ, which satisfy the conservation of S z and mirror-inversion symmetries. The ground state of the Heisenberg Hamiltonian has a matrix product state representation so that the results from quantum simulators can be certified by classical simulation. This allows us to enhance our confidence on the performance of quantum simulators, and paves the way for solving more complex problems which are not classically simulable, such as the ground state of 2-dimensional systems.

III. DIGITAL QUANTUM SIMULATOR
A digital quantum simulator is a quantum machine which manipulates the quantum state of its qubits through performing certain one-and two-qubit gate operations. In order to be a universal quantum computer, a digital quantum simulator needs to be capable of applying arbitrary single qubit unitary operations (which can be decomposed into rotations around x, y and z axes) on all qubits together with at least one two-qubit entangling gate between any pair of particles [36]. In this paper, we consider all the single qubit unitary operations to be a combination of the following gates P (θ) = 1 0 0 e iθ , and R α (θ) = e i θ 2 σα , with α = x, y, z, (2) where P (θ) is the phase gate and R α (θ) is the rotation around the α-axis. For a two-qubit entangling gate we use controlled-x (also known as CNOT) gate which is defined as By properly combining the above one-and two-qubit gates one can make any unitary operation. Nonetheless, there is not a straightforward recipe to design an efficient circuit to generate a specific target state. In particular, designing protocols for preparing an eigenstate of a Hamiltonian, on digital quantum simulators, has attracted significant attention in recent years [13, 22, 27-29, 31, 32]. We pursue two approaches for generating the ground state of many-body systems, namely adiabatic approach and VQE. In both methods, which will de discussed in the following sections, the entangling part of the circuit is a two-qubit unitary operator in the form of N (θ x , θ y , θ z ) = e i(θxσx⊗σx+θyσy⊗σy+θzσz⊗σz) (4) in which θ x,y,z are three angles. One can actually implement this unitary operator using a simple circuit shown in Fig. 1(a), in which three CNOT gates are combined by five local rotations around y and z axes [74]. Thanks to the symmetries of the Heisenberg Hamiltonian, throughout this paper, we consider θ x =θ y =θ z and thus we omit the dependence on three independent parameters and simply use N (θ).

IV. ADIABATIC EVOLUTION
Adiabatic theorem [53] is a well-known procedure for the preparation of the ground state of a complex Hamiltonian starting from a simpler one. The system is initially prepared in the ground state of a simple Hamiltonian H ad (0) which is then slowly varied in time to a more complex desired Hamiltonian H ad (T max ) after time t = T max . If the variation of the Hamiltonian is slow enough, the quantum state of the system always remains in the ground state of the instant Hamiltonian H ad (t) [75]. In this paper, we are focused in the preparation of the ground state of the Hamiltonian H in Eq. (1) using a digital quantum simulator. In Ref. [35], in an analog approach, a dimerized Heisenberg Hamiltonian is adiabatically evolved to the uniform Hamiltonian H according to where 0 ≤ t ≤ T max and (6) At t = 0 the Hamiltonian is fully dimerized with only odd couplings switched on and the ground state of the system is simply By choosing T max ∼ N 2 the condition of adiabatic theorem [75] is satisfied and at t = T max , where H ad (T max ) = H, the quantum state of the system becomes the ground state of the desired Hamiltonian H [35].
To obtain the time evolution of the system, governed by H ad (t), on a digital quantum simulator one has to adopt two essential steps: (i) discretize the time evolution into M ad steps over which the Hamiltonian remains fixed for a time interval of ∆t = T max /M ad ; and (ii) exploit the Suzuki-Trotter expansion for evolving the system in each time step. In each discretized time step k (with k going from 1 to M ad ) the Hamiltonian is considered to be fixed, namely H ad (k∆t). The time evolution operator at time step k is thus written as Notice that, in order to emulate the true time evolution of the H ad (t) the discrete time steps have to be small or equivalently the number of steps M ad has to be large. This will be discussed in more details in the following. As mentioned above, the second essential step for realizing the time evolution of a many-body system on a digital quantum simulator is to use Suzuki-Trotter expansion [76]. This will allow us to simulate the dynamics only through one-and two-qubit gate operations. We consider both the first and the second order realization of the Suzuki-Trotter expansion for the unitary operator U k (∆t), given in Eq. (7), at the time step k , namely While the first order Suzuki-Trotter operator U ST1 k (∆t) approximates the discrete time evolution U k (∆t) by a quadratic error as we present the quantum circuit for each time step of the first and the second order Suzuki-Trotter evolution, respectively.
In order to speed up the adiabatic protocol, we first determine the minimum time T max which is required to prepare the system in its ground state for a given threshold fidelity F = | Ψ(T max )|GS | 2 between the quantum state of the system |Ψ(T max ) and the target state |GS . By fixing the threshold fidelity F and system size N , one can solve the Schrödiner equation with time dependent Hamiltonian H ad (t) and find the minimum T max which results in a fidelity above the threshold F . As shown in Ref. [35], the minimum time needed scales as T max ∼ N 2 where the proportionality coefficient gets larger as the threshold fidelity increases. By specifying the minimum value of T max one can then apply the above adiabatic protocol on a digital quantum simulator to find the minimum required Suzuki-Trotter steps M * ad to achieve the threshold fidelity F . In Figs. 2(a)-(b), we plot M * ad as a function of N for various threshold fidelities when the first and the second order Suzuki-Trotter circuits are used respectively. Notably, for a given threshold fidelity, the number of layers in the the second order Suzuki-Trotter approximation is almost one order of magnitude smaller than the the first order one. However, the number of gates in each layer of the second order Suzuki-Trotter circuits are more than the number of gates in their corresponding first order. Hence, it is more meaningful to compare the number of gates needed for achieving the threshold fidelity rather than the number of layers. In practice, the two-qubit gates are far more challenging than the one-qubit gates. Therefore, in Figs. 2(c)-(d) we plot the number of CNOTS versus the system size N for various threshold fidelities in both first and second order Suzuki-Trotter circuits, respectively. Remarkably, the second order Suzuki-Trotter circuit shows a clear superiority over the first order by demanding significantly less number of CNOT gates for delivering the same fidelity. This clearly shows that the second order Suzuki-Trotter expansion is significantly more resource efficient than the first order for realizing the ground state of a many-body system on a digital quantum simulator.

V. VARIATIONAL METHOD
As shown in the previous section, the adiabatic approach for creating the ground state of many-body systems demands deep circuits, i.e. large M * ad , with a considerable number of gates. The near-term quantum simulators cannot provide such deep circuits due to their imperfect noisy quantum gates, in particular the twoqubit gates. Therefore, it is worthy of exploring the possibilities for the use of shallow circuits to simulate the ground state of many-body systems. Variational methods, such as variational quantum eigensolver (VQE) algorithm, provide an alternative to adiabatic evolution which can be realized in shallow circuits. In these methods, the post processed measurement outcomes of a fairly shallow quantum circuit has to be fed into a classical optimizer iteratively to find the optimal circuit parame- ters for realizing the ground state of a many-body system. In the VQE algorithm, the quantum device prepares a quantum state |ψ( θ) , which is often called the ansatz, with θ = (θ 1 , θ 2 , · · · , θ L ) being some controllable parameters in the circuit. Then one can measure the average energy H = ψ( θ)|H|ψ( θ) for this quantum state through appropriate measurements on the qubits of the simulator. For the Heisenberg model in Eq. (1), this reduces to measure correlation functions between adjacent qubits and then add them up, namely H = N −1 i=1 ψ( θ)|σ i · σ i+1 |ψ( θ) . The measured average energy can be used as an input for a classical optimizer to find its minimum by adaptively adjusting the parameters θ on the simulator. The minimization algorithm eventually finds the optimal parameters θ min for which |ψ( θ min ) is the ground state with minimum average energy. Therefore, the resources required for VQE is divided between the quantum simulator and the classical optimizer. In fact, in each Classical Optimization Iteration (COI) a new experiment on the quantum simulator is needed, making the VQE very time consuming. Thus, for the success of the VQE, two steps are essential: (i) a quantum circuit, as shallow as possible with minimal number of parameters; and (ii) an efficient classical optimizer which can find the minimum of average energy using least number of iterations.
Designing a quantum circuit, with minimal resources, which is capable of realizing the ground state of the Hamiltonian, here Heisenberg interaction, is a key step in a VQE algorithm [77]. The output of the circuit will be |ψ( θ) and by varying the parameters θ one can explore the Hilbert space for finding the quantum state with min- N=4 15  135  3  45  2  18  N=8 113  2373  18  594  3  63  N=10 247  6669  32  1344  3  81  N=16 770 34650  79  5451  5  225  N=20 1330 75810 132 11484  6  342   Table I: A comparison for the circuit depth M * between the adiabatic approach (with both first and second order Suzuki-Trotter approximation) and the VQE algorithm for various system sizes when the threshold fidelity is chosen to be F = 0.99.
imum average energy. Ideally, the quantum circuit should be able to explore the whole Hilbert space as θ vary. However, using symmetries of the system one can simplify the circuit, i.e. decreasing the number of parameters in θ, and instead of exploring the whole Hilbert space only focus on the relevant part of the Hilbert space. For instance, in the case of the Heisenberg Hamiltonian, considering the preservation of S z and the mirror-inversion symmetries can significantly reduce the number of parameters in θ which not only simplifies the quantum circuit but also makes the convergence of the classical optimization procedure faster. Considering these two symmetries, we suggest to use N (θ), introduced in Eq. (4), as the entangling gates of the circuit because they do not create any extra excitation in the system. These gates can be combined in a geometry similar to the first order Suzuki-Trotter steps. In order to be able to locally manipulate the phases of each qubit we can also add phase gates to every qubit of the circuit as shown in Fig. 3(a). The mirror-inversion symmetry implies that for even (odd) N the angles of the phase gates acting on site k and N − k + 1 are equal with the opposite (the same) signs. Therefore, one can easily show that the total number of parameters L in a VQE circuit is where, M V QE represents the number of layers. The proposed circuit is capable to explore a part of Hilbert space which includes the ground state. In principle, one can chose the parameters such that this circuit replicates the adiabatic circuit in Fig. 1(b)-(c), e.g. by putting the angles of the phase gates to zero. However, in the adiabatic circuit every layer is close to identity operator as it evolves the system only for small time ∆t. This results in a deep circuit demanding a considerable number of gates as shown in Fig. 2. The idea behind the VQE is to exploit an optimization algorithm to find the optimal set of parameters for which a shallow circuit can implement the ground state of a many-body system. In order to finalize the VQE circuit, apart from designing each layer, one has to also determine the number of layers M V QE . To find out the number of layers needed to successfully produce the ground state, one can use the fidelity between the quantum state of the circuit and the real ground state at the end of optimization procedure. The layers are increased one by one so that the attainable fidelity exceeds a threshold value. Note that fidelity is not an experimentally friendly quantity as one needs to know the quantum state. Therefore, in practice, instead of fidelity one can focus on the convergence of optimized average energy as the number of layers increases. For any system size N and a chosen value for the layer number M V QE we use 50L (for L being the number of parameters in θ) COI to see whether the obtained fidelity (or average energy) reaches the threshold value. If the optimization fails to reach the threshold fidelity then we increase the layer number M V QE by one unit until the algorithm is successful.
In every circuit, we initially start with a random guess for all the parameters in θ and use Adam algorithm (with amsgrad applied) [78,79], as an extension to stochastic gradient descent, to minimize our cost function, namely the average energy of the system. In Adam algorithm there are four different hyper parameters which have to be set as the input of the optimizer: (i) the learning rate α, which controls how fast the parameters are going to be updated in each iteration. (ii) β 1 , the exponential decay rate for the first moment estimation; (iii) β 2 , the exponential decay rate for the second-moment estimation; and (iv) a very small number to prevent any division by zero. We set the hyper parameters to be α = 0.01, β 1 = 0.9, β 2 = 0.999 and = 10 −8 .
In order to be independent of the initial random guess for the parameters θ we repeat the procedure for several (∼ 100) different random samples. Only if the fidelity of all of them exceeds the threshold after 50L iterations we call it a success. Similar to the adiabatic approach, the minimum number of layers which can successfully generate the ground state within the threshold fidelity is called M * V QE . In Table. I, we compare M * V QE with M * ad in both first and second order Suzuki-Trotter approximation and their corresponding CNOT gate numbers for different system sizes when the threshold fidelity is chosen to be F = 0.99. As it is clear from the Table, the number of required layers and the CNOT gate numbers in the VQE approach is significantly smaller than the adiabatic procedure. Even for a system of size N = 20, a fairly shallow circuit of M * V QE = 6 layers suffices to produce the ground state of the Heisenberg chain with the fidelity exceeding F = 0.99. This remarkable observation shows the superiority of the VQE over the adiabatic approach with respect to the complexity of the circuit. The simplification in quantum resources, i.e. the circuit, is achieved thanks to the employment of a classical optimizer. Thus, in order to truly investigate the complexity of the VQE approach, one has to consider the required resources for the classical optimization procedure as well. The qubit recursive algorithm for creating a large circuit using the optimized parameters of two subsystem with smaller sizes, here half size. Apart from the optimized gates, one entangling gate, shown in green, initialized by a random number is added between the two subsystems for connecting them. If the larger system needs more layers than the smaller ones then the extra layers are initialized with the parameters copied from the previous layer. (c) The layer recursive algorithm for making the circuit deeper step by step after optimizing each layer. The optimized parameters of the previous layers are used as the initial guess for the new circuit when a new layer is added. The parameters in the new layer are copied from the previous one.

VI. STRATEGIES FOR ACCELERATING THE CLASSICAL OPTIMIZATION
A conventional classical optimization procedure starts with a random guess for θ and iteratively minimizes the cost function, here the average energy. However, the random initial guess can be very far from the real minimum of the cost function. This may lead to two unwanted outcomes: (i) many optimization iterations are needed to get to the global minimum; and (ii) it is likely that the optimization gets trapped in a local minimum and thus fails to reach the right answer.  Therefore, it is of high interest to see whether one can start the optimization with a smart guess which is relatively close to the real minimum of the cost function. In the following, we first explain the random initial guess in details, as the first strategy, and then provide two different strategies with which we can speed up the convergence and improve the precision.
Strategy 1: random initialization. We first provide a detailed introduction for the basic random initialization strategy. A single circuit layer which is used in this strategy is shown in Fig. 3(a). One has to cascade several layers in order to achieve a high quality output. In a circuit of N qubits and M V QE layers there are L number of parameters, which is given in Eq. (9). In this strategy, the initial values of the parameters are sampled from a Gaussian distribution of mean 0 and variance 1. In Fig. 4(a) we plot the average energy as a function of optimization iterations for circuit of different layers. As the figure shows by increasing the iterations the average energy decreases steadily, however, if the number of layers are small the final converged value is far from the ground state energy. This shows that if the circuit is too shallow then it has no capacity to reach the ground state merely by classical optimization. As the number of layers increases the converged value is closer to the real ground state energy. In the Fig. 4(b) we give the corresponding obtainable fidelity averaged over 100 different random samples as a function of COI for circuits of different circuit layers M V QE . As the average energy approaches the real ground state the fidelity also goes to 1. In order to see the stability of the protocol, we repeat the optimization for 100 different random samples. In Fig. 4(c) we plot the average of obtained fidelities together with the error bars, given by the standard deviation, as a function of iterations for a system of length N = 20 and a circuit with M V QE = 6 layers. Even after ∼ 7500 iterations the average fidelity does not reach 0.95 showing that the random initial choice is not an efficient guess for the optimizer.
Strategy 2: qubit recursive. The first strategy to improve the classical optimization procedure is to use the optimized value of a smaller system as the initial guess for the larger one. Here, we focus on the case that we use the optimized parameters in a circuit of size N/2 as the initial guess for a system of size N . Let's first consider the case that the number of layers M V QE is identical for both the systems of size N/2 and N . In this case, we simply replicate the optimized circuit of size N/2 twice to make a larger system of size N . To have a complete circuit one has to add one entangling gate N (θ) to each layer at the connection between the two copies acting on odd couplings, as shown in Fig. 3(b). Therefore, in total, one has to add M V QE different entangling gates between the copies which are initialized with random guesses. Apart from these newly added gates, the rest of the one-and two-qubit gates are all initialized with the optimized value from the system of size N/2. The motivation behind this choice is that, the reduced density matrix of the subsystems from the ground state of the large system has large overlap with the ground state of the local system. Thus, using this strategy, one expects to converge to the real global minimum faster. In the case that the number of layers in the small chain is less than the layers required for the large system we simply repeat the same parameters of the last layer for the added layers. The convergence of the average energy versus the COI is shown in Fig. 5(a) for different circuit layers. In Fig. 5(b) the corresponding obtainable fidelity averaged over 100 different random samples are depicted as a function of COI. Indeed, the results show faster convergence for each choice of layers. In order to see the overall performance, we repeat the optimization for 100 different initial random guesses in a system of size N = 20 with M V QE = 6 layers. The average obtained fidelities, with corresponding error bars, are plotted in Fig. 5(c) as a function of iterations. The performance shows clear improvement in comparison with the random strategy as the average fidelity exceeds the threshold fidelity of F = 0.95 only after ∼ 3500 iterations.
Strategy 3: layer recursive. In this strategy, we start a circuit with only one layer. Obviously this circuit is too shallow to provide the ground state for large systems. Nonetheless, one can perform the optimization using random initial guess for minimizing the average energy. The minimization converges normally fairly quick, as the circuit is shallow and the number of parameters are not that large. However, the output state may have a small fidelity with the real ground state as the circuit is shallow. Hence, in the next step we add one more layer, which means that the circuit becomes deeper and more parameters exist to be optimized. The initial values for the parameters of the newly added layer are copied from the previous layer. At the end of the optimization, one can add a new layer and repeat this procedure. A schematic of the procedure is shown in Fig. 3(c). In Fig. 6(a) we plot the average energy versus the number of iterations for different layers. As the figure shows while the energy converges by increasing the iterations the final value reaches the ground state energy only for circuits with adequate depth. In order to see the overall performance,   we repeat the procedure 100 times and compute the final fidelities. In Fig. 6(b), we illustrate the obtainable fidelity averaged over 100 different random samples as a function of COI. As evident from the figure, this strategy shows faster convergence as contrasted with the previous ones. In Fig. 6(c) we plot the average fidelities, with the corresponding error bar, as a function of COI in a circuit of size N = 20 and depth M = 6. As the figure shows the speed of convergence is sensibly faster than the random initial guess, discussed in the strategy 1, and the average fidelity exceeds the threshold fidelity F = 0.95 only after ∼ 2500 iterations.
In order to compare the three VQE strategies, in Table II, we report the average fidelities obtained in each strategy for three different iterations in a system of size N = 16 and a circuit of M = 6 layers. The data clearly shows significant improvement of both qubit and layer recursive strategies over the random one, in particular, when the number of iterations are smaller. Interestingly, the layer recursive strategy also shows modest superiority over the qubit recursive strategy. This means that optimizing the system layer by layer and use their optimized data as the initial guess for building a deeper circuit significantly enhances the convergence of the classical optimization in a VQE algorithm. This is the main result of our paper.

VII. DECOHERENCE
In the NISQ era, the quantum simulators are noisy and their performance is affected by the presence of decoher-ence. There are several attempts to simulate the effect of noise in existing digital quantum simulators [32,66,[80][81][82][83][84][85]. Here, we consider two types of errors: (i) unitary ones due to imperfect gate pulses; and (ii) non-unitary errors due to decoherence. Regarding the gate errors, the operations of single qubit rotations are very precise. The performance of CNOT gates, however, may not be as perfect as it expected due to imperfect pulses which couples the two qubits. To see the quality of such errors we consider a particular realization of the CNOT gate through Ising type interaction of the two qubits 2 is the Hadamard gate, and φ should ideally be zero but due to imperfect pulses in the device is a random number uniformly distributed over the interval [−h, +h] with zero average. To see the impact of this error, we consider a VQE circuit for a chain of length N = 10 with M V QE = 3 layers using optimized parameters, obtained for the ideal case. In reality each CNOT gate will have a random phase φ in the circuit. We consider 100 different realization of such random phases and average the final fidelity for each value of h, as the noise strength. The average fidelity, using 100 different random realizations, versus noise strength h is plotted in Fig. 7(a). As the figures shows even for a relatively strong noise of h = 0.1 the fidelity still remains above 0.8.
A more serious form of noise in NISQ simulators is dephasing which appears as a non-unitary operation and destroys the superposition of quantum states. Dephasing cannot be digitalized and thus it is difficult to incorporate it in our circuit picture. However, in order to have an approximation of its effect, after each layer of gates we evolve the system according to a Lindblad master equation for a certain time ∆t which is equal to the time needed to perform the gates in the previous layer. Thus the output of the system after M layers will be where, ρ (M ) (ρ (M −1) ) represents the density matrix of the system at the output of the layer M (M − 1), U (M ) represents the unitary operation of all the gates in layer M and L represents the Lindblad master equation ∂ t ρ = L(ρ) with The coefficient γ shows the strength of the decoherence. In Fig. 7(b), we plot the fidelity as a function of γ∆t for a VQE circuit of length N = 10 and depth M V QE = 3 with all the parameters already optimized. As the figure shows up to γ∆t = 0.0125 the fidelity remains above 0.5 which reveals considerable contribution of the ground state in the output of the quantum simulator.

VIII. CONCLUSION
In this paper, we have investigated two different strategies, namely the adiabatic evolution and the VQE, for simulating the ground state of many-body systems on digital quantum simulators. Interestingly, simulating the ground state is more challenging in practice than the nonequilibrium dynamics. While we have focused on the simulation of the ground state of the Heisenberg model, our results can easily be generalized to any Hamiltonian which preserves the number of excitations and satisfies the mirror inversion symmetry, e.g. XXZ Hamiltonian. Considering the complexity of the quantum circuit, our results quantitatively show that the VQE approach demands shallower circuits with significantly less number of two-qubit entangling gates in compare with the adiabatic evolution. This makes the VQE more suitable for the NISQ simulators. The simplicity in the VQE circuit comes with the price of the necessity for a classical optimizer which may demand a significant amount of iterations to succeed. In each COI, one set of experiments have to be performed making the VQE very time con-suming. As our main results, we have developed two approaches, namely qubit and layer recursive methods, for accelerating the convergence of the classical optimizer. Both of these approaches, try to start the optimization procedure with an initial guess close to the global minimum and thus reduce the total number of iterations as well as the chance of getting trapped in local minimums. Our results show that, indeed, these strategies significantly reduce the COI required for the convergence of the VQE algorithms. Finally, we have also studied the effect of decoherence for the VQE algorithm. In particular, we have determined how the final fidelity is affected by imperfect unitary gate errors and the non-unitary dephasing effects.