Few-qubit quantum-classical simulation of strongly correlated lattice fermions

We study a proof-of-principle example of the recently proposed hybrid quantum-classical simulation of strongly correlated fermion models in the thermodynamic limit. In a"two-site"dynamical mean-field theory (DMFT) approach we reduce the Hubbard model to an effective impurity model subject to self-consistency conditions. The resulting minimal two-site representation of the non-linear hybrid setup involves four qubits implementing the impurity problem, plus an ancilla qubit on which all measurements are performed. We outline a possible implementation with superconducting circuits feasible with near-future technology.


Introduction
Using highly controllable quantum devices to study other quantum systems, i.e., quantum simulation [1,2,3,4], offers a means to tackle strongly correlated fermion models that are intractable on classical computers. This is vital for understanding complex quantum materials [5] with strong electronic correlations that exhibit a plethora of exciting physical phenomena of immediate technological interest. Examples of such effects include the Mott metal-insulator transition [6,7], colossal magnetoresistance [8], and high-temperature superconductivity [9,10].
Classical numerical methods have limited ability to study even significantly simplified toy models of strongly correlated fermions. For instance, exact diagonalization faces exponential scaling with the system size, while quantum Monte Carlo methods [11,12] are often crippled by the infamous fermionic sign problem [13]. Tensor network methods [14,15,16,17,18] are powerful in one spatial dimension where they track strong correlations accurately. However, in higher dimensional systems, correlations tend to grow more quickly with system size, making these methods computationally challenging.
Another well-established approach to the study of strongly correlated fermionic lattice systems is dynamical mean-field theory (DMFT) [19]. It reduces the complexity of the original problem, e.g., the Hubbard model [20] in the thermodynamic limit, by mapping it onto a simpler impurity problem that is subject to a selfconsistency condition relating its properties to those of the original model. Since an impurity problem is local, the mapping corresponds to neglecting spatial fluctuations. In the limit of infinite spatial dimensions this mapping is exact, but for arXiv:1606.04839v1 [quant-ph] 15 Jun 2016 finite dimensions it is an approximation. Nonetheless for lattice geometries with a large coordination number, self-consistently solving the impurity problem can yield an accurate approximate solution to the original Hubbard problem.
The 'impurity' itself consists of a single lattice site taken from the original problem, and so inherits on-site interactions from the Hubbard model. This impurity site is then immersed into a time-dependent, self-consistent mean-field with which it can dynamically exchange fermions. The mean-field thus attempts to model the rest of the lattice and by being dynamical can describe retardation phenomena. Overall the impurity problem can be represented by a Hamiltonian in which the interacting impurity site is coupled to a discrete set of non-interacting 'bath' sites. The bath sites represent the mean-field and if there is an infinite number of them then the self-consistency condition is guaranteed to be satisfied. However, in practical implementations only a finite number of bath sites are used, which restricts the frequency resolution of the bath so self-consistency condition can only be fulfilled approximately. Nevertheless, many strongly-correlated features, e.g., the Mott transition, are still be captured correctly [19].
Although DMFT maps a Hubbard model to an impurity model this is still a non-trivial quantum many-body problem to solve because of the interactions at the impurity site. It is usually solved by classical numerical methods, e.g., specialised versions of those used to tackle the original problem, which attempt to keep track of the quantum correlations between impurity and bath sites. Again this limits the number of bath sites that can be treated accurately.
Here, we consider an alternative approach where the impurity problem is solved with a quantum simulator, thus avoiding many issues that are inherent to the classical methods. Quantum simulation of fermionic models has so far been mostly restricted to the analogue paradigm, especially with ultracold atoms in optical lattices [21]. Digital simulation approaches, akin to universal quantum simulators [22], have started to emerge in recent years, for example based on superconducting circuits [23,24,25,26]. The number of qubits in these digital simulators is, however, presently rather small. A direct implementation of the Hubbard model would suffer from severe finite size effects. It is nevertheless still possible for a digital quantum simulator with a restricted number of qubits to describe fermionic models directly in the thermodynamic limit when the DMFT approach is adopted.
To demonstrate this method we focus on the minimal incarnation of DMFT, the so-called "two-site" DMFT [27], where the impurity model consists of one impurity site and only one bath site, both with local Hilbert space dimension four, subjected to two specially chosen self-consistency conditions. Since two-site DMFT considers only the smallest possible impurity model, the approach cannot match the accuracy of full DMFT, but it can still give a qualitatively correct description of the infinitedimensional Hubbard model, and its simplicity makes it a good starting point before advancing to more accurate schemes. For explicit details of two-site DMFT and its features compared to full DMFT we refer to Ref. [27].
The two-site system corresponds to four qubits, two for the impurity site and two for the bath site, while a fifth, ancillary qubit is used for measurements. This number of qubits is readily available in current digital quantum simulator platforms, with IBM having made a five-qubit quantum processor available to the public [28]. A nine-qubit processor has already been demonstrated in superconducting circuits [29,25,26]. Trapped-ion technologies also allow for digital quantum simulations with up to six qubits [30,31]. Being commensurate with current state-of-the-art technology is a further justification for studying this minimal model. Our scheme is readily generalisable to a larger number of qubits allowing for more accurate simulations and potentially offering an exponential speed-up over classical Hamiltonian-based DMFT methods [32].
The self-consistency conditions are taken care of iteratively in a classical feedback loop, which thus completes the non-linear, hybrid quantum-classical device we introduce. Dynamical mean-field simulations have already been proposed for such hybrid devices [33,32]. Quantum gates similar to the ones needed in the two-site scheme have been used in demonstrating digital quantum simulation of fermionic models with superconducting circuits [34,25]. We thus focus on superconducting circuits as a candidate platform, although, e.g., trapped ions [35,30,36,37] could also be considered.
This paper is organised as follows. In Section 2, we further elucidate the framework of DMFT applied to the Hubbard model in infinite dimensions. Section 3 introduces the two-site DMFT scheme in detail. Section 4 discusses the implementation of this two-site scheme with special attention to superconducting circuits. In Section 5, we show the results of our analysis. We end with a summary in Section 6 and give an outline of the single-qubit interferometry measurement scheme in the Appendix.

Hubbard model in infinite dimensions and dynamical mean-field theory
A standard model to describe strongly correlated electron systems in thermodynamic equilibrium is the Hubbard Hamiltonian In this model, electrons with spin projections σ =↓, ↑ 'hop' between adjacent lattice sites with tunnelling energy t. This process is described in the first term, where j, k denotes the sum over all nearest-neighbour sites j and k, andĉ † j,σ andĉ k,σ denote the fermionic creation and annihilation operators, respectively. The electrons interact with on-site Coulomb repulsion U > 0, described in the latter term by the product of the local number operatorsn j,↓ =ĉ † j,↓ĉ j,↓ andn j,↑ =ĉ † j,↑ĉ j,↑ . Here, we consider the paramagnetic Hubbard model in an infinite-dimensional Bethe lattice in the thermodynamic limit at zero temperature. This setup has very simple self-consistency relations, which makes it an ideal test-bed for a proof-ofprinciple demonstration of a hybrid quantum-classical scheme.
The DMFT approach [19] to solving this model consists in neglecting spatial fluctuations around a single lattice site and replacing the rest of the many-body lattice in the thermodynamic limit by a time-translation-invariant, self-consistent mean-field ∆(τ −τ ) (or ∆(ω) in the frequency domain), as illustrated in Fig. 1a. The isolated lattice site can dynamically exchange fermions with the mean-field at time instants τ and τ . This allows one to include retardation effects that are important is the impurity Green function and G R latt,jj (ω) is the local part of the lattice Green function. (b) In Hamiltonian-based DMFT methods, one considers an impurity model which describes the local part of the Hubbard model directly and represents the mean-field as a set of non-interacting bath sites that are connected to the central, interacting impurity site. (c) The minimal representation of DMFT involves the impurity site, with on-site interaction U and chemical potential µ, coupled via the hybridization energy V to only one bath site. The bath has on-site energy c that corresponds to the mean-field ∆(τ − τ ) and is subject to two self-consistency conditions.
in the presence of strong correlations. In short, the dynamical mean-field approach reduces the complexity of the full Hubbard model to an effective single-site system which is a slightly more benign many-body problem to solve. In infinite dimensions, DMFT becomes exact as the irreducible self-energy of the lattice model becomes strictly local in space, Σ latt,jk (ω) = δ jk Σ latt,jj (ω), and its skeleton diagrams agree with those of a single-site, or impurity, model [19].
The solution of the effective single-site, or impurity, problem also yields the solution of the infinite-dimensional Hubbard model due to the self-consistency condition. This leads to the retarded single-particle impurity Green function in the frequency domain being given by where µ is the chemical potential, and Σ imp (ω) denotes the impurity self-energy. We set = 1 throughout the paper. The impurity Green function describes the response of the many-body system after a localized removal or addition of a particle on the impurity site and is defined in the time domain and at zero temperature as where i is the imaginary unit, τ is real time, {·, ·} denotes the anticommutator, θ(τ ) is the Heaviside step function, and the average is computed in the ground-state |GS of the impurity model. The fermionic creation and annihilation operators are given in the Heisenberg picture. In the paramagnetic phase the Green function is spin symmetric and we therefore only need to work out G R imp (ω) for one spin configuration.
The initially unknown mean-field ∆(ω) has to be chosen such that G R imp (ω) matches the local part of the retarded lattice Green function G R latt,jj (ω), i.e., where j is the (randomly chosen) lattice site from which the removal or addition of a particle occurs in the translationally invariant lattice model. The DMFT selfconsistency condition Eq. (4) implies i.e., the impurity self-energy matches the local self-energy of the Hubbard model in the infinite-dimensional Bethe lattice.
In the general case, the DMFT self-consistency loop is iterated as follows (see also Ref. [19]). (i) First, guess the local self-energy Σ latt,jj (ω). (ii) The local lattice Green function can be computed as is the non-interacting density of states of a Bethe lattice. The constant t * emerges from the requirement that the Hubbard hopping needs to be scaled as t ∼ t * / √ z to avoid a diverging kinetic energy per lattice site in the limit of infinite coordination, z → ∞ [19]. (iii) With Eqs. (4) and (5), we obtain ∆(ω) from Eq. (2) and the impurity model is then defined. (iv) Compute the impurity Green function and obtain the impurity self-energy Σ imp (ω). There are several means to do this [19]. (v) Set Σ new latt,jj (ω) = Σ imp (ω). (vi) Check if the self-energy has converged. If not, go to step (ii) and repeat.
Once self-consistent, the solution of the impurity problem then gives access to local single-particle properties of the original lattice model. For example, the local lattice spectral function is given by where η is a positive infinitesimal. In Hamiltonian-based impurity solvers, one parameterizes ∆(ω) by a set of bath sites (see Fig. 1b). For any finite number of bath sites, the self-consistency condition (4) can only be approximately satisfied and in the extreme "two-site" DMFT it turns out to be more suitable to reformulate Eq. (4) in a manner specially focused on this minimal representation [27] (see Section 3). Note that two-site DMFT is only able to provide a qualitatively correct description of the Hubbard model even in infinite dimensions [27].

Quantum simulator based on two-site DMFT
In terms of the single-impurity Anderson model (SIAM), the smallest impurity problem involves one fermionic site corresponding to the impurity and only one fermionic site corresponding to the entire mean-field as described in the previous section. Since two qubits are needed to encode the local Hilbert space of a fermionic site, we only require four physical qubits to implement this representation in the lab. The SIAM Hamiltonian for only one bath site readŝ Here, U is the Hubbard interaction at the impurity site 1, and µ is the chemical potential that controls the electron filling in the grand canonical ensemble. Furthermore, c and V describe the on-site energy of the non-interacting bath site 2 and hybridization between the impurity and the bath site, respectively, and give the mean-field as See Fig. 1c for illustration of the two-site SIAM. The parameters c and V are initially unknown and they need to be determined iteratively such that two selfconsistency conditions are satisfied. For details of the derivation and motivation of these conditions we refer to Ref. [27].
The first condition is that the electron filling n imp of the impurity site and the filling n = n j↓ + n j↑ of the lattice model match, i.e., The second self-consistency condition is given by where quasiparticle weight reads In Eq. (9), is the second moment of the non-interacting density of states, and the final equality follows from the semicircular density of states of the Bethe lattice.

Two-site DMFT protocol
The hybrid quantum-classical device implementing two-site DMFT consists of a fewqubit digital quantum simulator in which the impurity Green function is measured and of a classical feedback loop in which the parameters of the two-site SIAM are updated. The two-site DMFT protocol is summarized in Fig. 2 and proceeds as follows (see also Ref. [27]).
1. First fix U and µ to the desired values in the SIAM and set the unknown parameters c and V equal to an initial guess. 2. Measure the interacting Green function iG R imp (τ ). This can be done using, e.g., single-qubit interferometry (see details in the Appendix). 3. After Fourier-transforming the impurity Green function, the impurity selfenergy is obtained classically from the Dyson equation Here, the non-interacting impurity Green function is given by From the derivative of the self-energy one obtains the quasiparticle weight Z which directly yields the updated hopping parameter V via Eq. (9). The update for c is found by minimizing the difference |n imp − n| [27]. 4. Steps 2 and 3 need to be repeated until V and c are self-consistent, and n imp = n. The self-consistent Green function G R imp (ω) and self-energy Σ imp (ω) thus obtained are used to calculate approximations to local single-particle properties of the Hubbard model.

Quantum algorithm for the single-impurity Anderson model with superconducting circuits
Here, we consider the quantum gates of the digital quantum simulator part in Fig. 2, with special focus on superconducting circuits as the platform of choice [34,25,26].

Jordan-Wigner transformation of the SIAM
To implement the two-site SIAM with qubits, the fermionic creation and annihilation operators need to be mapped onto tensor products of spin operators which then act on the qubits via quantum gates. In order to obtain as simple quantum gates as possible in Section 4.3 and in the Appendix, we consider an ordering of the qubits where the first two qubits encode the spin ↓ for both fermionic sites while the last two correspond to spin ↑. This is achieved via the Jordan-Wigner transformation given explicitly aŝ andĉ jσ = ĉ † jσ † . With this mapping the hybridization terms in the SIAM de- and The number operators becomê and thus the interaction term can be written as up to a constant. The total Hamiltonian then readŝ where we have dropped constant terms.

Quantum gates in superconducting circuits
We now consider how the Jordan-Wigner transformed SIAM in Eq. (24) can be implemented in an experimental arrangement based on superconducting circuits. We present two alternative approaches. The first one couples the qubits with a transmission line resonator, which leads to the so-called XY gate between the qubits. The second approach is the Controlled-Z φ (CZ φ ) gate, which can be obtained via a capacitive coupling of nearest-neighbour transmon qubits without using a resonator. These CZ φ gates have been implemented with high fidelities of above 99% for a variant of transmon qubits called 'X-mon' qubits [38].
XY gates with resonators -The basic Hamiltonian coupling a set of qubits to the resonator has the form of a detuned Jaynes-Cummings model. By adiabatically eliminating the resonator one obtains, when the resonator is in the vacuum state, the well-known XY model for a pair of qubits l and m aŝ Here, ∆ is the detuning between the qubit level-spacing and the resonator mode, g l is the coupling constant between qubit l and the resonator, andσ x andσ y are Pauli operators. The XY gate is universal for quantum computation and simulation in combination with single qubit gates, and is the natural interaction customarily employed in superconducting circuits.
CZ-φ gates with capacitive couplings -To perform the CZ-φ gate, one qubit is kept at a fixed frequency while the other carries out an adiabatic trajectory near an appropriate resonance of the two-qubit states. By varying the amplitude of this trajectory one can tune the conditional phase φ. The unitary for the CZ φ is given by

Quantum gate decomposition of the time-evolution operator
In order to use quantum gates for time-evolution, we utilize a Trotter decomposition of the time-evolution operator corresponding toĤ SIAM in Eq. (24). The first order Trotter expansion is given bŷ Here, N is the number of Trotter (i.e., time) steps and τ N is the size of the time step. In what follows, we use the two alternative approaches for quantum gates outlined in Section 4.2 to implement Eq. (27).
XY gates -As shown in Section 4.2, the XY gate, given by the expression XY = exp −i V 2 (σ x lσ x m +σ y lσ y m ) τ n , naturally appears when considering the use of a resonator quantum bus [39]. The quantum circuit for a single Trotter step with these gates is shown in Fig. 3a.
a. Method XY and quantum bus We consider a quantum algorithm based on the XY interaction or exchange gate, which naturally appears in cQED when considering the use of a quantum bus [1]. The quantum circuit for a single Trotter step is b. Method CZ-gate Moreover, it is also possible to describe the digital simulation in terms of Trotter steps involving the optimised gates (CZ-) [2]. We first consider the Hamiltonian in terms of zz interactions. We take into account that 3 where R i ↵ (✓) = exp( i ✓ 2 ↵ ) is the rotation along the ↵-axis of the ith qubit. The evolution operator of the Hamiltonian in Eq. (5) in terms of zz interactions is The sequence of gates for one Trotter step in the digital simulation of the SIAM model with four qubits is The gates A and B are two qubit gates in terms of zz interactions, A = exp i V A single Trotter step contains 5 ZZ two-qubit gates between nearest-neighbour qubits, 2 SWAP gates, and 20 single-qubit rotations. We note that a SWAP-gate amounts to three CZ-gates, and a ZZ-gate amounts to two CZ-gates, in general. This number can optimised further if we consider different orderings for odd and even Trotter steps, in which subsequent gates may be suppressed. This reorganisation of interactions do not affect in principle the digital theoretical error. That way, an odd Trotter step reads and one even Trotter step is CZ-φ gates -To be able to utilize the CZ-φ gates, we write the time-evolution operator in Eq. (27) in terms ofσ z -σ z (ZZ) interactions, taking into account that andσ y lσ where R (l) is the rotation along the α-axis of the lth qubit. Note that in the computational basis, one can write, e.g., where we have neglected global phases. Thus, we have the decomposition where the tunable CZ φ -gate is given by Eq. (26).
and one even Trotter step is The time-evolution operator in Eq. (27) in terms of ZZ interactions is given bŷ where R The sequence of gates for one Trotter step is depicted in Fig. 3b.
A single Trotter step contains 5 ZZ two-qubit gates (corresponding to the A and B gates in Fig. 3b) between nearest-neighbour qubits, 2 SWAP gates (for the B gate which acts on qubits 1 and 3), and 20 single-qubit rotations. We note that a SWAP-gate amounts to three CZ-φ gates, and a ZZ-gate amounts to two CZ-φ gates (see Eq. (31)). This number can be optimised further if we consider different orderings for odd and even Trotter steps as in Fig. 4, such that subsequent gates may be suppressed. This reorganisation of interactions does not in principle affect the Trotter error. Hence, for a pair of Trotter steps, the number of gates is reduced, and we may only consider 10 ZZ two-qubit gates between nearest-neighbour qubits, 4 SWAP gates, and 32 single-qubit rotations.

Results
We focus on the half-filled case, i.e., µ = U 2 and c = 0, which requires the least amount of quantum gates, since the C and D gates in Section 4 vanish. Note that since the value of c is fixed in this case, it need not be updated in the self-consistency loop. We use t * , the Hubbard hopping in infinite dimensions, as our unit of energy, hence time τ is measured in units of 1/t * . Note that τ refers here to the time in the evolution operatorÛ (τ ), not to the actual time to run the experiment. We show in Fig. 5 the state fidelities F = | Ψ(τ )|Ψ T (τ ) | 2 , where |Ψ(τ ) denotes the state obtained with exact time-evolution using the full, non-Trotterized oper-atorÛ (τ ) = exp(−iτĤ SIAM ) corresponding to the two-site SIAM in Eq. (6), and |Ψ T (τ ) is the state evolved using either the XY or CZ-φ quantum gates, for various Trotter steps N up to time τ = 6/t * . Note that the number of qubits corresponding to the two-site SIAM is fixed, leaving only N as the parameter to be varied for increased accuracy. We use the initial state |Ψ(τ = 0) =ĉ † 1↓ |GS /||ĉ † 1↓ |GS ||, where |GS is the ground-state of the two-site SIAM in Eq. (6), which is a relevant state for obtaining the impurity Green function at zero temperature (see Eq. (3)). As expected, using XY gates displays superior fidelities, since CZ-φ gates require an extra factorization of the hybridization term (see Section 4). For N = 24 steps, the state fidelity using XY gates remains over 99% throughout the evolution. In what follows, we use only XY gates for the time-evolution for concreteness.
As shown in Section 3, the main object of interest is the retarded impurity Green function. One possibility to measure iG R imp (τ ) is single-qubit interferometry (see the Appendix for details), which raises the total number of qubits in the experimental arrangement to five. In Fig. 6 we plot the impurity Green function obtained from evolving the state with XY gates compared to exact evolution of the two-site SIAM for different N . We see that the Green function from the XY approach starts to follow the curve of the exact Green function better for increasing N . In our Trotter steps up to time τ = 6/t * using the XY method (blue diamonds). Comparison is given to the exact Green function (red dashed line). We set U = 4t * and V = t * . subsequent analysis, we use N = 24 up to τ = 6/t * to study what two-site DMFT physics can be captured with the digital approach.
To obtain the impurity Green function in the frequency domain, we first consider some known and general analytic properties of the retarded Green function in Eq. (3). This Green function can be written as a sum of the particle and hole contributions as where |j is an eigenstate ofĤ SIAM with eigenenergy E j , and ω j = E j −E GS . In twosite DMFT, the interacting Green function is a four-pole function [27], which limits the number of terms in the above summation to four. Moreover, in the presence of particle-hole symmetry, we have j|ĉ † 1σ |GS 2 = j|ĉ 1σ |GS 2 , and Eq. (33) can be written as where α j = j|ĉ † 1σ |GS 2 . Thus, to obtain the impurity Green function in the frequency domain as we need to extract the unknown residues α j and poles ω j by fitting an expression of the form in Eq. (34) to the measurement data of iG R imp (τ ), as shown in Fig. 7a. This method to determine α j and ω j is far more reliable and requires fewer time steps than numerically Fourier-transforming the iG R imp (τ ) data. It can also be readily generalised to larger systems by including more terms in the sum in Eq. (33). Figure 7b shows the real part of the impurity Green function in the frequency domain, Re G R imp (ω + iη) , with residues and poles obtained from the fit in Fig. 7a, while in Fig. 7c we plot the real part of the impurity self-energy, Re [Σ imp (ω + iη)], obtained utilizing the Dyson equation (11). We clearly see the four-pole structure of the Green function, while the self-energy has two poles. The results are in excellent agreement with the exact solution of the two-site SIAM, with the poles of the selfenergy using fitted α j and ω j differing from the exact solution by 2%.  Figure 9 Quasiparticle weight as a function of interaction U . Self-consistent quasiparticle weight Z obtained from the XY method with 24 (blue diamonds), 36 (red circles), and 48 Trotter steps (yellow squares), compared to the exact solution of the two-site SIAM (purple stars). Inset: Same plot zoomed into the region around the critical interaction, Uc = 6t * .
Once we have obtained the impurity Green function, and thus the impurity selfenergy, we proceed according to the two-site DMFT protocol in Section 3 until self-consistency has been reached. In DMFT we are interested in the local lattice spectral function A latt,jj (ω) which, at self-consistency, is given by the impurity spectral function A imp (ω). In the paramagnetic phase of the infinite-dimensional Hubbard model, the spectral function has a three peak structure with an upper and a lower Hubbard band, corresponding to empty and doubly occupied sites, respectively, and a quasiparticle peak with integrated spectral weight Z between the bands [19]. In two-site DMFT, since the self-energy has two poles, this three peak structure can be qualitatively reproduced with the spectral function [27] where ρ 0 is the non-interacting density of states of the Bethe lattice. Figure 8 shows the spectral function in Eq. (36) where the impurity self-energy has been obtained both from the XY method and from exact numerics of the two-site SIAM using the interactions U = 5t * and U = 8t * . We notice that for U = 5t * , the Hubbard bands from the XY method are slightly dislocated and the quasiparticle peak is slightly narrower compared with the exact solution of the two-site SIAM, but the agreement is still very good. The overall shape of the spectral function from the XY method is unchanged compared to the exact case. This underestimation of the width of the quasiparticle peak stems from the fact that the fitting procedure in Fig. 7a causes the negative of the derivative of the self-energy in the XY method to be a bit larger than the exact value from the two-site SIAM, i.e., − , which leads to Z in Eq. (10) from the XY method to be slightly smaller than in the exact solution of the two-site SIAM, i.e., Z XY Z exact . For U = 8t * , the two spectral functions agree with maximum relative error of 10 −8 , since in this case V = 0 is found to be the self-consistent solution, whence the Trotterized evolution operator in Eq. (27) matches full evolution operator of the two-site SIAM, and thus there is no Trotter error. We observe that in Fig. 8 the central quasiparticle peak vanishes, which is characteristic of insulating behaviour. See Ref. [27] for a discussion of the artifacts of the spectral functions in two-site DMFT compared to full DMFT.
To study the transition between the two types of spectral functions in Fig. 8, we plot in Fig. 9 the self-consistent quasiparticle weight Z obtained from the XY method as a function of the interaction U for different Trotter steps N . We also show Z from the exact solution of the two-site SIAM for comparison. We see that the digital approach captures the correct trend of the curve, but in the metallic side underestimates to a small degree the values of Z for interactions close to U = U c = 6t * , which is the critical interaction for Mott transition in two-site DMFT at halffilling [27]. These results are consistent with the spectral functions in Fig. 8. The underestimation of Z can be diminished by increasing N , as shown in Fig. 9. It is noteworthy to mention that two-site DMFT overestimates the quasiparticle weight compared to full DMFT for interactions U < U c , as demonstrated in Ref. [27]. Above U c , we find Z = 0 to be the self-consistent solution, corresponding to the insulating phase.

Summary
We have proposed a quantum algorithm for two-site DMFT to be run on a small digital quantum simulator with a classical feedback loop, allowing the qualitative description of the infinite-dimensional Hubbard model in the thermodynamic limit. We have considered two alternative quantum gate decompositions consistent with state-of-the-art technology in superconducting circuits for the time-evolution operator. We found that an increasing number of Trotter steps improves the fidelity of our digital scheme to qualitatively describe the Mott transition. Our work therefore provides an interesting application for small-scale quantum devices. It also paves the way for more accurate quantum simulations of strongly correlated fermions in various lattice geometries, which are relevant to novel quantum materials, when the general self-consistency condition and larger number of qubits are used.