Finding the ground state of the Hubbard model by variational methods on a quantum computer with gate errors

A key goal of digital quantum computing is the simulation of fermionic systems such as molecules or the Hubbard model. Unfortunately, for present and near-future quantum computers the use of quantum error correction schemes is still out of reach. Hence, the finite error rate limits the use of quantum computers to algorithms with a low number of gates. The variational Hamiltonian ansatz (VHA) has been shown to produce the ground state in good approximation in a manageable number of steps. Here we study explicitly the effect of gate errors on its performance. The VHA is inspired by the adiabatic quantum evolution under the influence of a time-dependent Hamiltonian, where the - ideally short - fixed Trotter time steps are replaced by variational parameters. The method profits substantially from quantum variational error suppression, e.g., unitary quasi-static errors are mitigated within the algorithm. We test the performance of the VHA when applied to the Hubbard model in the presence of unitary control errors on quantum computers with realistic gate fidelities.


INTRODUCTION
Simulating systems of strongly correlated electrons, such as the iconic Hubbard model, is a key goal of condensed matter physics. But important effects, such as as high-T C superconductivity or detailed magnetic properties still pose serious computational challenges. The hope is that digital quantum computers or quantum simulators would bring the needed progress. The Hubbard model and spin models have been studied in several proposals and experiments, e.g., with ultra-cold gases [1][2][3] and trapped ions [4][5][6][7]. These experiments can be considered analog simulations, where the system to be studied is recreated by a well controllable artificial one. The goal is to simulate systems which are beyond the reach of classical computations. But so far classical simulations can match all existing fermionic analog simulators, and the experiments on the fermionic Hubbard model -while representing impressive technological advances -are still at the proof-of-principle state. One of the problems is that analog simulators based on fermions are limited to high temperatures as compared to the intrinsic coupling strengths [8].
In recent years, systems with increasing numbers of high-fidelity and fully controllable Josephson qubits have become available, and they were integrated in a single processor. This opens the perspective of simulating the Hubbard model, e.g., its time evolution or correlation functions, using a gate-based approach [9][10][11]. Qubits with fidelities at the threshold for the implementation of quantum error correction have been demonstrated [12,13]. However, for the near-term prospects the number of qubits required for full quantum error correction is prohibitively large [14][15][16]. Hence, for meaningful nearterm applications it is crucial to estimate the effects of errors [17]. For certain situations, methods to verify the performance of quantum simulators with errors have been suggested [18,19], and some proposals for error reduction exist [20,21].
For one of the important goals, the simulation of the ground state of a quantum system, it has been suggested and demonstrated in few-qubit experiments that variational algorithms require only a relatively low number of gates and, in addition, variational methods intrinsically suppress the impact of errors [22][23][24]. In general, variational approaches apply a unitary operator to an initial state |ψ 0 that is easy to prepare. The unitary operator U (θ) depends on a set of parameters θ that is varied to minimize the energy where H is the Hamiltonian of the system of interest. In this paper we study explicitly the effect of gate errors on a variational algorithm for finding the ground state of the Hubbard model. We will use a specific variational ansatz, namely the variational Hamiltonian ansatz (VHA) [25]. It is inspired by the adiabatic ground state evolution as explained in more detail below. Specifically we address the following questions: How close can E(θ) get to the ground state energy E g of the Hamiltonian H, and how close can U (θ)|ψ 0 approximate the true ground state of H, if gate errors occur during the implementation of the unitary operator U (θ).
Generally the specific nature of gate errors is not known, therefore we work with a simple but representative model. Every gate can be interpreted as a rotation of the qubit register. In our model gate errors are modelled as over-rotations (or under-rotations). As discussed in earlier work [17] the over-rotation angle δϕ can be related to the minimal gate fidelity F min of the gate via F min = cos(δϕ). 1 Because of the vanishing slope of the cosine at the maximum (i.e., |δϕ| = arccos(F min ) ≈ 2(1 − F min ) for F min ≈ 1) gate fidelities need to get very close to 100 % to significantly limit the magnitude of the over-rotations.
Below we also compare the VHA to the adiabatic state preparation based on the Trotter expansion. We find that the VHA produces a better approximation to the ground state with far fewer steps, and therefore gates, than adiabatic state preparation. For adiabatic state preparation, even for weak gate errors, upon increasing the number of steps, the states created have decreasing overlap with the actual target ground state. In contrast, the VHA achieves high overlap with the exact ground state; even with gate errors. This is due to the error mitigation capabilities of variational approaches. For the (still smallsize) Hubbard model considered as an example in the following we find that a gate fidelity of F min = 99.9 % is sufficient for a meaningful simulation.
where θ collects all the variational parameters θ α,k . The optimization criterion is the minimization of the energy expectation value (1) of the final state |ψ f = U (θ)|ψ 0 with respect to the n · N variation parameters θ (the ground state is, per definition, the state with minimal energy). The ansatz (2) is inspired by the adiabatic time evolution under the influence of the Hamiltonian H = H 0 + V composed of, e.g., a non-interacting part H 0 and the interaction V . If the interaction is turned on slowly on the time scales given by the inverse energy scales of the Hamiltonian, the initial ground state |ψ 0 of H 0 develops adiabatically into the ground state |ψ g of H. To simulate this evolution in a Trotter expansion the time τ of the evolution is divided into a large number n of Trotter time steps τ /n, each shorter than the inherent time scales, leading to During each of the short time steps one further decomposes the Hamiltonian into sub-operators. In a simulation using an available quantum computer the suboperators are chosen such that the short time evolution can be realized by the available gate operations. The similarity between the operators (2) and (3) justifies the expectation that the VHA can produce the evolution from a ground state |ψ 0 of the non-interacting system to the ground state |ψ g of the full Hamiltonian. In addition, by introducing variational parameters the VHA can deviate from the adiabatic path and follow, through optimization, a more efficient one. Having a more efficient evolution via VHA allows for greatly reducing the necessary number of steps n, as compared to the number of Trotter steps in an adiabatic evolution, while still achieving high accuracy. 2 Moreover, by optimizing the variational parameters one also mitigates the error introduced by faulty gates, an effect which had been termed variational error suppression.

MODEL HAMILTONIAN, ITS DECOMPOSITION, AND MAPPING TO QUBITS
The model system we investigate in this paper is the Hubbard Hamiltonian of spin-1 2 fermions with hopping amplitude t between nearest neighbors j, j , on-site energy U , and c j,s being the annihilation (creation) operator of a fermion on site j with spin s. In the following we consider two-dimensional square lattices and focus on the parameter values U = 2t with repulsive on-site interaction, U > 0.
For the implementation of the variational unitary operation of Eq. (2) we separate the Hamiltonian into N = 5 parts: The non-interacting part is split into four terms, H 1 , . . . , H 4 , as illustrated by Fig. 1. We distinguish between hopping terms in horizontal and vertical direction The arrows indicate the hopping terms between neighboring sites. For the simulation we divide them into four sets: First in horizontal (solid lines) and vertical (dashed lines) directions. Then we subdivide for each direction into even (green) and odd (orange) terms. Summing up all hopping terms for each of the four sets yields the sub-Hamiltonians H1, . . . , H4. Note that the hopping terms within a set commute among each other. and for each direction we group even and odd terms, i.e., every other term in each direction of the 2D system. The on-site interaction terms are collected in H 5 . Note that all terms collected within one H α commute among each other. Hence, the execution of an exponential e iθ α,k Hα in Eq. (2) can be performed exactly by sequentially applying the gates that account for the individual terms, without introducing an error associated with the Trotter expansion.
As written explicitly above we introduce a manageable number of 5 variational parameters per step and -as the example below shows -we reach high-quality results already for 10 steps or less.
The mapping onto a qubit system is performed via the Jordan-Wigner transformation. Gates performing on-site interaction terms are realized by ZZ-type interactions between qubits, hopping terms require XX + Y Y interactions. In addition the transformation introduces Jordan-Wigner strings which are implemented by additional chains of controlled Z gates. A precise description of the mapping and gate sequence of the algorithm is presented in the Appendix.

ERROR MODEL AND PROCEDURE
The goal is to prepare the ground state |ψ g of the Hubbard model (4) on a quantum computer. We start from the ground state |ψ 0 of the non-interacting model (U = 0), 3 which -in principle -can be prepared effi-ciently on a quantum computer [26], and we apply the variational Hamiltonian ansatz in order to evolve this state towards |ψ f , which should be close to the ground state |ψ g . We modeled the algorithm and the quantum gates including the gate errors on a classical computer. For the system sizes considered, |ψ 0 and |ψ g (without errors) can also be found exactly through classical numerical diagonalization. This allows us to test the quality of the results.
The unitary evolution is implemented by a gate based algorithm described in more detail in the Appendix. Each gate can be written in the form e iϕA with a real angle ϕ and an operator A composed of Pauli operators. Hence, it can be interpreted as a rotation of the quantum state. We model unitary gate errors by overrotations δϕ (which may be positive of negative), such that the faulty gate reads e i(ϕ+δϕ)A . The magnitude of the random δϕ is given by the minimal gate fidelity F min , where F min = cos(δϕ) [17]. Performing a sequence of gates with random normally distributed over-rotations with zero mean and variance Var(δϕ) one finds -for weak over-rotations, i.e., fidelities close to one -an averaged minimal gate fidelity F min = 1 − Var(δϕ) 2 /2. When we introduce gate errors in the following, we assume a certain gate fidelity F min and add random over-rotations to each gate according to the above relation. However, once the over-rotation for a specific gate is chosen, this value is kept constant during the consecutive stages of the optimization process. This accounts for quasi-static errors, which are considered an appropriate noise model for superconducting qubits, where the noise spectrum is dominated by low frequencies [27].
Finally, for given gate fidelities, we measure the performance of the VHA by evaluating the final state fidelity. This quantity is defined as the absolute value of the overlap | ψ g |ψ f | between the exactly known ground state |ψ g and the final state |ψ f = U (θ)|ψ 0 of the VHA according to Eq. (2), after the optimization of the parameters θ. 4

VHA VERSUS ADIABATIC EVOLUTION
We first study the quality of the variational Hamiltonian ansatz, Eq. (2), in comparison to the adiabatic evolution, Eq. (3). In both cases we start from the same initial state |ψ 0 . For the comparison, the same number of steps n (steps of the VHA or Trotter steps) and the same gate sequence is used (with appropriately different parameters). For low n, the gate sequence introduces an error while implementing the exponential e −i τ n H0 for the adiabatic evolution, since the summands of H 0 do not commute. But for the comparison with the VHA at equal gate count we did not introduce finer Trotter time steps. On the other hand, the time τ in the adiabatic evolution was optimized for the given number of Trotter steps n. Table I shows the results for a two-dimensional Hubbard lattice with 2×2 sites. We present the final state fidelity | ψ g |ψ f | in percent after performing, on one hand, the adiabatic evolution and, on the other hand, the VHA for different numbers of (Trotter) steps n and different averaged minimal gate fidelities F min . The algorithm we used requires 20 two-qubit gates per (Trotter) step, i.e., in total we perform from 40 to 100 two-qubit gates.
Even for perfect gate fidelities of 100 % the adiabatic evolution does not improve the state fidelity significantly if we increase n. (Note that the initial state fidelity | ψ g |ψ 0 | = 98.87 % is already high due to the small system size; no improvement could be achieved for n = 2.) For a lower gate fidelity of 99.99 % the adiabatic evolution barely increases the final state fidelity above 99 %. For still lower gate fidelity of 99.9 % the adiabatic evolution fails to improve the final overlap altogether. This was to be expected from the results of our previous work [17], where we provided estimates for the maximum number of gates that can be handled for a given gate fidelity in a quantum simulation with fixed parameters. Indeed for a gate fidelity of 99.9 % the gate count for the adiabatic evolution exceeds this limit.
On the other hand, we observe significantly better performance for the VHA. For perfect gates two steps give already a very high final state fidelity; only three steps are necessary to achieve essentially a perfect result (an error of about 10 −12 was observed, which is within numerical inaccuracies). Even with gate errors present we achieve high final state fidelities. The numbers clearly show that introducing more steps helps suppressing the quasi-static errors considered here.
The rightmost column, labeled by 99.90 * , illustrates the error mitigation provided by the VHA. For the data in this column we took the optimized parameters for perfect gates and used them in the evolution according to Eq. (2) with faulty gates with F min = 99.9 % without any further optimization. This procedure does not take advantage of the potential of the VHA for error mitigation. The low performance of this method illustrates the power of variational error suppression.
The set of gate errors are chosen random but fixed (static) for both methods. However, in different runs they are chosen independent corresponding to the given gate fidelities. The results of Table I are averaged over many runs with different sets of errors. We can add that the error suppression of the VHA also reduces the standard deviation of the results for different error sets significantly. For the adiabatic evolution and lower gate fidelities we needed of the order of 10 5 runs in order to reach the shown accuracy of the average. The VHA, on the other hand, needs only a low number of runs for larger n, even for low gate fidelities, to achieve the same accuracy, and even a single run is already quite reliable. 5

SCALING UP
Next we extend the analysis of the variational Hamiltonian ansatz to larger systems. Table II shows the final state fidelity | ψ g |ψ f | for various step numbers n and gate fidelities F min , now for the VHA applied to a 3×2 and a 3×3 Hubbard model. Here we show the results obtained for a single realization of the gate errors for each gate fidelity and step number. Statistical fluctuations are reduced due to the large number of gates per step. Averaging would introduce only small differences to the data presented. It can be ignored, especially for n > 6 where the variational error suppression is strong. For the 3×2 and 3×3 systems the algorithm requires 44 and 81 two-qubit gates per iteration step, respectively, i.e., overall up to 810 two-qubit gates were applied to the initial state.
We find again that the VHA produces very high final state fidelities. However, we also notice how high gate fidelities are necessary for a good performance of the algorithm. For a gate fidelity of F min ≤ 99.9 % which should be accessible in the next few years, we did not reach final state fidelities above 99 %. Further investigations showed that the limited final state fidelities are not necessarily a flaw of the VHA itself but rather of the required optimization. For increasing system size we found the results to be more and more sensitive to the choice of start parameters for the optimization problem. This suggests that the optimizer does not find the global optimum for the parameters but gets trapped in local extrema. The data of Table II were obtained with a rather limited set of start parameters (see the Appendix for further details).
To substantiate this conclusion we performed the VHA for some values of F min ≤ 99.9 % and n ≥ 6, exploring a larger set of start parameters (see the Appendix for more information). Table III displays the final state fidelity of the 3 × 2 and 3 × 3 systems for the optimized start parameters, showing a significant improvement of the final state fidelity over the results of Table II. We emphasize that this was not because of a favorable set of random gate errors; once better start parameters were found the results are changing little with varying the gate errors.
Another measure of the performance of VHA is to look at the value of the ground state energy. For the 3 × 3 Hubbard model with U = 2t > 0 the exact value is E g = −9.67 t. For the chosen initial state the state fidelity is already 96.18 %, but the expectation value of the Hamiltonian is only ψ 0 |H|ψ 0 = −9.29 t. After 10 steps of the VHA for a gate fidelity of 99.90 % the final state fidelity has improved close to 99 %, and the expectation value of the Hamiltonian reaches ψ f |H|ψ f = −9.60 t.

CONCLUSION
In this paper we studied in detail the quantum simulation of Hubbard models of small size and with a specific type of gate errors, but our analysis still allows drawing several conclusions: (i) The variational Hamiltonian ansatz (VHA) produces the ground state wave function of the Hubbard model in good approximation with a number of steps which is much lower than the number of Trotter steps needed in an adiabatic approach.
(ii) The effects of (static) gate errors are strongly mitigated by the variational methods.
(iii) For the considered system size, gate fidelities of the order of 99.9 %,which should be within reach for stateof-the-art digital quantum computers, allow preparation of the ground state with a final state fidelity above 99 %.
(iv) This performance can be reached with a low number of variational parameters per step (5 in our case for the 2D Hubbard model).
It is clear that introducing more variational parameters, up to one parameter per gate, would enhance the variational error suppression of quasi-static errors. Eventually it leads to approaches like the variational quantum eigensolver [24]. However, introducing more parameters poses a challenge to the classical optimization routines. We found that even for our small set of parameters, the emerging optimization problem poses a substantial obstacle, since with growing system size the gradients with respect to the variation of the parameters decrease [28]. The difficulties with the optimization algorithms appear a stronger limitation of the performance of variational algorithms than a limited set of variational parameters.
Better optimization algorithms could help with further issues. One could consider non-static portions of gate errors and statistical measurement errors. Such fluctuating errors are difficult to manage for optimizers, particularly for gradient based optimization protocols. We also noted the need for a good choice of the initial set of variational parameters, as well as of the initial state. The latter could be obtained, e.g., from mean field theory. Finally more advanced quantum gate sequences implementing the terms of the Hamiltonian can lower the gate count and reduce the impact of gate errors. We applied up to 810 two-qubit gates in our examples of rather small systems. Algorithms with superior scaling behavior for Hubbard models [11,29] should be considered to improve the performance of the VHA for larger systems.

Gate sequence
For illustration we show the gate sequence producing the unitary transformation e iθ α,k Hα from Eq. (2). The different H α , introduced in the main text, contain either hopping terms −t(c † j,s c j ,s + c † j ,s c j,s ) or on-site interactions U c † j,↑ c j,↑ c † j,↓ c j,↓ . We assume that the hardware of the quantum computer allows for ZZ-like and XX + Y Y interaction and, for simplicity, unrestricted connectivity between the qubits. The terms summarized in each of the H α commute among each other making their ordering irrelevant.
For the following discussion it is convenient to absorb the spin index in a consecutive numbering of the lattice sites via (j, ↑) → j and (j, ↓) → j + M where M is the total number of sites. With this notation the Jordan-Wigner transformation becomes c j = j−1 l=1 (−σ z l )σ − j , which involves the Jordan-Wigner string In the variational approach the gate e iθU σ + j σ − j σ + j+M σ − j+M needs to be implemented with some parameter θ, which embodies a ZZ interaction (up to some single-qubit phases). To account for the gate errors we add an over-rotation δθ to the parameter, the size of which depends on the assumed gate fidelity.
The hopping terms (where, say, j < j ) transform to −t σ + j σ − j j −1 l=j+1 (−σ z l ) + H.c. . We do not assume the hardware to allow for more than two-qubit interactions. A hopping term, can be modeled by the XX + Y Y interaction, i.e., e −iθt(σ + j σ − j +σ + j σ − j ) , which we denote as t gate. Again, gate errors lead to an over-rotation δθ to be added to the parameter θ. Residual Jordan-Wigner strings can then be implemented by sandwiching the t gate with controlled Z gates (CZ) as displayed to Fig. 2. The CZ gate has again ZZ-like interaction, and we implement it as e for the control qubit j and the target qubit j . Over-rotations are introduced as an addition δϕ to the angle π.

Choice of start parameters
When using the VHA, in order to minimize the energy according to Eq. (1), one has to start with an initial guess for the parameters θ. In Table I and Table II for each value of n and F min we tried three different sets of start parameters, motivated by some physical reasoning: 1. Since the VHA is inspired by the adiabatic evolution  2. A t gate implementing an XX + Y Y interaction between qubits j and k is nested between controlled Z gates.
The CZ gates introduce a Jordan-Wigner string such that the gate sequence implements a fermionic hopping term between the orbitals j and k.
one of our choices of parameters θ α,k from Eq. (2) was to mimic the adiabatic evolution. We set θ α,k = 1 t for α ∈ {1, . . . , 4} (i.e., the hopping elements) and θ 5,k = k n 1 t (i.e., the interaction terms), with t being the hopping energy of the Hamiltonian (4). This represents the adiabatic evolution (3) during time τ = n t . 2. We chose a parameter set where not only the interaction but also the hopping is turned on gradually, i.e., θ α,k = k n 1 t for all α. 3. We noticed that the optimized parameters of the VHA usually do not resemble an adiabatic path or show steady growth, rather they are more evenly distributed. Hence, we chose a set where θ α,k = 1 n 1 t for all α. This set represents a Trotter expansion with n steps for a time evolution for the duration τ = 1 n . For all choices of initial parameters, after the optimization the final values were vastly different. We conclude that for the larger systems one often gets stuck in a local optimum. For this reason, we tried further sets of initial parameters with results shown in Table III: Firstly, we added initial parameters recreating an adiabatic evolution similar to point 1 described above, but for τ = 1 t . Secondly, we chose an even distribution, similar to point 3, but with θ α,k = r t where we varied r ∈ {0.1, 0.2, . . . , 1.0}. (We noticed that usually the magnitude of the final parameters were between zero and ∼ 1 t .) These additional start parameters improved the data of Table II towards the results of Table III. A deeper understanding how to deterministically find suitable initial parameters should be acquired. However, this is out of the scope of this work. For this paper the tested parameter sets already helped to show the capabilities of the VHA, achieving remarkable results in the considered small systems.