## Classical Benchmarks for Variational Quantum Eigensolver Simulations of the Hubbard Model

7

Antonios M. Alvertis, 1, 2 Abid Khan, 3 Thomas Iadecola, 4, 5 Peter P. Orth, 4, 5, 6 and Norm Tubman 1 KBR, Inc., NASA Ames Research Center, Moffett Field, CA 94035, United States 2 Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States 3 Department of Physics, University of Illinois Urbana-Champaign, Urbana, IL, United States 61801 4 Department of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA 5 Ames National Laboratory, Ames, IA 50011, USA 6 Department of Physics, Saarland University, 66123 Saarbr¨cken, Germany u 7 NASA Ames Research Center, Moffett Field, CA 94035, United States

(Dated: August 5, 2024)

Simulating the Hubbard model is of great interest to a wide range of applications within condensed matter physics, however its solution on classical computers remains challenging in dimensions larger than one. The relative simplicity of this model, embodied by the sparseness of the Hamiltonian matrix, allows for its efficient implementation on quantum computers, and for its approximate solution using variational algorithms such as the variational quantum eigensolver. While these algorithms have been shown to reproduce the qualitative features of the Hubbard model, their quantitative accuracy in terms of producing true ground state energies and other properties, and the dependence of this accuracy on the system size and interaction strength, the choice of variational ansatz, and the degree of spatial inhomogeneity in the model, remains unknown. Here we present a rigorous classical benchmarking study, demonstrating the potential impact of these factors on the accuracy of the variational solution of the Hubbard model on quantum hardware. We find that even when using the most accurate wavefunction ans¨tze for the Hubbard model, the error in its ground a state energy and wavefunction plateaus for larger lattices, while stronger electronic correlations magnify this issue. Concurrently, spatially inhomogeneous parameters and the presence of off-site Coulomb interactions only have a small effect on the accuracy of the computed ground state energies. Our study highlights the capabilities and limitations of current approaches for solving the Hubbard model on quantum hardware, and we discuss potential future avenues of research.

## I. INTRODUCTION

Understanding strongly-correlated electron systems is one of the central challenges of condensed matter physics. Materials with strong electronic interactions exhibit a wealth of interesting properties, including superconductivity [1, 2], Mott insulating behavior [3, 4], excitonic ground states [5, 6], non-Fermi liquid behavior [7], competing and intertwined orders [8-10], and non-trivial magnetism [11, 12]. One of the most famous models that is capable of capturing such phenomena is the FermiHubbard model [13-22]. In its canonical form, its Hamiltonian is written as servables. However, in recent years several approaches have been put forward in order to derive material-specific Hubbard models with additional complexity to that of Eq. (1), including ab initio downfolding [23-25], embedding techniques [26, 27] and more. The general framework of these schemes entails electronic structure calculations within, e.g. , density functional theory to derive the parameters entering the Hubbard model for a given material, allowing also for multiple orbitals per site, spatial anisotropy, long-range interactions, etc. , as a means of capturing the complexity of real materials [28, 29]. For example a recently proposed approach to ab initio downfolding generates the following Hamiltonian for a material [23]:

$$H = - t \sum _ { \sigma } \sum _ { \langle \mathbf R R ^ { \prime } \rangle } a _ { \mathbf R } ^ { \sigma \dagger } a _ { \mathbf R ^ { \prime } } ^ { \sigma } + U \sum _ { \mathbf R } n _ { \mathbf R \uparrow } n _ { \mathbf R \downarrow }, \quad ( 1 )$$

where R is the lattice vector corresponding to a site within the studied lattice, σ the electron spin, t a hopping integral between nearest neighbors ( ⟨·⟩ denotes nearestneighbor summation only), U a repulsive on-site interaction, a σ † R ( a σ R ) an operator creating (destroying) an electron of spin σ at site R , and n R ↑ = a ↑† R a ↑ R the number operator for spin-up electrons at site R , and similar for spin-down.

While the Hubbard model can qualitatively capture the behavior of strongly-correlated systems, it does not generally offer quantitative predictive accuracy of ob-

$$H & = \sum _ { \sigma } \sum _ { \mathrm R R ^ { \prime } } t _ { \mathrm R R ^ { \prime } } a _ { \mathrm R } ^ { \sigma \dagger } a _ { \mathrm R ^ { \prime } } ^ { \sigma } \\ & + \frac { 1 } { 2 } \sum _ { \sigma \rho } \sum _ { \mathrm R R ^ { \prime } } U _ { \mathrm R R ^ { \prime } } a _ { \mathrm R } ^ { \sigma \dagger } a _ { \mathrm R ^ { \prime } } ^ { \rho \dagger } a _ { \mathrm R ^ { \prime } } ^ { \rho } a _ { \mathrm R ^ { \prime } } ^ { \sigma },$$

̸

in the case of a single band. In this representation of a system, the hopping and Coulomb terms may vary across the lattice. This Hamiltonian represents the system within a low-energy subspace, typically around the Fermi level. Here R = R ′ represents on-site interactions, whereas R = R ′ refers to longer-range interactions. While here we will restrict ourselves to including nearest-neighbor terms, one could in principle account

for longer-range interactions as well.

Solving for the ground and excited states of the Hubbard model is a quantum many-body problem, which generally scales exponentially with the lattice size when solved exactly on classical computers. As a result, one is restricted to studying finite clusters of relatively small sizes [30], although there exist several approximate methods and numerical techniques which have allowed studying larger systems [31, 32]. On the other hand, the Hamiltonians of Eq. (1) and Eq. (2) are amenable to simulation on quantum computers, for example using variational quantum algorithms (VQAs), where a linear scaling with system size is possible in the ideal scenario [33]. Moreover, several other potential benefits of solving the Hubbard model on quantum computers, as these technologies progress into the fault-tolerant era, have been discussed [34]. These facts make quantum computing an attractive avenue to pursue for tackling such problems. Indeed, there have been several studies solving different flavors of the Hubbard model on quantum computers or classical simulators thereof [35-43]. Some of these studies have shown that VQAs can reproduce qualitative features of the Hubbard model such as its magnetic properties [39, 40], and it has been shown that for a 1 × 8 Hubbard chain increasing the ansatz complexity can lead to well-converged values for energies and site occupancies, while spin-spin correlation functions are more challenging to capture [40].

At this point it remains unclear whether the qualitative success of VQAs in describing features of the Hubbard model will also translate into the quantitative reproduction of ground state properties such as the energy, and whether this will depend on the system size and other factors. It is therefore important to systematically benchmark the impact of varying lattice sizes on the quantitative features of the solutions obtained with VQAs, particularly given the importance of extrapolations to the thermodynamic limit for drawing comparisons to established methods [31]. Moreover, it is crucial to establish the impact on VQA simulation accuracy of electronic correlations, the chosen wavefunction ansatz, and the presence of off-site Coulomb interactions and spatially inhomogeneous parameters.

Here we present a detailed classical benchmarking study of the variational quantum eigensolver [44, 45] (VQE), a promising algorithm for quantum simulation on noisy intermediate-scale quantum (NISQ) hardware, applied to single-band Hubbard models of varying degrees of complexity. Our focus lies on comparing the quality of different wavefunction ans¨tze when scaling to a larger system sizes. We perform classical tensor-network simulations of VQE calculations for Hubbard models at half filling on 1D and 2D lattices of up to 16 sites, corresponding to 32 qubits. We consider different onsite interaction strengths U/t = 2 8 and the inclusion of nearest-, neighbor interactions V and spatial inhomogeneities. We show that the so-called number preserving (NP) ansatz, designed specifically for Hubbard models [35], outper- forms other popular ans¨tze, but still incurs errors and a converges extremely slowly for larger lattices with strong electronic correlations. Moreover, we demonstrate that VQE with the NP ansatz can successfully describe Hubbard models with nonuniform parameters and interactions beyond the on-site term, which is encouraging for the utility of this approach for solving Hamiltonians representing complex materials. Additionally, we show that an optimization based on maximizing the overlap with a wavefunction computed with a classical reference method can result in a small reduction to the energy error of VQE in some cases, but it can substantially improve the wavefunction fidelity and capture long-range correlations. Our study serves as a reference point along the path to accurate simulation of quantum materials on near-term quantum computers.

## II. METHODS

## A. Classical simulation of VQE

Our aim is to obtain the ground state of generalized Hubbard models with Hamiltonians of the form in Eq. (2). Here, we focus on models with and without nearest-neighbor repulsion and spatially anisotropic hopping parameters. To do so, we simulate the VQE classically by representing the variational ansatz state as a matrix product state (MPS) within a recently proposed variational tensor network eigensolver (VTNE) approach [46]. Specifically, following Ref. [46], we start from a product state | ψ 0 ⟩ in a checkerboard configuration, i.e. , with alternating spin up and down. We then generate a variational ansatz state via a parameterized quantum circuit (PQC) as follows:

$$\left | \psi _ { \text{PGQC} } ( \theta ) \right \rangle = U _ { n } ( \theta _ { n } ) \dots U _ { 1 } ( \theta _ { 1 } ) \left | \psi _ { 0 } \right \rangle. \quad \ \ ( 3 )$$

The precise form of the operators U i ( i = 1 , . . . , n ) is determined by the choice of ansatz used in our simulations, more details on which are given below. Each of these operators takes as arguments a set of parameters θ i , which are initialized randomly. Some ans¨ atze also include additional parameter-free gates as specified below. The PQC is represented as an MPS | ψ χ ( θ ) ⟩ with bond dimension χ , and we can therefore compute the energy expectation value

$$E _ { \chi } ( \theta ) = \left \langle \psi _ { \chi } ( \theta ) \right | H \left | \psi _ { \chi } ( \theta ) \right \rangle, \text{ \quad \ \ } ( 4 )$$

with the Hamiltonian H represented as a matrix product operator (MPO). Within our optimization scheme, we vary the parameters θ of the PQC in order to minimize the energy in Eq. (4). Unless otherwise explicitly stated, we use the maximum bond dimension χ max = 2 n / q 2 , where n q = 2 N N x y is the number of qubits needed to simulate the Hubbard model on an N x × N y square lattice. Since this bond dimension is sufficient to exactly

represent an arbitrary wavefunction on n q qubits, the MPS representation yields the energy of the exact PQC. For all tensor operations in this work we have used the ITensor software package [47, 48]. Our VQE results are benchmarked against those obtained within the density matrix renormalization group (DMRG) [49] solution of the studied systems, as implemented within ITensor.

Similar to Ref. [46], we first optimize the noninteracting ( U = 0) case, and then use the resulting parameters θ to initialize the optimization of the full interacting case. The random initial parameters θ are obtained from a Gaussian distribution N (0 10 , -5 ) with zero mean and variance σ 2 = 10 -5 . The energy minimization is terminated when one of three conditions is satisfied: the energy tolerance [50] reaches 10 -7 , the energy gradient reaches 10 -6 , or the optimization reaches 1000 steps. For every value of U/t and for each lattice size, we perform ten independent optimizations (using the L-BFGS method [51]), starting from different random parameters θ , and we take the minimum value among these as our estimate of the ground state energy. In this manner, we reduce the chances of the system becoming trapped in a local minimum.

For some cases studied here, we compare the results from this energy-based optimization to those from an overlap-based optimization which minimizes the following loss function, which is defined as the logarithm of the infidelity of the VTNE wavefunction | Ψ VTNE ⟩ , with respect to the DMRG wavefunction | Ψ DMRG ⟩ , which we consider to be the ground truth:

$$f = \log _ { 1 0 } ( 1 - | \left < \Psi _ { V T N E } | \Psi _ { D M R G } \right > | ^ { 2 } ). \quad ( 5 ) \quad ^ { \text{ters} } _ { \text{inal} }$$

Comparing the performance of the energy- and overlapbased approaches helps to disentangle the effects of cost function landscape and ansatz expressivity in capturing the physics of the Hubbard model. Of course the ultimate goal of VQAs is to obtain the ground state properties of systems for which classical methods would struggle, which would in turn constitute such an overlap-based optimization impractical. However, even in those cases one could perform an overlap-based optimization using an approximate, but still reasonably accurate classicallycomputed ground state for the system of interest. This would in turn provide a good starting point for subsequent optimization based on energy minimization, making it important to benchmark the performance of this overlap-based approach.

Before we proceed, it is worth emphasizing that, while the VTNE approach is generally meant as a preoptimization for VQE [46], our simulations use the maximal bond dimension χ max (unless otherwise stated, such as for some of our 32 qubit results) and therefore constitute an exact simulation of the VQE algorithm for a given PQC. We will therefore simply refer to our VTNE optimization as VQE.

## B. Ans¨ atze

Below we outline the ans¨tze for the PQC in Eq. (3) a which are used throughout this manuscript.

## 1. Number preserving (NP) ansatz

The number preserving (NP) ansatz [35] is a generalization of the Hamiltonian Variational Ansatz (HVA) [52]. As discussed in Ref. [35], for a 2D squarelattice system, a layer of the HVA applies a unitary operator of the form

$$e ^ { i t _ { v _ { 2 } } H _ { v _ { 2 } } } e ^ { i t _ { h _ { 2 } } H _ { h _ { 2 } } } e ^ { i t _ { v _ { 1 } } H _ { v _ { 1 } } } e ^ { i t _ { h _ { 1 } } H _ { h _ { 1 } } } e ^ { i t _ { H _ { o } } H _ { o } }, \quad ( 6 )$$

where t H o , t h 1 , t h 2 , t v 1 , t v 2 are variational parameters. Additionally, H o consists of the on-site terms of the Hamiltonian. The vertical hopping terms in the Hamiltonian are partitioned into components H v 1 , H v 2 such that each component consists of mutually commuting terms, and the same is done with the horizontal hopping terms, where the partitions are denoted as H h 1 , H h 2 . The NP ansatz replaces each hopping and on-site term within H o , H v 1 , H v 2 , H h 1 , and H h 2 , with a general two-qubit NP operator, which takes two parameters θ, ϕ :

$$U _ { \text{NP} } ( \theta, \phi ) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos \theta & i \sin \theta & 0 \\ 0 & i \sin \theta & \cos \theta & 0 \\ 0 & 0 & 0 & e ^ { i \phi } \end{pmatrix}. \quad \quad ( 7 )$$

This generalized HVA introduces independent parameters for each term in the Hamiltonian, unlike the original HVA [52], which introduces a single parameter for each class of terms, e.g. , for all horizontal hopping terms. Therefore, while for the NP ansatz used here the number of parameters in the optimization scales with the system size, this is not the case in the original HVA. Since sites that are adjacent in the 2D lattice system are not necessarily adjacent in the 1D qubit encoding, a series of fermionic SWAPs are needed to bring together these sites in the circuit before applying the two-qubit NP operator. After we apply the NP operator, we apply fermionic SWAPs to send the qubits back. As discussed in Ref. [35], a simplified circuit has hopping terms between vertically adjacent qubits implemented locally using an NP gate by applying a parameterless fermionic SWAP gate

$$F S W A P = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & - 1 \end{pmatrix}.$$

A single layer of this ansatz consists of applying the U NP gate to each pair of qubits coupled by the Hamiltonian of Eq. (1), and requires 10 N N x y -4 N x -4 N y parameters for an N x × N y lattice. Note that the number of parameters in the NP ansatz scales to leading order as O ( ℓN N x y ),

where ℓ is the number of HVA layers. In contrast, the number of parameters scales as O ( ℓ ) in the original HVA of Ref. [33] (see Eq. (6)). Similar to Ref. [46], we apply R z ( θ ) gates to each qubit prior to the application of the NP ansatz, which was found to improve optimization.

## 2. Excitation preserving (EP) ansatz

We also employ an excitation preserving (EP) ansatz, which is very similar to the NP ansatz outlined above. As implemented within Qiskit [53], we apply the two-qubit gate

$$U _ { \text{EP} } ( \theta, \phi ) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos ( \theta / 2 ) & - i \sin ( \theta / 2 ) & 0 \\ 0 & - i \sin ( \theta / 2 ) & \cos ( \theta / 2 ) & 0 \\ 0 & 0 & 0 & e ^ { - i \phi } \end{pmatrix}. \quad ( 9 ) & \kappa ^ { \cdot }$$

This is very similar to the NP operator of Eq. (7). However, unlike the NP case, we do not perform fermionic swaps between the qubits, but rather apply the gate to neighboring qubits, where for two-dimensional lattices we adopt a one-dimensional indexing scheme, with the index of a site ( x, y ) given by ( x -1) N y + y . Therefore, for one-dimensional lattices, the NP and EP ans¨tze are a essentially identical. However, while for two-dimensional lattices the lack of fermionic swaps makes the EP ansatz more straightforward to use than the NP, it also means that vertically adjacent qubits are not coupled directly. As a result, despite using two-qubit gates which are closely related to those of the NP ansatz, the EP ansatz is less tailored to the the physics of the Hubbard model Eq. (1), as will also be reflected by our results on twodimensional lattices.

## 3. Unitary coupled cluster (UCC) ansatz

The unitary coupled cluster (UCC) ansatz consists of applying excitation operators to the initial wavefunction | ψ 0 ⟩ as follows:

$$\left | \Psi _ { \text{UCC} } \right \rangle = \exp \left ( \hat { T } - \hat { T } ^ { \dagger } \right ) \left | \Psi _ { 0 } \right \rangle, \quad \ \ ( 1 0 ) \quad \text{the}$$

$$\hat { T } = \sum _ { i } ^ { \text{occ} } \sum _ { a } ^ { \text{vir} } \theta _ { i } ^ { a } \hat { a } _ { a } ^ { \dagger } \hat { a } _ { i } + \sum _ { i j } ^ { \text{occ} } \sum _ { a b } ^ { \text{vir} } \theta _ { i j } ^ { a b } \hat { a } _ { a } ^ { \dagger } \hat { a } _ { b } ^ { \dagger } \hat { a } _ { j } \hat { a } _ { i } + \cdots \. \ \ ( 1 1 ) \quad \ \text{then} \quad \ \text{begin} { } \\ \text{region} \\ \text{we nl.}$$

Here the creation and annihilation operators ˆ a † and ˆ act a on the occupied ( i, j, . . . ) and virtual ( a, b, . . . ) orbitals of the initial wavefunction, respectively.

We employ the factorized form of the UCC ansatz,

$$| \Psi _ { U C C } \rangle = \prod _ { i j \cdots \ a b \cdots } ^ { \text{occ} } \prod _ { \ } \hat { U } _ { i j \cdots } ^ { \text{ab} \cdots } | \Psi _ { 0 } \rangle \,, \quad \quad ( 1 2 ) \quad \text{gen} \quad \text{and}$$

for which the individual gates are defined as

$$\hat { U } _ { i j \cdots } ^ { a b \cdots } = \exp \left ( \theta _ { i j \cdots } ^ { a b \cdots } ( \hat { a } _ { i j \cdots } ^ { a b \cdots } - \hat { a } _ { a b \cdots } ^ { i j \cdots } ) \right ) \\ \hat { a } _ { i j \cdots } ^ { a b \cdots } = \hat { a } _ { a } ^ { \dagger } \hat { a } _ { b } ^ { \dagger } \dots \hat { a } _ { j } \hat { a } _ { i }.$$

We only include single- and double-excitation operators, ˆ U a i and ˆ U ab ij , within the so-called unitary coupled cluster ansatz with singles and doubles (UCCSD). The UCCSD ansatz, which we will refer to simply as the UCC ansatz, is number preserving, similar to the NP and EP ans¨tze, a however, unlike the linear scaling of the number of parameters of those ans¨tze with system size (per layer), a UCCSD scales quadratically [54], becoming more complex for larger systems.

## III. RESULTS

We now systematically study the impact of varying the number of variational parameters and the ansatz, the strength of the electronic interactions, and the strength of disorder and off-site Coulomb interactions, on the quality of the VQE solutions of the Hubbard model with respect to DMRG. Moreover, we benchmark the performance of using an overlap-based optimization, compared to the more conventional energy-based optimization. The optimized VQE energies for all cases presented below are given in the Supplemetary Material.

## A. Impact of the ansatz and the number of variational parameters

We start our analysis by describing the VQE result for the Hubbard model of different system sizes in d = 1 and d = 2 square lattice geometries. We employ open boundary conditions throughout. We will first explore how the NP ansatz, which was designed specifically for solving the Hubbard model on quantum hardware, performs as a function of lattice size and number of ansatz parameters. Our benchmarks make use of PQCs of up to thirteen layers, with each layer contributing 10 N N x y -4 N x -4 N y parameters.

In Fig. 1a we plot for several one-dimensional lattices the difference between our VQE energy and the DMRG benchmark per lattice site, for U/t = 2. This interaction strength lies in the weakly- to moderately-correlated regime [32, 55]. As our ground state energy from VQE we plot our best estimate, obtained as the lowest value from ten independent optimizations with random starting points. The shaded regions in Fig. 1a indicate the full range of energies obtained over all optimizations. While it is clear that the VQE optimization converges to the ground state energy even for the larger lattices as more parameters are added to the PQC, the convergence becomes increasingly slow for larger lattice sizes, and even plateaus to errors of order 10 -2 in those cases.

a

FIG. 1. Energy difference per lattice site, between VQE and DMRG ground state energies when using the NP and EP ans¨tze, a for 1D and 2D lattices, on a logarithmic scale. The shaded region indicates the difference between the minimum and maximum ground state energies obtained using VQE with a total of ten different random starting points. The inset of panel b shows a log-log plot of the minimum number of parameters needed to achieve an error per lattice site ∆ = 0 01, when using the NP . ansatz for 2D lattices.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-fdc17f70-c3ed-4aaf-a850-fcf702f485e2.jpg)

For the 1 × 12 lattice with 13 layers of the NP ansatz we have also performed an extended optimization, where we performed a hundred independent optimizations (rather than ten) starting from different random parameters. Nevertheless, even when we explore a larger part of phase space, we do not find any appreciable reduction in the ground state energy produced by VQE. Fig. 1b visualizes the performance of VQE on a logarithmic scale for two-dimensional lattices at U/t = 2. We again see that the convergence of VQE using the NP ansatz slows down with increasing system size, consistent with our benchmarks of the one-dimensional case.

To further understand the convergence behavior of the VQE energy towards the DMRG value, we plot in the inset of Fig. 1b the minimum number of parameters that is required in order to reach an error value ∆ = E VQE -E DMRG N of 0 01 with the NP ansatz and for . 2D lattices. We see that the necessary number of parameters to reach ∆ = 0 01 scales as ( . N x × N y ) n , with n = 2 13. . A similar power-law behavior with an exponent n in the range of 1 6 . -2 is found for 1D systems depending on the value of ∆.

We now turn to results obtained with the EP ansatz.

FIG. 2. Energy difference per lattice site, of the VQE ground state from DMRG for a 2 × 3 lattice, as a function of the number of parameters in the NP and EP ans¨tze (corresponding a to 1 to 13 ansatz layers), and also when using the UCC ansatz.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-8db068bb-b370-4af8-aaa8-698b0b3e7954.jpg)

Figs. 1c and 1d visualize the error per lattice site of VQE with respect to DMRG for U/t = 2 and for one- and twodimensional lattices, respectively. As for the case of the NP ansatz, we optimize PQCs of up to thirteen layers.

For one-dimensional lattices, the performance of VQE based on the EP ansatz is very similar to that of VQE based on the NP ansatz. However, for two-dimensional lattices, it is immediately evident that the VQE energy converges to the DMRG one at a much slower rate than when using the NP ansatz. Moreover, for larger lattices we observe that the VQE energy plateaus at substantial error values; when using 13 layers of this ansatz for a 3 × 3 lattice for example, the error per lattice site of the VQE energy compared to the DMRG value is 0 019-. 0 029, compared to 0 0003-0 001 for the NP ansatz. . . . We attribute this behavior to the fact that the EP ansatz does not directly couple vertically adjacent qubits in twodimensional lattices, due to the lack of fermionic swaps, hence missing a key part of the physics of the Hubbard model in which entanglement is generated efficently and directly between neighboring lattice sites.

The UCC ansatz is a popular choice for VQE calculations, particularly for chemistry applications. Given this popularity, we benchmark in Fig. 2 the performance of the NP, EP and UCC ans¨ atze for the solution of a 2 × 3 Hubbard model with U/t = 2. Unlike the layered NP and EP ans¨ atze, the UCC ansatz contains a set number of oneand two-electron excitation operators, and a number of parameters that is determined only by the lattice size. We see in Fig. 2 that for our case study, the UCC ansatz requires significantly more parameters (quadratic scaling, see Section II B) compared to NP or EP PQCs of up to thirteen layers (linear scaling), thus leading to a higher computational cost. Indeed we were not able to study lattices larger than 2 × 3 using the UCC ansatz due to the steep increase in the number of optimization parameters - for a 3 × 3 lattice, the ansatz has 3213 parameters, more than three times that of the NP ansatz at thirteen layers. Moreover, we see in Fig. 2 that despite the larger number of optimization parameters associated with the UCC ansatz, there is no significant improvement of the results compared to the EP ansatz, while the NP ansatz significantly outperforms it. In the future, it will be worthwhile to explore ways to optimize the UCC ansatz for lattice model simulations such as by adaptive approaches or by truncating excitation operators through a ranking scheme of their relative importance, which has been implemented for molecular systems through an MP2 pre-optimization scheme [56]. Another striking observation when comparing the layered ans¨tze is that, for small parameter numbers a ≲ 150, the EP ansatz actually slightly outperforms NP. However the rate of convergence with the number of parameters slows down markedly for the EP ansatz as this parameter number threshold is surpassed, and the NP ansatz systematically yields smaller errors. We find that it is a general feature across different lattices that for smaller number of parameters the NP and EP ans¨tze perform similarly, a however the energy error obtained from optimizing the EP PQC plateaus for larger numbers of parameters.

a

FIG. 3. Energy difference per lattice site, of the VQE ground state from DMRG 1D (panel a ) and 2D (panel b ) lattices, when using the NP ansatz and for U/t = 8. In panel c we compare the energy error per lattice site of the VQE ground state between the cases of U/t = 2 (solid lines) and U/t = 8 (dashed lines) for 2D lattices.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-fe5750ae-0b80-47ff-887f-0bdfdf7b9096.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-5fb3db14-44da-4b16-989e-1d96e8a21385.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-3f0a8698-2e76-4a72-b95b-d33170ec2d56.jpg)

## B. Impact of interaction strength U/t

Up to this point, we have focused on the case with interaction strength U/t = 2. However, for several applications of interest in condensed matter physics, stronger electronic correlations become important, and previous results on 1 × 8 Hubbard chains indicate that this regime

FIG. 4. Comparison of the energy difference per lattice site, of the VQE ground state from DMRG for 2D lattices, when using the NP ansatz and for U/t = 2, in the presence of disorder d (panel a ) and off-site Coulomb interactions V (panel b ). Solid lines represent the case of no disorder/off-site Coulomb interaction d = 0 , V = 0, dotted lines represent d = 0 2 and . V = 0 2 . respectively, while dashed lines represent d = 0 8 and . V = 0 8 respectively. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-efc91e39-4548-4761-8870-0952929082fa.jpg)

might be more challenging to capture with VQAs [40]. To address the quality of VQE simulations in this strong coupling regime, we set U/t = 8. Motivated by the results of Fig. 1, we adopt the NP ansatz in the following. Figs. 3a and 3b demonstrate the convergence of the VQE results towards the DMRG ones as a function of the number of parameters in the PQC for various lattice sizes in 1D and 2D, respectively. It becomes clear from these results that stronger electronic correlations pose more of a challenge to accurate VQE simulations of the Hubbard model, even when using the NP ansatz. To emphasize this, we plot in Fig. 3c the VQE energy error per site for 2D lattices, for U/t = 2 and U/t = 8. It is also worth noting that similar to U/t = 2 and Fig. 1b, we find for U/t = 8 that there is a power-law dependence of the minimum number of parameters needed to achieve an accuracy of ∆ = 0 01, on the lattice size, with an exponent . of n = 2 18 for 2D systems and . n = 1 61 for 1D systems. .

The challenge posed by systems with greater electronic correlation motivates a search for an ansatz that is tailored to the strong coupling regime, for example, one that exploits the well-known mapping of the Hubbard to the Heisenberg model at strong coupling. Another possibility to improve the ground state energies for strongly correlated Hubbard models could be to use solutions to the Heisenberg model as starting points for VQE, as was recently proposed [57].

## C. Impact of disorder and off-site Coulomb interactions

An important factor towards increasingly realistic materials simulations on quantum hardware is spatial inhomogeneity not present in the simple Hubbard model of Eq. (1). Indeed, the material-specific Hubbard models of the form (2), derived through methodologies such as ab initio downfolding [23-25] and embedding [26, 27], allow for spatial variations in the hopping and interaction parameters, e.g. , due to static disorder present in the material.

We therefore now consider the impact of such disorder on the accuracy of VQE simulations of Hubbard models. To do so, we study the ground state properties of the Hamiltonian

$$\frac { \text{a} } { \colon d } \quad H = - \sum _ { \sigma } \sum _ { \langle \text{R} \text{R} ^ { \prime } \rangle } t _ { \text{R} ^ { \prime } } a _ { \text{R} } ^ { \sigma \dagger } a _ { \text{R} ^ { \prime } } ^ { \sigma } + U \sum _ { \text{R} } n _ { \text{R} \uparrow } n _ { \text{R} \downarrow }, \quad ( 1 4 )$$

̸

where we now have disorder in the nearest-neighbor hopping ( R = R ′ ) and chemical potential ( R = R ′ ) terms, t RR ′ = t · δ ⟨ RR ′ ⟩ + δt RR ′ . The function δ ⟨ RR ′ ⟩ is equal to one for nearest neighbors, and zero otherwise, since the chemical potential in our system is set to zero. The disorder enters via δt RR ′ = d · N (0 1), with , d a user-defined parameter for the magnitude of the disorder and N (0 , 1) a function drawing a random number from a Gaussian distribution of zero mean and unit variance. In this way we ensure that the values of t RR ′ of nearest neighbor hopping are centered around t , and the chemical potential is centered around zero, both with a standard deviation of d . In Fig. 4a we plot the error per lattice site of the VQE ground state energies with respect to the DMRG ones, for different lattices and as a function of the number of optimization parameters in the ansatz. We use disorder values of d = 0 2 and . d = 0 8. . We observe a small increase in the error of the VQE energies as compared to the homogeneous case ( d = 0), however this increase is not monotonic with the value of d . Despite this small increase in the VQE energy error, these results suggest that even for disordered systems with nonuniform parameters, VQE using the NP ansatz can recover the ground state of the Hubbard model with reasonable accuracy, within different regimes of disorder.

So far we have focused on Hubbard models where Coulomb interactions are present only on-site. However,

## Energy- vs overlap-based optimization

a

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-415233e9-3cf3-440a-b5c7-ecd5886abf01.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-5b511bfe-2076-4b08-a2e8-2d56204a5724.jpg)

c

FIG. 5. Comparison of the energy difference per lattice site, of the VQE ground state from DMRG, as obtained from an energy-based optimization (solid lines) and an overlap-based optimization (dashed lines), when using the NP ansatz. The different panels visualize the results for 1D and 2D lattices, and for the cases U/t = 2 and U/t = 8.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-c34cf753-cfb3-4c60-b130-e62f153b3b69.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-3ae91ce3-aea8-4403-bd8a-fd396dbc71b1.jpg)

in realistic material systems there are also longer-range interactions present, which motivates us to investigate the impact of including a finite nearest-neighbor repulsive Coulomb interaction V in our model, i.e. , to consider models of the form

$$H = - t \sum _ { \sigma } \sum _ { \langle \mathbf R R ^ { \prime } \rangle } a _ { \mathbf R } ^ { \sigma \dagger } a _ { \mathbf R ^ { \prime } } ^ { \sigma } + U \sum _ { \mathbf R } n _ { \mathbf R \uparrow } n _ { \mathbf R \downarrow } \quad \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } \\ + V \sum _ { \sigma \sigma ^ { \prime } } \sum _ { \langle \mathbf R \mathbf R ^ { \prime } \rangle } n _ { \mathbf R \sigma } n _ { \mathbf R ^ { \prime } \sigma ^ { \prime } }. \quad \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{\ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{\ } \\ & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{  & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & \text{$$

In Fig. 4b we visualize the energy error per lattice site of the VQE solutions of the Hubbard model on 2D lattices, with respect to the DMRG energy, as a function of the number of parameters used in the optimization, and for different values of V , specifically V = 0 (solid lines), V = 0 2 (dotted lines) and . V = 0 8 (dashed lines). . Similar to the case of Fig. 4a of having finite disorder in the lattice, we see that while indeed having off-site interactions generally leads to slightly larger errors in the VQE energy, the increase of the error is non-monotonic with V . The fact that VQE with the NP ansatz seems capable of accurately describing spatially inhomogeneous systems with Coulomb interactions beyond on-site is encouraging for the simulation of increasingly complex and inhomogeneous systems on quantum hardware. We relate this capability to the fact that, unlike in the original HVA ansatz [33], the NP ansatz contains independent variational parameters on every link that appears in the Hamiltonian.

## D. Impact of using an overlap-based optimization

Here we compare the performance of conventional VQE, where we drive the optimization by minimizing the energy expectation value of the wavefunction, to overlapbased optimization, where we minimize the logarithm of the infidelity with respect to the DMRG ground state (eq. (5)). In Fig. 5 we visualize the energy error per lattice site of VQE based on overlap optimization (dashed lines), alongside the values previously obtained from energy optimization (solid lines), for 1D and 2D lattices, and for the cases of U/t = 2 and U/t = 8. We find that in the one-dimensional case the overlap-based optimization

U/t=2 -2D

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-f460e208-5fe4-4757-aa54-abdb4319e737.jpg)

U/t=8 -2D

FIG. 6. Comparison of the fidelity of the VQE ground state with respect to DMRG, as obtained from an energy-based optimization (solid lines) and an overlap-based optimization (dashed lines), when using the NP ansatz, and for 2D lattices with U/t = 2 (panel a ) and U/t = 8 (panel b ).

leads to a small improvement of the VQE energies, for both U/t = 2 and U/t = 8. This suggests that the residual energy error for larger lattices is a result of the complex cost function landscape rather than the expressivity of the NP ansatz. However, for two-dimensional lattices we find no substantial differences between the two methods, and the energy-based optimization can even slightly outperform the overlap-based one in some cases.

It is also worth investigating how the energyand overlap-based optimization perform in terms of producing high-fidelity VQE ground states. Figs. 6a and 6b visualize the fidelity F = | ⟨ Ψ VQE | Ψ DMRG ⟩ | 2 of the VQE ground state, when using the NP ansatz for 2D lattices, and for U/t = 2 and U/t = 8. We see that with the addition of more parameters to the optimization, both methods converge towards the limit of F = 1, except for the largest system size 3 × 3. However, the convergence of the overlap-based optimization is generally faster, especially for larger lattices, a feature which is particularly prominent in the case of U/t = 8, where strong electronic correlations are present. The fact that the overlapbased optimization performs better in terms of yielding higher-fidelity ground states is not surprising, given that this method is tailored specifically to minimize the loss function of eq. (5), i.e. , to minimize the infidelity of the VQE ground state. However, it is worth emphasizing that the higher fidelity solutions produced via overlapbased optimization do not always possess lower energies compared to the lower-fidelity wavefunctions obtained via energy-based minimization. An example of this is seen for 3 × 3 lattices, and for both U/t = 2 and U/t = 8. While overlap-based optimization produces ground states of higher fidelities as seen in Fig. 6, the energy error of these states is in fact larger than that from energy-based optimization, as seen in Fig. 5.

FIG. 7. Converged energy difference of VQE and DMRG energy per lattice site, for a 4 × 4 Hubbard model with U/t = 2, and for different starting points of the optimization. Both levels of theory employ a bond dimension of χ = 512. Here we use 12 layers of the NP ansatz, corresponding to 1,568 optimization parameters within VQE.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-4b59c8cf-5d48-4002-889b-5d3a42670349.jpg)

## E. Scaling up and extracting correlation functions

Our results so far have demonstrated that convergence of the VQE energy to the DMRG value with respect to the number of the optimization parameters becomes slower for larger lattice sizes. It is therefore interesting to explore this convergence in more detail for a larger lattice of size 4 × 4 for U/t = 2, to examine whether this trend holds even for larger systems. Since each lattice site can be occupied by a spin-up and/or a spin-down electron, this corresponds to 4 × × 4 2 = 32 qubits. Due to the large computational cost associated with simulating a lattice of this size, we restrict the bond dimension of the MPS representation within our VTNE approach, and also within

&gt;

&lt;c1j

&gt;

&lt;c1j

FIG. 8. Spin correlation function ⟨ C 1 j ⟩ (panel a ) and ⟨ C ij ⟩ (panel b at the DMRG level, panel c within energy-based VQE, and panel d within overlap-based VQE), for a 1 × 6 Hubbard model with U/t = 8. Panels e h -visualize the same quantities for a 1 × 12 Hubbard chain instead.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00836v1-figure-8b90a9a1-f4a1-4218-aaaa-c000f313d49e.jpg)

DMRG, to χ = 512. Here we also set d = V = 0. As in previous cases, we initialize ten VQE optimizations from different, random starting points. Two of these failed to converge, and for the remaining eight we visualize their converged energy errors per lattice site with respect to DMRG in Fig. 7. The lowest energy obtained through VQE optimization, with 12 NP ansatz layers corresponding to 1,568 parameters, has a relative error of 2 2% in . its energy, compared to the DMRG value. This demonstrates that even for these larger lattices, reasonable accuracy for the ground state of the Hubbard model can be achieved using VQE, as long as enough parameters are included in the optimization.

Overall, our work so far has focused on the quantitative benchmarking of VQE simulations of the Hubbard model. However, qualitative features of the model can be important to reproduce, and indeed VQAs have been shown in various examples to yield ground states that display the well-known antiferromagnetic behavior of the Hubbard model in the strong coupling regime [31, 40]. Our VQE implementation based on tensor networks can recover such behaviors; in Fig. 8, we visualize the spin correlation function

$$\langle C _ { i j } \rangle = \langle S _ { i } ^ { z } S _ { j } ^ { z } \rangle - \langle S _ { i } ^ { z } \rangle \langle S _ { j } ^ { z } \rangle \quad \quad ( 1 6 ) \quad \sin _ { - 1 _ { - } }$$

in the case of U/t = 8 for a 1 × 6 and 1 × 12 chain, within VQE based on energy- and overlap-based optimizations, and within DMRG. To obtain the VQE results here we have used 13 layers of the NP ansatz. We observe that for the smaller chain, the spin properties from both VQE approaches and DMRG are in near perfect agreement. For the larger 1 × 12 chain, VQE based on an energy minimization fails to fully capture longer-range spin correlations, in agreement with previous findings [39, 40]. However, the VQE ground state obtained through overlapbased optimization is in near-perfect agreement to the DMRG one in terms of its long-range spin correlation, highlighting the utility of this approach for computing such features.

## IV. CONCLUSIONS AND OUTLOOK

In this work we have presented a rigorous benchmark of the quantitative features of the solution of the Hubbard model using variational quantum eigensolvers. We find that the NP ansatz of Ref. [35] yields energies with the smallest error relative to a DMRG reference, compared to the EP and UCC ans¨tze. However, even when a using the NP ansatz, the error in the ground state energy increases for larger lattices, and ultimately plateaus with increasing number of optimization parameters. Moreover, we have found that systems with stronger electronic correlations pose an even greater challenge in terms of their accurate quantitative description, and require more optimization parameters. At the same time, spatial inhomogeneity in the parameters of the Hubbard model does not seem to strongly affect the convergence rate of the results, nor does the addition of nearest-neighbor Coulomb interactions. This is encouraging as we are moving towards the description of more complex versions of the Hubbard model, representing real materials through schemes such as ab initio downfolding. We have also shown, in agreement with previous work, that the qualitative characteristics of the ground state wavefunction of the one-dimensional Hubbard model reproduce the expected magnetic behavior depending on the U/t ratio. We found that while for smaller lattices the VQE solutions obtained via energy minimization can give quantitatively accurate results for spin correlation functions, for larger systems it struggles to capture long-range correlations. However, a VQE optimization based on maximizing the overlap to the solution of a classical reference method succeeds in describing these long-range phenomena. While of course ultimately it is the goal of

the community to utilize quantum computers to obtain the properties of Hubbard model ground states for systems beyond the reach of classical methods, where computing the overlap to the 'true' ground state will not be possible, this result still highlights the potential of an overlap-based optimization using an approximate classical reference state to provide a strong starting point for energy-based VQE. The overlap-based optimization results also highlight that the NP ansatz is expressive enough to capture features like long-range correlations, suggesting that improvements to the energy-based optimization approach may be able to remedy these issues.

Moving forward, it will be important to extend the applicability of current approaches to the solution of Hubbard models with multiple electronic bands, as these multi-band models describe several physical phenomena of relevance to applications, including superconductivity and excitonic insulating behavior [19, 58-60]. Given the challenge with systematically converging to the global energy minimum with current variational approaches as demonstrated here, even in the one-band case and particularly for the strongly correlated regime, it will be important to develop flexible approaches to achieve this. One

- [1] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, Strongly Correlated Superconductivity, Science 296 , 2364 (2002).
- [2] D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Superconductivity in an infinite-layer nickelate, Nature 572 , 624 (2019).
- [3] D. Lee, B. Chung, Y. Shi, G.-Y. Kim, N. Campbell, F. Xue, K. Song, S.-Y. Choi, J. P. Podkaminer, T. H. Kim, P. J. Ryan, J.-W. Kim, T. R. Paudel, J.-H. Kang, J. W. Spinuzzi, D. A. Tenne, E. Y. Tsymbal, M. S. Rzchowski, L. Q. Chen, J. Lee, and C. B. Eom, Isostructural metal-insulator transition in VO , 2 Science 362 , 1037 (2018).
- [4] S. Grytsiuk, M. I. Katsnelson, E. G. C. P. van Loon, and M. R¨ osner, Nb3Cl8: a prototypical layered MottHubbard insulator, npj Quantum Materials 9 , 8 (2024).
- [5] Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. J¨ck, a S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak, N. Fishchenko, X. Liu, G. Farahi, F. Xie, Y. Xu, K. Watanabe, T. Taniguchi, B. A. Bernevig, R. J. Cava, L. M. Schoop, A. Yazdani, and S. Wu, Evidence for a monolayer excitonic insulator, Nature Physics 18 , 87 (2022).
- [6] L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak, and J. Shan, Strongly correlated excitonic insulator in atomic double layers, Nature 598 , 585 (2021).
- [7] G. R. Stewart, Non-fermi-liquid behavior in d -and f -electron metals, Rev. Mod. Phys. 73 , 797 (2001).
- [8] E. Dagotto, Complexity in strongly correlated electronic systems, Science 309 , 257 (2005), https://www.science.org/doi/pdf/10.1126/science.1107559.
- [9] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Collo-

possibility is to explore other ansatz generation schemes, such as adaptive ans¨tze [43, 61-65] or problem-inspired a ans¨tze tailored for the strongly-correlated regime. a

## V. ACKNOWLEDGMENTS

The authors acknowledge useful discussions with Bryan Clark and Erik Gustafson. This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Superconducting Quantum Materials and Systems Center (SQMS) under contract No. DE-AC02-07CH11359. We are grateful for support from NASA Ames Research Center. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award ASCR-ERCAP0024469.

- quium: Theory of intertwined orders in high temperature superconductors, Rev. Mod. Phys. 87 , 457 (2015).
- [10] S. Sachdev, Colloquium: Order and quantum phase transitions in the cuprate superconductors, Rev. Mod. Phys. 75 , 913 (2003).
- [11] H. L. Zhuang and R. G. Hennig, Stability and magnetism of strongly correlated single-layer vs 2 , Phys. Rev. B 93 , 054429 (2016).
- [12] H. Guo, Z. W. Li, L. Zhao, Z. Hu, C. F. Chang, C.-Y. Kuo, W. Schmidt, A. Piovano, T. W. Pi, O. Sobolev, D. I. Khomskii, L. H. Tjeng, and A. C. Komarek, Antiferromagnetic correlations in the metallic strongly correlated transition metal oxide LaNiO3, Nature Communications 9 , 43 (2018).
- [13] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, The Hubbard Model, Annual Review of Condensed Matter Physics 13 , 239 (2022).
- [14] H.-C. Jiang and T. P. Devereaux, Superconductivity in the doped Hubbard model and its interplay with nextnearest hopping t , Science ′ 365 , 1424 (2019).
- [15] E. Gull, O. Parcollet, and A. J. Millis, Superconductivity and the pseudogap in the two-dimensional hubbard model, Phys. Rev. Lett. 110 , 216405 (2013).
- [16] S. S. Kancharla and S. Okamoto, Band insulator to mott insulator transition in a bilayer hubbard model, Phys. Rev. B 75 , 193103 (2007).
- [17] H. Lee, Y.-Z. Zhang, H. O. Jeschke, and R. Valent´ ı, Competition between band and mott insulators in the bilayer hubbard model: A dynamical cluster approximation study, Phys. Rev. B 89 , 035139 (2014).
- [18] J. E. Hirsch and S. Tang, Antiferromagnetism in the twodimensional hubbard model, Phys. Rev. Lett. 62 , 591 (1989).
- [19] T. Kaneko, K. Seki, and Y. Ohta, Excitonic insulator

- state in the two-orbital hubbard model: Variational cluster approach, Phys. Rev. B 85 , 165135 (2012).
- [20] K. Sugimoto, S. Nishimoto, T. Kaneko, and Y. Ohta, Strong coupling nature of the excitonic insulator state in ta 2 nise 5 , Phys. Rev. Lett. 120 , 247602 (2018).
- [21] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68 , 13 (1996).
- [22] Q. Si and G. Kotliar, Fermi-liquid and non-fermi-liquid phases of an extended hubbard model in infinite dimensions, Phys. Rev. Lett. 70 , 3143 (1993).
- [23] K. Nakamura, Y. Yoshimoto, Y. Nomura, T. Tadano, M. Kawamura, T. Kosugi, K. Yoshimi, T. Misawa, and Y. Motoyama, RESPACK: An ab initio tool for derivation of effective low-energy model of material, Computer Physics Communications 261 , 107781 (2021), arXiv:2001.02351.
- [24] H. Zheng, H. J. Changlani, K. T. Williams, B. Busemeyer, and L. K. Wagner, From real materials to model Hamiltonians with density matrix downfolding, Frontiers in Physics 6 , 1 (2018), 1712.00477.
- [25] R. Arita, H. Ikeda, S. Sakai, and M.-T. Suzuki, Ab initio downfolding study of the iron-based ladder superconductor bafe 2 s 3 , Phys. Rev. B 92 , 054515 (2015).
- [26] C. Vorwerk, N. Sheng, M. Govoni, B. Huang, and G. Galli, Quantum embedding theories to simulate condensed systems on quantum computers, Nature Computational Science 2 , 424 (2022).
- [27] T. Zhu and G. K.-L. Chan, Ab initio full cell gw +DMFT for correlated materials, Phys. Rev. X 11 , 021006 (2021).
- [28] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Frequencydependent local interactions and low-energy effective models from electronic structure calculations, Phys. Rev. B 70 , 195104 (2004).
- [29] P. R. C. Kent and G. Kotliar, Toward a predictive theory of correlated materials, Science 361 , 348 (2018).
- [30] E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc, and J. Riera, Static and dynamical properties of doped Hubbard clusters, Phys. Rev. B 45 , 10741 (1992).
- [31] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jim´ enez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof'ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5 , 041041 (2015).
- [32] M. Qin, T. Sch¨ afer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard Model: A Computational Perspective, Annual Review of Condensed Matter Physics 13 , 275 (2022).
- [33] D. Wecker, M. B. Hastings, N. Wiebe, B. K. Clark, C. Nayak, and M. Troyer, Solving strongly correlated electron models on a quantum computer, Phys. Rev. A 92 , 062318 (2015).
- [34] A. A. Agrawal, J. Job, T. L. Wilson, S. N. Saadatmand, M. J. Hodson, J. Y. Mutus, A. Caesura, P. D. Johnson, J. E. Elenewski, K. J. Morrell, and A. F. Kemper, Quantifying fault tolerant simulation of strongly correlated systems using the fermi-hubbard model (2024),
- arXiv:2406.06511 [quant-ph].
- [35] C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Strategies for solving the Fermi-Hubbard model on nearterm quantum computers, Physical Review B 102 , 1 (2020), arXiv:1912.06007.
- [36] P.-L. Dallaire-Demers and F. K. Wilhelm, Method to efficiently simulate the thermodynamic properties of the fermi-hubbard model on a quantum computer, Phys. Rev. A 93 , 032303 (2016).
- [37] J. M. Reiner, F. Wilhelm-Mauch, G. Sch¨ on, and M. Marthaler, Finding the ground state of the Hubbard model by variational methods on a quantum computer with gate errors, Quantum Science and Technology 4 , 10.1088/2058-9565/ab1e85 (2019), 1811.04476.
- [38] Z. Cai, Resource Estimation for Quantum Variational Simulations of the Hubbard Model, Phys. Rev. Appl. 14 , 014059 (2020).
- [39] S. Stanisic, J. L. Bosse, F. M. Gambetta, R. A. Santos, W. Mruczkiewicz, T. E. O'Brien, E. Ostby, and A. Montanaro, Observing ground-state properties of the FermiHubbard model using a scalable algorithm on a quantum computer, Nature Communications 13 , 5743 (2022).
- [40] B. Anselme Martin, P. Simon, and M. J. Ranˇi´, Simulatc c ing strongly interacting Hubbard chains with the variational Hamiltonian ansatz on a quantum computer, Phys. Rev. Res. 4 , 023190 (2022).
- [41] L. Mineh and A. Montanaro, Solving the Hubbard model using density matrix embedding theory and the variational quantum eigensolver, Phys. Rev. B 105 , 125117 (2022).
- [42] Y. Yao, F. Zhang, C.-Z. Wang, K.-M. Ho, and P. P. Orth, Gutzwiller hybrid quantum-classical computing approach for correlated materials, Phys. Rev. Research 3 , 013184 (2021).
- [43] A. Mukherjee, N. F. Berthusen, J. C. Getelina, P. P. Orth, and Y.-X. Yao, Comparative study of adaptive variational quantum eigensolvers for multi-orbital impurity models, Commun Phys 6 , 1 (2023).
- [44] M. Cerezo, K. Sharma, A. Arrasmith, and P. J. Coles, Variational quantum state eigensolver, npj Quantum Information 8 , 113 (2022).
- [45] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The Variational Quantum Eigensolver: A review of methods and best practices, Physics Reports 986 , 1 (2022).
- [46] A. Khan, B. K. Clark, and N. M. Tubman, Preoptimizing variational quantum eigensolvers with tensor networks (2023), arXiv:2310.12965.
- [47] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
- [48] M. Fishman, S. R. White, and E. M. Stoudenmire, Codebase release 0.3 for ITensor, SciPost Phys. Codebases , 4 (2022).
- [49] F. Verstraete, T. Nishino, U. Schollw¨ck, M. C. Ba˜uls, o n G. K. Chan, and M. E. Stoudenmire, Density matrix renormalization group, 30 years on, Nature Reviews Physics 5 , 273 (2023).
- [50] The energy tolerance is defined by the absolute difference between the energy at the final step and the energy at the penultimate step.
- [51] J. Nocedal and S. J. Wright, Numerical optimization (Springer, 1999).

- [52] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Phys. Rev. A 92 , 042303 (2015).
- [53] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).
- [54] J. R. McClean, J. Romero, R. Babbush, and A. AspuruGuzik, The theory of variational hybrid quantumclassical algorithms, New Journal of Physics 18 , 023023 (2016).
- [55] M. Qin, H. Shi, and S. Zhang, Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method, Phys. Rev. B 94 , 085103 (2016).
- [56] J. W. Mullinax and N. M. Tubman, Large-scale sparse wavefunction circuit simulator for applications with the variational quantum eigensolver (2023), arXiv:2301.05726.
- [57] B. Murta and J. Fern´ndez-Rossier, From heisenberg to a hubbard: An initial state for the shallow quantum simulation of correlated electrons, Phys. Rev. B 109 , 035128 (2024).
- [58] J. Kuneˇ s and P. Augustinsk´, y Excitonic instability at the spin-state transition in the two-band hubbard model, Phys. Rev. B 89 , 115134 (2014).
- [59] C. De Franco, L. F. Tocchio, and F. Becca, Metalinsulator transitions, superconductivity, and magnetism
- in the two-band hubbard model, Phys. Rev. B 98 , 075117 (2018).
- [60] L.-H. Hu and C. Wu, Two-band model for magnetism and superconductivity in nickelates, Phys. Rev. Res. 1 , 032046 (2019).
- [61] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10 , 3007 (2019).
- [62] D. Claudino, J. Wright, A. J. McCaskey, and T. S. Humble, Benchmarking Adaptive Variational Quantum Eigensolvers, Front. Chem. 8 , 10.3389/fchem.2020.606863 (2020).
- [63] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, QubitADAPT-VQE: An Adaptive Algorithm for Constructing Hardware-Efficient Ans \ 'atze on a Quantum Processor, PRX Quantum 2 , 020310 (2021).
- [64] N. Gomes, A. Mukherjee, F. Zhang, T. Iadecola, C.-Z. Wang, K.-M. Ho, P. P. Orth, and Y.-X. Yao, Adaptive Variational Quantum Imaginary Time Evolution Approach for Ground State Preparation, Advanced Quantum Technologies 4 , 2100114 (2021).
- [65] J. S. Van Dyke, K. Shirali, G. S. Barron, N. J. Mayhall, E. Barnes, and S. E. Economou, Scaling adaptive quantum simulation algorithms via operator pool tiling (2023), arxiv:2206.14215 [cond-mat, physics:quant-ph].
Antonios M. Alvertis Abid Khan Thomas Iadecola Peter P. Orth Norm Tubman
Classical Benchmarks for Variational Quantum Eigensolver Simulations of the Hubbard Model
Simulating the Hubbard model is of great interest to a wide range of
applications within condensed matter physics, however its solution on classical
computers remains challenging in dimensions larger than one. The relative
simplicity of this model, embodied by the sparseness of the Hamiltonian matrix,
allows for its efficient implementation on quantum computers, and for its
approximate solution using variational algorithms such as the variational
quantum eigensolver. While these algorithms have been shown to reproduce the
qualitative features of the Hubbard model, their quantitative accuracy in terms
of producing true ground state energies and other properties, and the
dependence of this accuracy on the system size and interaction strength, the
choice of variational ansatz, and the degree of spatial inhomogeneity in the
model, remains unknown. Here we present a rigorous classical benchmarking
study, demonstrating the potential impact of these factors on the accuracy of
the variational solution of the Hubbard model on quantum hardware. We find that
even when using the most accurate wavefunction ans\"{a}tze for the Hubbard
model, the error in its ground state energy and wavefunction plateaus for
larger lattices, while stronger electronic correlations magnify this issue.
Concurrently, spatially inhomogeneous parameters and the presence of off-site
Coulomb interactions only have a small effect on the accuracy of the computed
ground state energies. Our study highlights the capabilities and limitations of
current approaches for solving the Hubbard model on quantum hardware, and we
discuss potential future avenues of research.