Gutzwiller Hybrid Quantum-Classical Computing Approach for Correlated Materials

Inspired by the fast-pace evolution of noisy intermediate-scale quantum (NISQ) computing technology, novel resource-efficient hybrid quantum-classical approaches are under active development to address grand scientific challenges faced by classical computations. Proof-of-principle applications of NISQ technology in quantum chemistry have been reported in solving ground state properties of small molecules. While several approaches have also been proposed to address the long-standing strongly correlated materials problem, a complete calculation of periodic correlated electron model systems on NISQ devices, as a crucial step forward, has not been demonstrated yet. In this paper we showcase the first self-consistent hybrid quantum-classical calculations of the periodic Anderson model on Rigetti quantum devices, within the Gutzwiller variational embedding theoretical framework. It maps the infinite lattice problem to a quasi-particle Hamiltonian coupled to a quantum many-body embedding system, which is solved on quantum devices to overcome the classical exponential scaling in complexity. We show that the Gutzwiller hybrid quantum-classical embedding (GQCE) framework describes very well the quantum phase transitions from Kondo insulator to metal and from metal to Mott insulator in the correlated electron lattice model, with critical parameters at the phase boundaries accurately determined. The GQCE simulation framework, equipped with a full arsenal of evolving quantum algorithms to solve ground state and dynamics of quantum chemistry or equivalently finite embedding systems, is a well-adapted approach toward resolving complicated emergent properties in correlated condensed matter by exploiting NISQ technologies.


I. INTRODUCTION
Quantum computing holds the promise to revolutionize modern high-performance computations by potentially providing exponential speedups compared to currently known classical algorithms for a variety of important problems such as Hamiltonian simulation [1,2]. Accurately predicting the properties of competing phases or simulating the dynamics of interacting quantum mechanical many-body systems directly addresses grand challenges in quantum chemistry and materials science, with far reaching consequences on the science, technology and health sectors of our digital society [3].
While not being fully fault-tolerant, the currently available noisy intermediate-scale quantum (NISQ) hardware [4] is still extremely powerful as recently demonstrated by the Google team [5]. As the number of coherent gate operations is limited, however, the development of resource efficient algorithms with sufficiently short quantum circuits is crucial in order to be able to tackle open scientific problems on NISQ devices. One such example is the variational quantum eigensolver (VQE) algorithm that addresses the eigenvalue problem [6,7]. It was successfully implemented on NISQ technology to yield the ground state energy of small molecules such as H 2 , HHe + , LiH and BeH 2 [6,[8][9][10]. The algorithm is a hybrid quantum-classical approach that combines a quantum computation of a suitable cost function such as the Hamiltonian with a classical optimization routine. Instead of adiabatic state preparation followed by quan- * Electronic address: ykent@iastate.edu tum phase estimation [11,12], which requires deep circuits, VQE employs a shallow variational circuit to evolve a chosen initial state into a target state. The objective cost function is then measured as a weighted sum of expectation values for associated Pauli terms. The variational parameters are classically optimized to minimize the cost [6,7]. Different forms of the variational circuit have been discussed in the literature, for example, based on the unitary coupled cluster ansatz (UCC) [7,13,14], or a trotterized adiabatic preparation ansatz [15]. Common to all of them is that the number of variational parameters rapidly increases with the number of orbitals, which makes the generally non-convex classical optimization problem increasingly difficult to solve. This is further complicated by the presence of noise on real NISQ devices. While directly applicable to small molecules, simulating infinite periodic quantum materials on NISQ devices therefore requires further algorithmic development.
Various proposals of efficient quantum algorithms for correlated materials have been discussed in the literature. For example, using an adiabatic VQE approach in a dual plane wave basis leads to particularly favorable scaling of the circuit depth and the number of qubits required for periodic systems [16]. This work proposed to simulate the Jellium model as a benchmark on near-term devices. Alternatively, there exists a long tradition in correlated materials theory to map infinite periodic systems onto effective impurity models. For example, dynamical mean-field theory (DMFT) has been established as the state-of-the-art classical computational approach for correlated materials [17][18][19][20]. In a pioneering work quantum algorithms based on adiabatic state preparation and phase estimation were introduced to solve the impu-2 rity Green function repeatedly, upon reaching the convergence with the local lattice Green function [21]. A hybrid quantum-classical approach based on a simplified twosite version of DMFT has also been proposed [22], and recently implemented using a generalized VQE method to find both the ground and excited states of the impurity model [23]. None of the algorithms proposed in these influential works have yet been fully demonstrated on a real NISQ device, since they require resources that are still beyond the current technology [24].
Here, we propose and demonstrate a new resourceefficient hybrid quantum-classical algorithm to simulate correlated materials on NISQ devices: the Gutzwiller quantum-classical embedding (GQCE) simulation framework. We have implemented the algorithm on Rigetti's quantum cloud service using PyQuil [25,26], and have performed the first fully self-consistent calculations of an infinite periodic correlated electron model on a quantum computer. Our GQCE approach is based on the powerful Gutzwiller variational embedding theory that is known to be able to capture many of the phenomena associated with strong local correlations such as Mott-Hubbard transitions [27]. When combined with density-functional theory (DFT), it has been shown to successfully address ground state properties of real correlated materials [27][28][29][30][31][32][33][34].
Similar to DMFT, the Gutzwiller embedding method maps the infinite interacting lattice model onto an effective impurity problem consisting of a cluster of correlated orbitals that are embedded in a self-consistent medium. In contrast to DMFT, which requires to solve for the fully frequency dependent impurity self-energy, the Gutzwiller theory uses only the ground state single-particle density matrix of the correlated impurity cluster. In practice, the approach amounts to finding a self-consistent solution of a set of coupled eigenvalue equations. The Gutzwiller embedding method is therefore ideally suited to be formulated as a hybrid quantum-classical algorithm, where the ground state of the correlated impurity cluster is determined using VQE. The GQCE approach is illustrated in Fig. 1.
As a first non-trivial benchmark study we have performed fully self-consistent GQCE calculations of the periodic Anderson model on Rigetti's Aspen-4 quantum device. Our results show that GQCE correctly describes the ground state phase diagram of the correlated model, which contains the Kondo insulator, correlated metal, and Mott insulator phases [17,35,36]. In contrast to Hartree-Fock theory, the critical parameters for the associated quantum phase transitions are also accurately determined using GQCE. Our work demonstrates the current capabilities of NISQ devices in the simulation of correlated materials. More importantly, it identifies GQCE as a promising framework for performing practical quantum computations of correlated materials in the near term, where NISQ devices may offer a quantum advantage.

II. RESULTS
The hybrid quantum-classical Gutzwiller embedding framework (GQCE) is based on the Gutzwiller variational embedding approach that has been developed in References [27,31]. It was shown to be equivalent to the rotationally invariant slave-boson method in the saddlepoint approximation. We review the formalism in the supplemental materials (SM), and focus in the main text on the part of the algorithm where quantum computers may offer a quantum advantage. We then present results of GQCE calculations for the infinite PAM lattice model. These fully self-consistent calculations were performed on Rigetti's Aspen-4 quantum computer.

A. GQCE approach to periodic Anderson model
To perform a first non-trivial benchmark study of GQCE for infinite systems, we consider the periodic Anderson model on the Bethe lattice in infinite dimension, as illustrated in Fig. 1(a). It is specified by a Hamiltonian composed of an itinerant c-band, a local interacting d-orbital and their onsite hybridization [36], The itinerant c-band center and the correlated d-orbital level at site are given by c and d , respectively. U is the intra-orbital Hubbard interaction parameter on the correlated d-level site, and V determines the onsite hybridization strength between c-and d-electrons.
On a Bethe lattice in infinite dimensions or with infinite nearest-neighbour connectivity, the conduction band has a density of states in the semicircular form ρ c ( ) = where D is the half band width, and hereafter set to be 1. The model hosts a diversity of paramagnetic electronic phases, including metal, band insulator, Kondo insulator and Mott insulator and the quantum phase transformations, which have been extensively studied in the literature [17,35,36]. Highly accurate numerical results have been obtained within the framework of DMFT [17][18][19][20], which becomes exact for systems in infinite dimension. This renders the PAM with Bethe lattice an ideal testbed for hybrid quantum-classical calculations for infinite systems on NISQ devices.
In this work, we choose the particle-hole symmetric point ofĤ d with Fermi level at 0, i.e., d = −U/2, and fixed V = 0.4, U = 2. In this parameter space, the system starts with a Kondo insulating (KI) phase at c-band center c = 0, transforms to metallic (M) phase with rising c , and finally enters the Mott-Hubbard insulating (MI) phase [36]. The zero temperature critical KI−M c at KI and M phase boundary is 0.07, and the critical M −M I c leading to M to MI transformation is 1.08, obtained from the DMFT calculations with the numerical renormalization group (NRG) as the impurity solver [36,37], which are numerically exact for the model in the current study.
We consider a Gutzwiller variational ansatz with the correlator acting on Hilbert space spanned by the single particle d-orbital only. This leads to a set of coupled eigenvalue problems of the Gutzwiller embedding Hamiltonian, which provides an accurate description of the local degrees of freedom, and the effective mean-field quasi-particle Hamiltonian, as shown in Fig.1 Apart from a constant, the Gutzwiller embedding Hamiltonian can be written aŝ Here D is the coupling strength between the d−orbital and the bath orbital f , at an orbital level −λ c , as the input from calculations of the effective quasi-particle Hamiltonian (see Eq. S14 and S21). The general local one-body matrix, such as the kinetic energy renormalization matrix R introduced in Eq. S6 and coupling matrix D, bears a 2 × 2 diagonal form with degenerate diagonal elements due to spin-symmetry in the paramagnetic state.
The output R and λ matrices are calculated from ground state solution of the embedding Hamiltonian, more specifically, through the one-particle density matrix, together with the ground state energy for total energy evaluation of the whole system. To solve the ground state of the embedding Hamiltonian using VQE on quantum computers, we first translate the fermionic Hamiltonian into molecular orbital representation, obtained from a spin-restricted Hartree-Fock (HF) calculation. It is subsequently transformed into a bosonic qubit representation via parity mapping [38,39]. As the ground state is restricted to total S = 0 at half-filling N e = 2, the embedding Hamiltonian can be represented in a two-qubit basis by the Z 2 symmetry aŝ Here X i , Y i and Z i are the Pauli operators acting on qubit i, and {g} are determined by embedding Hamiltonian parameters in Eq. (5), as well as the HF molecular orbitals. The asymmetric two-site embedding Hamiltonian is slightly more complex than that of the hydrogen dimer H 2 , which is a widely used example for the application of VQE in quantum chemistry [9]. To find the ground state energy and single-particle density matrix, we use VQE with an unitary coupled-cluster ansatz [6]. For a two-electron system, the UCC ansatz at double excitation level is known to be exact. And the single-excitation has no contribution to the ground state energy according to Brillouin theorem [40]. The UCC ansatz can be reduced to a simple form in two-qubit representation using parity transformation [38,39] where θ is the variational parameter bounded by [−π, π]. |01 is the spin-restricted HF ground state wave function, which is obtained by standard self-consistent calculations using PySCF package efficiently run on classical computers [41].
To get the expectation value of the embedding Hamiltonian Eq. (6) under the UCC wave function on quantum computers, we group Pauli terms that are diagonal in a common tensor-product basis. A typical quantum circuit, composed of the initial HF state preparation, UCC ansatz, and a measurement of Pauli term X 0 X 1 , is shown in Fig.1(c). In addition to Pauli terms in the embedding Hamiltonian, Pauli term Y 0 is also measured with the optimized UCC ansatz to derive the one-particle density matrix of the embedding system. The VQE code is developed based on a quantum computing library pyQuil [25,26], where we use a simultaneous perturbation stochastic approximation algorithm to optimize the noisy objective function on real devices [42]. The GQCE calculation amount to finding a fixed-point solution of the Gutzwiller nonlinear Lagrange equations (S14-S20) presented in the Supplementary Material (SM). These equations correspond to a set of coupled eigenvalue problems.
When the error vector functions F, which is defined in equations (S19-S20) of the SM can be evaluated accurately enough, the modified Powell hybrid method [43] is the method of choice to reach the solution. It employs information about the numerical Jacobian. Since the noise level of the real quantum device due to gate infidelities and decoherence is significant [44], however, F cannot be accurately calculated within GQCE. We thus employ the "exciting-mixing" method, which uses a self-tuned diagonal Jacobian approximation, instead of the Powell method. It is implemented in the SciPy library [45] and is well suited to solve the root problem of noisy nonlinear equations. The iterative procedure starts with the solution of quasi-particle Hamiltonian defined by initial guess of {R, λ} (Eq. S14), It leads to the matrices {D λ c }, which defines the embedding Hamiltonian (Eq. S18). The embedding Hamiltonian is solved on quantum computers to get the ground state single-particle density matrix. The vector error function (Eq. S19 and S20) is evaluated, based on which the matrices {R, λ} is adjusted by the "exciting-mixing" algorithm. The procedure continues with the updated quasi-particle Hamiltonian, until convergence is reached.

B. Quantum simulation results
The GQCE calculations on the PAM are carried out in two ways. First, we use a statevector simulator, which represents an ideal fault-tolerant quantum computer with an infinite number of measurements. Second, we use Rigetti's Aspen-4 quantum device, which contains 13 qubits in total. Figure 2 demonstrates the convergence of total energy, kinetic energy renormalization factor Z ≡ R † R, and maximal element of the error vector F in our GQCE calculation on Aspen-4 as a function of iteration number. The iterative nonlinear solver starts from the Hartree-Fock mean-field solution and reaches convergence after about 20 iteration steps. The remaining steps are used to estimate the error bars. The calculation results on the real quantum device closely follow that of noiseless simulations. with fluctuations that result from the device noise. The maximal absolute value of the error vector elements in Fig.2(c) levels near 0.01, which coincides with the scale of the two-qubit CZ-gate fidelity of the device. Because of the stochastic nature of quantum computing on real devices, hereafter we report results by mean values with errors estimated from the final 20 iterations. The standard deviation is estimated to be 0.05(3%) for total energy and 0.03(5%) for Z-factor.
In Fig. 2, the center of the conduction band is set to zero, c = 0. The system is then in the Kondo insulator phase, where the local correlated d-orbital, which is located at zero energy as well, hybridizes with the cband and opens a Kondo gap. The Z-factor in Fig. 2(b) shows appreciable amount of reduction from unity, manifesting local on-site Coulomb interaction effects. The kinetic energy (or hybridization) is effectively lower than the Kondo energy scale and the Kondo gap.
Let us now consider the quantum phase diagram as we tune the position of the conduction band c . Even in this restricted parameter space, where all other parameters are held fixed, the PAM goes through a series of quantum phase transition from Kondo insulator to metal and from metal to Mott insulator. We compare our GQCE findings to the numerically exact phase boundaries at zero temperature that have been determined by DMFT calculations using NRG as the impurity solver [36,37]. To extract the phase boundaries, we calculate the change of total energy E, renormalization Z-factor and total electron filling n as a function of c . Results are shown in the upper panels of Fig.3, which also includes the numerically exact phase boundaries from DMFT+NRG for comparison.
As seen in Fig. 3(a), the total energy from noiseless simulations monotonically increases with increasing c and reaches a constant as the system crosses the metal-Mott insulator transition. The GQCE calculations on Aspen-4 follows closely the exact energy curve along the phase transformation path, yet with a sizable error bar that originates from the noise of the device. The error bar drops by several orders of magnitude in the Mott insulator phase due to an efficient treatment of the (orbitalselective) Mott insulating phase within Gutzwiller theory, which exploits that Mott localized quasi-particle bands are pinned at the chemical potential with integer filling [27]. The embedding Hamiltonian in the Mott phase has a doubly degenerate ground state, which can be written as tensor product states |00 and |11 in the two-qubit parity basis. In practice, we choose one of the states to evaluate the energy and one-particle density matrix, followed by a symmetrization in the spin-sector to recover spin-symmetry. As only shallow circuits with measurements for tensor product state are necessary, the GQCE calculations generates accurate results due to the high fidelity of the single-qubit gates on the real devices.
We compare GQCE to HF calculations, where the embedding Hamiltonian solver is chosen to be at HF meanfield level (green curves in Fig. 3). Within HF, the total energy is monotonically increasing and significantly larger than the GQCE result. Crucially, it bears no signature of the metal-insulator phase transitions. The important physical phenomenon that is not captured within HF theory is the suppression of energetically unfavorable doubly occupied sites in the Hilbert space of the correlated d-orbitals.
In Fig. 3(b), we show the kinetic energy renormaliza- tion Z-factor [46], which is a key physical concept captured by Gutzwiller theory. When the conduction band center rises above the zero chemical potential, the renormalization Z-factor drops gradually and vanishes at the metal to Mott-insulator transition. Remarkably, for the parameters shown, GQCE predicts a metal-Mott insulator transition phase boundary that is in perfect agreement with the numerically exact value obtained from DMFT. The Z-factor obtained from GQCE calculations on the Aspen-4 quantum device closely follows the exact simulated data, with sizable error bar in the Kondo insulating and metallic phase. Again, the error drops significantly in the Mott phase for reasons described above. Within the HF approximation, the renormalization factor remains constant, Z HF = 1, demonstrating that the metal-Mott insulator transition is beyond the description of HF theory.
Finally, in Fig. 3(c), we show that the variation of the total electron filling is an effective way to locate the phase boundaries. In the Kondo (Mott) insulator phases, the electron filling is equal to two (one), while it is in between the two values for the correlated metal phase. The electron filling obtained from GQCE calculations on Aspen-4 agrees well with the exact statevector simulations. Its behavior can be used to locate the phase boundaries and the obtained critical parameter values are in decent agreement with the numerically exact ones. In contrast, the HF approach can only identify the transition from the Kondo insulator to the metal. As the correlation-induced renormalization of the hybridization is not captured within HF theory, the Kondo energy scale is overestimated and the Kondo insulator phase incorrectly persists up to larger values of c , (see dotted line in Fig.3(c)).

III. DISCUSSION
The Gutzwiller method adopts a Jastrow-type variational wave function, which describes the ground state properties of a correlated model beyond an effective single-particle mean-field theory. Although there is no efficient way to evaluate the associated full Green's function, the coherent part of it can be straightforwardly calculated [31]. The resulting coherent spectral density of states (DOS), which includes coherent quasi-particle excitations, can be used to distinguish the different quantum phases in the model. It is shown in the lower panels of Fig. 3(d-f), which correspond to Kondo insulator, correlated metal and Mott insulator phases. Data from GQCE calculations on Rigetti's Aspen-4 device are shown to be in excellent agreement with exact simulation results.
In the Kondo insulator phase (Fig. 3(d)), the center correlated d-orbital hybridizes with the conduction band, resulting in a finite hybridization gap. The inset shows that the hybridization gap from GQCE calculations agrees well with the exact simulation result and is significantly reduced compared with the HF mean-field value, due to the correlation-induced renormalization of the hybridization strength. As the conduction band is lifted to c = 0.8, the system is situated in a metallic phase. The hybridization gap is still present but moves to higher energy, and the chemical potential is located at the sharp quasi-particle resonance peak. The total coherent spectral weight decreases in accordance with the smaller quasi-particle weight Z as shown in Fig. 3(b). At c = 1.3, the coherent spectral weight completely vanishes as the d-orbital becomes Mott localized at halffilling. In the Mott phase, the incoherent lower and upper Hubbard bands, together with the conduction c- band, define the band gap size and distinguish between a Mott-Hubbard versus charge-transfer insulator phase. Although the GQCE calculations at this level cannot explicitly generate the Hubbard bands [47], the band gap size and characteristics can still be resolved by varying the chemical potential and monitoring the electron filling [48].
The GQCE framework is built on the open-source "CyGutz" package, which is an implementation of the Gutzwiller embedding approach in classical computer [34]. We have developed the quantum computing module of GQCE using IBM Qiskit and Rigetti's Forest SDK [25,26,49], which will also be released as opensource. The statevector simulator in IBM Qiskit and the wavefunction simulator in Forest SDK have been employed for the noiseless simulations. The GQCE calculations on the real quantum device are conveniently performed through the quantum cloud service (QCS) provided by Rigetti. It provides a quantum machine image that is co-located with the quantum infrastructure, which allows fast virtual execution of hybrid quantum-classical programs at low latency cost. As previously mentioned, the quantum processing unit (QPU) used in this study is Aspen-4. The device contains 13 qubits in total, among which we choose qubit 0 and 1 for the calculations. The associated two-qubit CZ-gate, which is one controlling factor for the noise level of the calculation results, has a fidelity of about 95.8%. Platform-level optimizations of parametric compilation and active qubit reset, which dramatically reduce the latency in QCS platform, have been utilized in our GQCE calculations. We adopt a resource-efficient readout error mitigation technique for the measurements [26]. It first measures the expectation value ξ of each Pauli term of the objective function in the associated tensor product ground state, which basically characterizes the read-out error rate. The expectation value of the same observable in the VQE ground state is then divided by ξ, correcting the bias due to measurement errors. More advanced error mitigation approaches, such as Richardson extrapolation techniques, have been proposed and experimentally realized recently [44]. It is of great interest to explore these techniques to see how feasible and efficient they are to further improve the accuracies of GQCE calculations.
To conclude, we have successfully implemented a novel hybrid quantum-classical simulation framework for correlated materials, which is based on the Gutzwiller variational embedding theory. In combination with density functional theory, this GQCE approach can describe ground state properties of correlated multi-orbital quantum materials. The GQCE method lends well to NISQ technology as it maps the infinite lattice system to an effective interacting impurity model, which is selfconsistently coupled to a non-interacting fermionic bath. To obtain a self-consistent solution of a set of coupled eigenvalue equations, we employ VQE with a standard UCC ansatz that runs on quantum computing devices. Using Rigetti's quantum cloud service, we have performed the first fully self-consistent hybrid quantumclassical calculation of an infinite correlated model. As a non-trivial benchmark calculation we have investigated the ground state quantum phase diagram of the periodic Anderson model and found good agreement with numerically exact results.
Going forward, GQCE calculations share the favorable polynomial system size scaling of VQE approaches in solving the correlated embedding Hamiltonian. Therefore, GQCE promises to be able to consider larger embedding cluster Hamiltonians, which take multi-orbital or spatial correlations into account. These are necessary to describe the impact of short-range non-local fluctuations [50] and non-local order parameters such a d-wave superconductivity [51][52][53] or vestigial orders [54]. In the near term, a robust VQE solution of a 28-qubit Hubbardtype Hamiltonian, which is equivalent to a Gutzwiller embedding Hamiltonian of a single f -orbital site in rareearth and actinide materials, would bring the capabilities of GQCE calculations on NISQ devices to the verge of what is currently possible on classical computers, thus demonstrating practical quantum advantage. approximation (GA) [S1, S5-S7], which becomes exact for systems in infinite dimension or with infinite coordination number [S4]. It is equivalent to rotationally invariant slave-boson method with saddle-point approximation [S2, S8, S9]. For nonlocal one-particle density operator or hopping operatorĉ † RiαĉR jβ , which cannot be fully represented in the reduced local Hilbert space spanned by the correlated orbitals at a single site, the expectation value has a closed form as that of the noninteracting wavefunction subject to additional renormalization with R the renormalization matrix.f operator is introduced to distinguish it from the physicalĉ operators defined in the Hamiltonian. It describes a physical local electron correlation-induced renormalization effects on the kinetic energy of the systems. For any local observableÔ i [{ĉ † iα }, {ĉ iα }], which is defined in a local correlated Hilbert space, the expectation value can be rigorously evaluated through the local reduced many-body density matrix, or equivalently the ground state wavefunction Φ i of the Gutzwiller embedding HamiltonianĤ emb i at site i with the embedding Hamiltonian of the following form The embedding Hamiltonian essentially describes a physical impurityĤ loc i coupled to a quadratic bath λ c of the same orbital dimension, with a coupling matrix D.
The total energy of the system per unit cell can now be written as withT It consists of contribution from the expectation value of an effective quasi-particle Hamiltonian and that of the quantum embedding Hamiltonians.

III. LAGRANGE EQUATIONS
The total energy in Eq. S9 is to be minimized in the parameter space defined the noninteracting wavefunction Ψ 0 and local Hilbert space of each nonequivalent embedding Hamiltonian, subject to the Gutzwiller constraints The constrained minimization can be conveniently formulated with the following Lagrange function