Simulating strongly interacting Hubbard chains with the Variational Hamiltonian Ansatz on a quantum computer

Hybrid quantum-classical algorithms have been proposed to circumvent noise limitations in quantum computers. Such algorithms delegate only a calculation of the expectation value to the quantum computer. Among them, the Variational Quantum Eigensolver (VQE) has been implemented to study molecules and condensed matter systems on small size quantum computers. Condensed matter systems described by the Hubbard model exhibit a rich phase diagram alongside exotic states of matter. In this manuscript, we try to answer the question: how much of the underlying physics of a 1D Hubbard chain is described by a problem-inspired Variational Hamiltonian Ansatz (VHA) in a broad range of parameter values ? We start by probing how much does the solution increases fidelity with increasing ansatz complexity. Our findings suggest that even low fidelity solutions capture energy and number of doubly occupied sites well, while spin-spin correlations are not well captured even when the solution is of high fidelity. Our powerful simulation platform allows us to incorporate a realistic noise model and shows a successful implementation of noise-mitigation strategies - post-selection and the Richardson extrapolation. Finally, we compare our results with an experimental realization of the algorithm on IBM Quantum's ibmq_quito device.


I. INTRODUCTION
Quantum computers are a potentially disruptive technology [1][2][3][4] promising chemical precision solutions of problems where electron-electron correlations play a major role. A number of algorithms have been proposed for both the currently non-existing quantum error-corrected devices [5,6] and the currently available NISQ devices [7]. The Variational Quantum Eigensolver (VQE) [8] and Quantum Imaginary Time Evolution (QITE) [9,10] are two NISQ front runners for solving problems in chemistry and condensed matter physics. The Hubbard model [11] is one of the most studied problems in condensed matter physics. Elegant and seemingly simple it has been the starting point in studying high-temperature superconductivity [12] and metalto-insulator transitions [13], just to name a few. Although the Hubbard model is analytically solvable in one dimension with the famous Bethe ansatz [14], a solution in two and three dimensions is generally unknown and has been the focus of half of century of numerical studies [15,16].
In this manuscript, we use an 8-site 1D Hubbard model as a benchmarking tool for contemporary quantum computing algorithms, or to be more precise the VQE with the Variational Hamiltonian Ansatz (VHA) [17,18]. The VHA leverages an ansatz that is inspired by the Hubbard Hamiltonian itself and an adiabatic way of thinking to construct the trial wavefunction. Adopting a condensed matter physics perspective, we try to quantify how much of the underlying physics is preserved with an ansatz that minimizes the energy and respects all physical symmetries of the system. Our advanced simulation platform allows us to incorporate a realistic noise model, to quantify how much is the VHA influenced by pure dephasing, spin relaxation, and gate errors and to com-pare with experimental results from the IBM Quantum's ibmq quito device. Furthermore, we attempt an error mitigation strategy called the Richardson extrapolation [19,20] -a measurement of an observable with artificially modulated noise levels followed by an extrapolation to the zero-noise case, alongside with post-selection.
The adiabatic evolution inspired VHA was first introduced in Ref. [17] and have been investigated for Hubbard models in Refs. [18,21,22]. Furthermore, studies of spin systems were conducted in Ref. [23] as well as optimization studies with the well-studied Quantum Approximate Optimization Algorithm (QAOA) which share the same structure (proposed before VHA in Ref. [24]). The study conducted here advances that of Ref. [18] with a more elaborate noise model and discussions of potential noise-mitigation strategies. Furthermore, spin correlations, noise mitigation, and a broader range of U are discussed as compared to Ref. [21] and focuses on how well certain variables are captured as compared to resource estimation as in Ref. [22].
Our findings suggest that different quantities of the physical system require different levels of ansatz complexities. Namely, even for an ansatz of relatively low depth doubly occupied site number and energy are predicted at a satisfactory level. On the other hand, long-range spin correlations are much more sensitive to the infidelities of the trial solution -even a 90% fidelity is insufficient to capture the correct long-range spin-spin behaviour. Finally, Richardson extrapolation strongly mitigates the effects of quantum noise, although we are only able to validate this for a 2-site Hubbard model (4 qubit simulation).
This paper is organized as follows, in Section II we describe the methodology, in III we show noiseless results for our estimation of energy, fidelity, doubly occupied site number and spin-spin correlations. In Section IV we discuss the noise model, post-selection, the Richardson extrapolation and experimental results before concluding in Section V.

II. METHODS
Variational Quantum Eigensolver (VQE) has been increasingly popular in quantum computing. First proposed and implemented in Ref. [8], VQE offers a method to obtain approximate groundstates of electronic systems described by an Hamiltonian H. The VQE algorithm commonly has the following structure: 1. An initial state |Ψ 0 is prepared on the quantum computer.
4. Based on the energy, a classical minimizer will optimized the parameters θ until it converges to a minimal energy.
VQE is particularly suited to NISQ devices as it allows to obtain results with short circuits but also to mitigate errors thanks to the flexibility of the variational parameters. The choice of the parametrized circuit is essential in the VQE procedure. Two key elements need to be taken into account when selecting the type of the parametrized circuit: the ability to produce the groundstate with sufficiently high fidelity, and the implementability of the circuit on contemporary hardware. These two requirements are intertwined due to the limited length, connectivity and qubit number current-day quantum circuits can have, but also the number of variational parameters a classical optimizer can handle. We can usually distinguish between two strategies when designing an ansatz: 1) the so-called hardware efficient ansatz prioritizes on two-qubit gates which can be naturally performed in the system due to a specific qubit topology [25,26]. Although providing satisfying results on problems that require a small number of qubits, this approach often fails to scale to larger problems and can suffer from optimization barren plateaus [27]. 2) Numerous approaches proposed problem-inspired ansätze which produce variational states from relevant subspaces [28][29][30][31][32][33]. Built from physical arguments, they often require more costly circuits. Furthermore, adaptive schemes have been proposed to reduce the circuit length of those methods [34,35].
Let us now present the thinking behind constructing a VHA trial solution inspired in the Hamiltonian of the problem itself. Observe two Hamiltonians, H 0 and H f . H 0 corresponds to an eigenproblem which is easily solvable and H f to a problem which is difficult to solve. The basic idea in adiabatic evolution is to prepare the groundstate |Ψ 0 corresponding to the initial Hamiltonian H 0 . Then, one can evolve the initial state in time T described by the following time-dependent Hamiltonian H(s) = (1 − s)H 0 + sH f , with s(t = 0) = 0 and s(t = T ) = 1. If the evolution rate is slow enough with respect to the gap of the system, Fock's adiabatic theorem guarantees that |Ψ(t) = T e −i t 0 H(s(τ ))dτ |Ψ 0 (where T is the time-ordering operator) remains in the groundstate of the instantaneous Hamiltonian H(s).
In general, and on a quantum computer, the timeevolution operator can be approximated with a first order Trotter-Suzuki approximation by dividing the integration over time into N time step of duration ∆τ so that N n=1 e −i(H0(1−s(n∆τ ))+H f s(n∆τ ))∆τ . Secondly, we decompose H into non-commuting parts {H α }: H = α H α . The terms contained in H α should commute among each other in order to implement them simultaneously. By using again the Trotter-Suziki approximation, we obtained the following circuit: However, to respect the adiabatic criterion long evolution times are required leading to long circuits beyond the reach of NISQ devices. In order to shortcut adiabatic paths, VHA replaces the time steps by variational parameters, with the hope of reducing the circuit depth needed to obtain satisfying states. The ansatz can be formulated by the following: where, N L is the number of layers, and θ are the variational parameters. Let us consider now the case of the 1D Hubbard model, described by the following Hamiltonian For the Hubbard Hamiltonian and in the context of trying to determine its groundstate with adiabatic methods one can choose H 0 to be the free fermion Hamiltonian H t = −t iσ (c † iσ c i+1σ + c † i+1σ c iσ ) and the final Hamiltonian to be the full Hubbard Hamiltonian, which ↑ Uu(θ u n ) corresponds to the situation of slowly increasing the interaction strength U in the system. While a one-body state can be in principle efficiently prepared on a quantum computer [36][37][38], the variational circuit can be formulated by: with Fig. 1 shows a schematic of one layer of the ansatz.

III. RESULTS
We start by comparing different ways of preparing the initial state of the Hubbard model on a quantum computer. We will restrict ourselves to one-body state as initial states due to the simplicity of their preparation. We investigate two natural choices: a free fermion U = 0 groundstate and a mean-field Hartree-Fock (HF) solution as the latter is commonly utilized in VQE studies in the chemistry domain with widely studied Unitary Coupled Cluster Ansatz (UCC) [28,39,40] (recently schemes to prepare Gutzwiller wave function was also proposed in Refs. [41,42]). In similarity to problems in chemistry the HF state provides a good approximation to the groundstate energy (see Fig. 2a).
However, often in the study of condensed matter systems, one is interested in more quantities than just the mere energy such as long-range correlations or the expectation value of the square of the total spin S 2 . The most comprehensive measure of the quality of the obtained state is the fidelity F = | Ψ(θ)|Ψ ex | 2 , where Ψ(θ) is the approximate solution and Ψ ex the exact solution. The fidelity will be 0 if the approximated solution is totally different from the exact one and 1 if these two solutions match.
The mean-field solution is obtained by transforming the two-body term in the Hamiltonian responsible for interaction into a one-body term thanks to the following approximation: i n i↑ n i↓ i ( n i↑ n i↓ + n i↓ n i↑ ). Starting from initial guess { n iσ 0 }, the average densities { n iσ } are tuned self-consistently to minimize the energy E HF = Ψ HF |H HF ({ n iσ })|Ψ HF where |Ψ HF is the groundstate of the mean-field Hamilto- iσ n iσ . The total spin operator S 2 can be express as x, y, z and i the position index. One can first express those operators in terms of fermionic operators:

the spin raising and lowering operators S
Thanks to the chosen fermion-to-qubit mapping described in the Appendix V A, one can write those fermionic operators into qubit operators which can in principle be measured on a quantum hardware.
In terms of fidelity, the free fermion solution is much closer to the exact solution in the regime of 0 ≤ U/t ≈ 10 Fig. 2b. Moreover, the free fermion solution respects the S 2 symmetry of the problem. On the other hand, the Hartree-Fock solution produces non-zero values for the total spin S 2 . Note that a Restricted Hartree-Fock calculation will return the free fermion groundstate as the Hubbard model contains only two-body interactions between opposite spin electrons from the same site (i.e. same spatial orbital). As VHA is Hamiltonian-based, the symmetries are preserved throughout the circuit (all the gates commute with the particle number operator N and the spin operators S z and S 2 ). Therefore choosing an initial state with a right set of quantum numbers avoids the variational ansatz to span non-physically relevant parts of the Hilbert space. Other strategies aiming at including symmetries in the VQE process have been also proposed in [32,[43][44][45][46].
We use an iterative approach consisting in initializing our simulations by targeting an initial value of Coulomb repulsion from the intermediate regime: U 0 = 5 (we choose t = 1). We try many random initial parameters contained in [−π/5, π/5], keep the lowest energy solution, and use the optimized parameters as initial parameters to simulate the case U 1 = U 0 ± δU . Throughout the simulations, the optimized parameters for the case U N −1 = U 0 ± (N − 1)δU are used as initial parameters for U N = U N −1 ± δU . This method is motivated by the fact that the choice of initial parameters has a significant impact on the optimized solution (as detailed in the Appendix V B), and by trying multiple starting parameters, we mitigate the effect of initial parameters choice, which improves the result for the whole range of U values and speed up the optimization routine. Fig. 3 shows the results for a N = 8 Hubbard chain. Whereas energy and fidelity errors remain small in weakly interacting regime U 2 even for limited numbers of ansatz layers N L , the non-interacting initial state U/t already providing a close match, quality of results decreases as U increases. This can be understood as the result of lower fidelity of the initial state (Fig. 2). Fig. 3 shows the effect of the number of layers N L on the fidelity of the optimized states at different values of U . To simulate a Hubbard chain in a low interacting regime, only a few layers are required to reach relatively high fidelities > 95%. However, as U increases, the overlap between the initial and the exact solution becomes smaller. Adding more layers to the parametrized circuit leads to better results, although not systematically -highlighting the complex nature of the resulting variational state. This is shown for example that the N L = 3 ansatz produces higher energy states than the N L = 5 case but leads to similar fidelity results in the region U 12. This shows that the ansatz produces polluting states that have low energies but are further away from exact solutions. Similarly, Fig. 4 shows that N L = 4 ansatz produces better fidelities than N L = 5, 6 simulations. One can note a kink in the N L = 10 results at U = 5. This shows that the multi-start parameters still resulted in a false minima value with slightly higher energy and lower fidelity. However, solutions were improved when simulating the next case U ± δU .
Despite the fact that low depth variational circuits fail to produce the exact groundstate of the 1D Hubbard chain, we can still wonder how good are those approximate solutions describing some particular expected be- haviour of our system. Fig. 3b shows the fidelity as a function for U for different depths and its corresponding values for n ↑ n ↓ = 1 N i n i↑ n i↓ , N being the number of sites. Despite having a low overlap with the exact solution, the optimized states recover pretty well the transition from delocalized to localized electrons as the on-site Coulomb repulsion increases. This shows that the ansatz builds up non-trivial correlations to a noninteracting state, and captures this basic manybody phenomenon.
Secondly, we take a look at long-range spin correlations in our variational states, defined by C z ij = S z i S z j − S z i S z j . Because two electrons of opposite spin will repel each other on the same site, increasing U will induce an antiferromagnetic order in the chain. As a result, we expect that the spin value of one end of the chain is correlated with the sites along the chain. This is shown in Fig. 5 where we get sign-alternating values for C z ij | i=0 . We observe that as U increases, the ansatz fails to capture long-range correlations.

IV. NOISY SIMULATIONS AND EXPERIMENT
VQE algorithms were introduced to extend the capabilities of current hardwares, strongly limited by noise and decoherence. A natural question when investigating the performance of an ansatz is how robust it is when noise occurs in the hardware. Here we investigate how noise degrades the output of the algorithm. To do so, we include idle noise and gate errors in our simulation. Idle noise will occur on inactive qubits that undergo spontaneous relaxation from |1 to |0 and loss of quantum information (ie loss of the phase difference knowledge) between |0 and |1 states (thus making the qubit more and more "classical"). The processes will be characterized by the typical timescales T 1 and T 2 . Gate errors will be modelled by a depolarizing channel motivated from randomized benchmarkings [47] which assume a probability of error p 1 (p 2 ) for single-qubit (two-qubit) gates. The quantum channels used to model noise are detailed in the Appendix V C.
To set the inactive duration, we assign an average duration of t 1 = 60 ns for single-qubit gates and t 2 = 425 ns, as well as T 1 = 120 µs and T 2 = 85 µs. For gate errors, we set p 1 = 0.3% and p 2 = 1% which corresponds to current-day quantum computers performances from IBM [48].
The simulations were based on density matrix calculations. Fig. 6 shows the obtained results for 2-site Hubbard system and simulations incorporating noise being compared with noiseless simulations. Whereas for the noiseless case exact solutions are obtained with only one layer of ansatz for all value of U , noise affects greatly the quality of the results. While idle noise induces an infidelity of about 2%, adding gate errors degrade significantly the results, indicating that it is the major source of error.
In an attempt of mitigating those errors, we test two error mitigation techniques: post-selection and zero-noise extrapolation. As described in Appendix V D, postselection uses the a priori known symmetries of the target solution to discard any measurement output that doesn't respect those physical constraints. In this work, we rejected the shots not respecting the number of particules and S z symmetries. To further improve results, we simulate the same circuit with error rate p 2 and 2p 2 and perform Richardson extrapolation as proposed by [19] and experimentally implemented in [20]. The idea is to increase on purpose the noise level on quantum hardware. An observable quantity Ô can be written as Ô = O * + k c k k , being the quantity representing to noise level and O * the targeted noiseless value (see Appendix V E). Therefore, by doubling the error rate of the gates, one can perform a first-order extrapolation. This gives an approximation of the noiseless value of the energy E * at each iteration of the minimization process, but also for other quantities like n ↑ n ↓ measured from the optimized state. Fig. 6 shows that both techniques can significantly improve the results in terms of energy and number of doubly occupied site values.
Finally, we test those results against real implementation of the algorithm on IBM's hardware ibmq quito [48] as shown in Fig. 6. Although the benchmark parameters were close to the noise parameters ( exp = 102µs) at the time of the experiment, the experimental results seem more affected by noise. This can be explained by the absence of read-out and shot noise in our model (the experiment's number of shots was set to 20,000), but also by the fact that additional noise occurs in the hardware like crosstalk noise. Finally, additional SWAP gates (decomposed with 3 CNOTs) are needed because of the hardware's topology. However, the experimental results are in qualitative agreement with the noisy simulations as both energy and double occupancies follow the same behaviour. While we didn't implement zero-noise extrapolation as it represents a difficult experimental challenge, post-selection improves overall the results for the energy values as well as the number of doubly occupied sites, where the steep decrease for small values of U is better captured.

V. CONCLUSION
We investigated the Variational Hamiltonian Ansatz (VHA) [17] through the example of the one-dimensional Hubbard model. This ansatz is inspired by the natural idea of adiabatic evolution, in which a free fermion initial state is driven towards the correlated groundstate by slowly increasing interaction in time. Each time steps are replaced by variational parameters in the hope to find an optimized path to the groundstate. We perform classical simulations of the VQE algorithm for an 8-site Hubbard chain which corresponds to 16 qubit simulations. We have also carried out simulations up to 12 sites (24 qubits) which are giving qualitatively equivalent results. However, we have chosen to present the 8-site chain results as it represents a good trade-off between system size and computational time.
By calculating the energy error and the fidelity of the optimized states, we quantify the performances of the algorithm. Although low numbers of layers give satisfying results in weakly correlated regimes, higher circuit depths are required to target strongly interacting groundstates. We study how well the approximated groundstates capture the physics of the Hubbard model by investigating correlations built upon the free fermion Slater determinant. We observe that short-length circuits can capture well the localization of electrons as U increases. Longrange antiferromagnetic correlations are however harder to obtain despite high fidelity states and require a larger number of ansatz layers.
Finally, we test the ansatz against noise for a 2-site Hubbard model. We include noise models in our simu-lations based on density matrix calculations and compare with experimental results from IBM Quantum's ibmq quito device. Because of two-qubit gate errors results are greatly degraded, which indicates the practical limit of Trotterized-like ansätze. Zero-noise extrapolation enables us to qualitatively improve our noisy simulation results for both the energy and the number of doubly occupied sites. Further improvements should be aimed at both proposing relevant initial states requiring short circuits still exhibiting the right symmetries as well as bridging this scheme towards hardware efficient approaches. ACKNOWLEDGMENT We are grateful to the Association Nationale de la Recherche et de la Technologie for funding. The computations were carried out on the Atos Quantum Learning Machine (QLM). We acknowledge IBM Quantum for the access of their quantum devices [48]. The quantum circuit schematics were drawn thanks to Quantikz package [49]. Note added: After submitting our manuscript, an independent preprint [50] investigated experimentally spin and charge properties of the Hubbard model by means of a 1-layer VHA circuit. The reported results are in qualitative agreement with our study.

A. Qubit mapping and gate decomposition
In order to express the ansatz in terms of a qubit circuit, one needs to map the fermionic operators onto qubit operators. Among the different existing mappings, we choose the Jordan-Wigner transformation. To implement this, one first needs to choose an ordering for the different one-body states. For a Hubbard chain of length N , we choose to index the states by site index, the N first being attributed to spin-up, the N second for spin-down states, so that c † i↑ = c † i and c † i↓ = c † i+N . We express the creation and annihilation operators as follows: Having express fermionic operators into Pauli strings, one can rewrite the operator contained in the ansatz: In order to prepare one-body states as initial states, we use a method based on QR decomposition [36,37] implemented in the OpenFermion package [38] that decompose a basis transformation into a set of Givens rotations, that  8. Schematic of the quantum circuit U 1 h (θ) for N = 6 site chain for spin α (the same circuit will be applied to spin β qubits indexed from N to 2N − 1). A gate e −i(XX+Y Y )θ/2 is applied between each pair of qubits 2i and 2i + 1 for i ∈ [0, N/2 − 1].
can be expressed as: The last step is to find a gate decomposition that suits one's hardware specification. Aiming for universality, we choose in our simulations to decompose those unitary operations in terms of single-qubit gates and CNOTs. Figs. 10, 11 show the decomposition for the unitary transformation associated to the hopping terms defined in Eq. 7 (proposed in Ref. [51]) and the interaction terms defined in Eq. 8. Figs. 7, 8, and 9 show circuits associated to the different parts of the Hamiltonian. Fig. 12 displays the gate decomposition for Givens rotation (Eq. 9) taken from [30]. In this configuration, each element of the ansatz is decomposed with 2 CNOTs. For an N-site chain (N > 2) and a ladder qubit topology, this leads to 6N -4 CNOTs per layer of the ansatz. 9. Schematic of the quantum circuit U 2 h (γ) for N = 6 site chain for spin α (the same circuit will be applied to spin β qubits indexed from N to 2N − 1). A gate e −i(XX+Y Y )γ/2 is applied between each pair of qubits 2i + 1 and 2i + 2 for i ∈ [0, N/2 − 2].

B. Optimization details
We use a multi-start procedure for initialization simulation for U 0 = 5, which consists in trying different random initial parameters and keeping the lowest energy solution after minimization. To show how sensitive are the results with respect to the initial parameters, we try 100 random initial parameters contained in [− π 5 , π 5 ] and compare the obtained results. Fig. 13 shows the average optimized energies and their associated fidelities with standard deviation.

C. Noise model
In order to model noise occurring in real devices, we use the Kraus operators representation to describe quan- tum channels associated with noise models (Eq. 10). We combine noise acting on idle qubits as well as gate errors. During their inactive periods, the qubits undergo simultaneously phase damping described by (Eq. 11) which corresponds to the loss of quantum information without loss of energy and amplitude damping described by (Eq. 12) which corresponds to spontaneous relaxation of qubits to |0 .
The relaxation probability p AD is defined by typical time scale T 1 such that p AD (t) = 1 − e −t/T1 , t being the idle time. Similarly, phase damping is associated to an error probability p P D with a typical time scale Moreover, the gate imperfections are modelled by a depolarizing channel acting after each gate [47]. The one-qubit channel is described by Eq. (13) (the two-qubit version being the tensor product of Eq. (13)). As most benchmark specifications of quantum devices provides only errors rate for 1-qubit gate 1 and 2-qubit gate 2 , we approximate the depolarizing probabilities by p 1 = 1 and p 2 = 2 . In Eq. (13), {X, Y, Z} are the three Pauli matrices:

D. Post-selection
To mitigate errors coming from noise, we post-select the measurements that respect the number of particles and S z symmetries. To do this, we first need to choose a measurement scheme to evaluate the Hamiltonian. We use the following scheme proposed in [21]: • The ZZ terms are measured in the computational basis without need for change of basis circuit.
• To efficiently measure the hopping terms 1 2 (XX + Y Y ), we implement a change of basis B that diagonalize the operator in the computational basis such that B † 1 2 (XX +Y Y )B = D, D being diagonal. The circuit shown in Fig. 14 transforms a hopping term into D = |01 01| − |10 10|.
As B also conserves the number of particles for each spin, we can reject any shots that don't contain the right number of spin up and down electrons, as they are only produced by errors occurring in the device. For example, for a 2-site model, any measurement (when evaluating the hopping terms or the onsite interaction terms) that is not contained in {0101, 1010, 1001, 0110} can be discarded. This method allows mitigating amplitude damping, read-out, and bit-flip errors.

E. Richardson extrapolation
This error mitigation method was introduced first in Ref. [19]. Let us suppose an open system composed of qubits initially prepared in ρ 0 . We apply a quantum circuit described by a time-dependent driving Hamiltonian K(t) = α J α (t)P α , where P α are Pauli strings and J α (t) real coefficients. The system will be described by Eq. (14): ∂ ∂t ρ(t) = −i[K(t), ρ(t)] + λL(ρ(t)).
L(ρ) is a noise generator that is invariant under time rescaling and independent from K(t), which can be experimentally quite demanding. Then, an observableÔ can be measured on the state ρ λ (T ) drived by K for a duration T , such that O K (λ) = tr(Ôρ λ (T )), which can be expressed for λ 1 (λ being the experimental noise strength) by the following: O * being the noise-free expectation value: O * = tr(Ôρ 0 (T )). Let assume one can scale up the noise level λ by a factor c. By running the quantum circuit for different noise levels λ i = c i λ for i = 0, ..., n (c 0 = 1, c i > 1 for i > 0) one can obtain an improved estimate of O * up to precision O(λ n+1 ) expressed as: with the coefficients {γ i } being solutions of i γ i = 1 and i γ i c k i = 0 for k = 1, ..., n. In the case n = 1, this is equivalent to a simple linear extrapolation: