Hybrid quantumclassical approach to correlated materials
Abstract
Recent improvements in control of quantum systems make it seem feasible to finally build a quantum computer within a decade. While it has been shown that such a quantum computer can in principle solve certain small electronic structure problems and idealized model Hamiltonians, the highly relevant problem of directly solving a complex correlated material appears to require a prohibitive amount of resources. Here, we show that by using a hybrid quantumclassical algorithm that incorporates the power of a small quantum computer into a framework of classical embedding algorithms, the electronic structure of complex correlated materials can be efficiently tackled using a quantum computer. In our approach, the quantum computer solves a small effective quantum impurity problem that is selfconsistently determined via a feedback loop between the quantum and classical computation. Use of a quantum computer enables much larger and more accurate simulations than with any known classical algorithm, and will allow many open questions in quantum materials to be resolved once a small quantum computer with around one hundred logical qubits becomes available.
I Introduction
The current workhorse for materials simulation is densityfunctional theory (DFT) [1]. DFT circumvents the exponential scaling of resources required to directly solve the electronic quantum manybody Hamiltonian by mapping the problem of finding the total energy and particle density of a system to that of finding the energy and particle density of noninteracting electrons in a potential that is a functional only of the electron density, and requiring selfconsistency between the density and potential. The noninteracting electron and selfconsistency problems are manageable on classical computers. While the universal functional that gives the dependence of the potential on the electron density cannot be efficiently evaluated [2], several approximations to it are known empirically to be good; for example the local density approximation (LDA) [3]. LDA and related approximations yield reliable results for many weakly correlated materials, such as band insulators, metals, semiconductors and classes of biomolecules. However, many effects, including the intriguing physics of strongly correlated transition metal materials near a Mott transition [4], the phenomenon of hightemperature superconductivity [5], the properties of heme and other molecular complexes involving transition metals [6], and the complex physics of the actindes [7] are beyond the scope of DFT.
To go beyond DFT, one can attempt to exactly solve the quantum manybody problem using methods such as full configuration interaction (FCI), however because the required computational effort on classical computers scales exponentially in the number of orbitals, the method is limited to relatively small systems. In contrast, it has been shown that quantum computers can in principle solve certain electronic structure problems [8; 9] in polynomial time. Recent advances towards building a small quantum computer [10; 11; 12] have led to increasing interest in what a small quantum computer could realistically simulate, and it has been shown that the simulation of small molecules [13; 14; 15; 16; 17] and simplified model Hamiltonians [18] is within reach. However, in part because of the multiplicity of relevant interaction terms, the scaling of the currently known algorithms is not benign enough to allow naive direct simulations of complex correlated materials, for which thousands of electrons would have to be considered.
To make progress, the problem must therefore be simplified. One approach is to approximate the material with idealized model Hamiltonians, such as the Hubbard model [19], which have a sufficiently small number of interaction terms that they can easily be studied on quantum computers [18]. While capturing qualitative phenomena, such simple models do not offer a quantitative description of real materials.
Ii A hybrid quantumclassical approach
We thus propose a hybrid approach, combining classical and quantum algorithms within the framework of the DFT plus dynamical mean field theory (DMFT) embedding approach. In this framework a computationally inexpensive DFT calculation is used to define a set of orbitals and determine the electronic structure of the majority of the orbitals, while a more expensive manybody method (here, DMFT) is used to solve a reduced model involving a much smaller set of correlated orbitals. The simplification lies in that one must deal only with a small set of correlated orbitals; the tradeoffs are that one requires a relatively complete solution of the correlated problem (full frequency dependence of the Green’s function) and one obtains only an approximation to the true answer. The required highly accurate solution of a small problem is an ideal application for a small quantum computer where it enjoys maximal benefits over all known classical algorithms.
This approach is being successfully employed in calculations of properties of correlated materials [20]. DMFT provides an approximation to the solution of a full correlated problem by leveraging the solution of a quantum impurity problem, in which a finite cluster of interacting orbitals is selfconsistently coupled to a bath of noninteracting electrons. DMFT becomes exact in an infinite lattice coordination limit or when the momentum dependence of the self energy may be neglected [21], or when the number of orbitals in the impurity model becomes equal to the number of orbitals in the original problem to be solved [22]. Practical implementations require impurity models with a relatively small number of orbitals, thus providing only an approximate solution to the full model; however in many cases the approximation is reasonably good. DMFT has been very successful at qualitatively describing the Mott transition [23; 24] and its ‘cluster’ extensions [22] have produced important results for model systems (in particular the two dimensional Hubbard model). The combination with DFT is quantitatively explaining some properties of correlated materials [25; 20].
However, within a classical computational framework the complexity of the impurity model scales exponentially with the number of orbitals, placing severe limitations on the types of materials that can be tackled and restricting most realmaterials DFT+DMFT simulations to the “singlesite” DMFT approximation involving an impurity model representing a single correlated atom and neglecting all momentum dependence of the electron self energy. The restriction to just a small set of correlated orbitals also means that the method cannot be used to test the embedding hypothesis by systematically increasing the number of kinds of orbitals treated as “correlated”. We show here that these limitations can be overcome by using a quantum computer to solve the impurity problem. Even an impurity problem with only degrees of freedom would enable the study of fundamentally new problems by allowing materials with multiple correlated atoms per unit cell to be considered, allowing cluster DMFT calculations of real materials and (given some further development of DMFT methodology) examination of the embedding hypothesis and the closely associated “double counting correction”.
The Hamiltonian of the embedded impurity problem in DMFT may be written
(1) 
where creates a fermion in one of the spin orbitals of the interacting system labeled by a combined spin and orbital index ; creates a fermion on one of the bath sites. is finite by construction and , but many approaches approximate the problem using a finite number of bath sites.
While the hopping integrals and interaction integrals are directly given by the underlying material, the bath coupling and bath energies are determined from a selfconsistency condition involving the Green’s function of the impurity model and appropriate matrix elements of the Green’s function of the lattice model. The selfconsistency condition is typically solved iteratively by repeating the following steps: (i) Starting with an initial guess for the bath parameters and , solve for the ground state of Eqn. (1) and extract the impurity Green’s function . (ii) From , calculate the selfenergy and an updated noninteracting Green’s function . (iii) Determine the discrete bath parameters and to closely match the desired noninteracting Green’s function. In practical calculations about twenty solutions of the impurity model are required. For further details, we refer to Appendix C.
The impurity solver is by far the most computationally demanding step in this loop. A variety of approaches exist [24; 26]. Viewed as a generic quantum problem, the model has an interesting sparsity structure: the interactions only couple the impurity states, while the bath states are noninteracting. This sparsity structure is exploited by the widelyused continuoustime quantum Monte Carlo (QMC) [26] methods, in which the bath is integrated out leaving an action involving degrees of freedom. However, these QMC algorithms [27; 28] scale exponentially in and currently remain limited to . Also they suffer from a severe sign problem in low symmetry situations where there is no choice of basis that diagonalizes the hybridization function at all frequencies.
Exact diagonalization (ED) solvers [29] approximate the continuous bath by a finite number of bath orbitals and do not take advantage of the sparsity structure. Therefore on a classical computer they have a cost that is exponential in the total system size . This and the need to obtain the full Green’s function means that practical calculations are limited to , in other words to five or fewer correlated orbitals, often corresponding to just a single correlated atom within a unit cell, and with a very small number of bath sites per correlated orbital. Recent developments [30; 31; 32] based on quantum chemical methods to define reduced basis sets for the ED calculation permit inclusion of somewhat larger numbers of bath orbitals, but at least as presently formulated these methods work in a natural orbital basis which strongly mixes the bath and correlated orbitals, so that the sparsity structure mentioned above cannot be exploited. In a parallel development, ideas to solve the impurity problem using tensor networks [33] have recently started to show great promise [34; 35; 36].
Iii Quantum algorithm for the impurity solver
Significantly larger problems can be tackled when solving the impurity problem on a quantum computer. The key points are that the wavefunction of the impurity problem requires only logical qubits and that the quantum computation takes advantage of the sparsity structure mentioned above. In particular, the number of bath sites affects the number of required qubits, but does not have a strong effect on the computation time. We leverage standard quantum algorithms discussed previously for quantum chemistry and model Hamiltonian applications to first obtain a quantum representation of the ground state of (1), and then measure the Green’s function.
To obtain the ground state, we combine adiabatic state preparation and quantum phase estimation (QPE) [37; 38]. We start from the easily prepared ground state of a simple Hamiltonian and evolve it under a timevarying Hamiltonian that adiabatically interpolates from to the desired Hamiltonian (1). Changing the parameters slowly compared to the inverse spectral gap of ensures that the wave function always remains close to the ground state. Possible choices for the initial Hamiltonian could be either the atomic limit of turning off all hopping terms, such that the ground state becomes a simple product state of occupied and unoccupied spin orbitals, or turning off interactions such that the initial state is a Slater determinant which can be efficiently prepared using techniques discussed in Ref. 18. At the end of the adiabatic process, QPE can be used to measure the energy of the state and collapse the wavefunction into an eigenstate of the Hamiltonian. QPE relies on the ability to apply to the state, and avoids measuring the individual terms of the Hamiltonian separately, since such measurements do not commute among themselves and with and would thus destroy the state. In contrast, quantum phase estimation performs a measurement that is diagonal in the energy basis and which will project a state close to the ground state onto the ground state with high probability. This is achieved with an accuracy , where is the total computation time. Details of the method are discussed in Appendix B. State preparation has succeeded if we measure the ground state energy ( then is the corresponding ground state wave function), but needs to be repeated if the measured energy corresponds to an excited state.
Due to the constraint of unitary evolution on a quantum computer, we can only measure the Green’s function in real time (or real frequencies [18]). For , the particle and hole Green’s functions in real time are defined as:
(2) 
In the following, we present the details involved in extracting the Green’s function in Matsubara frequencies from our quantum impurity solver. We will suppress the orbital indices for notational simplicity and implicitly assume that the Green’s function is treated as a matrix. As the self consistency condition is best enforced in imaginary frequencies, we perform a Hilbert transformation to find
(3) 
where , and is a fictitious inverse temperature. To obtain groundstate properties, is chosen sufficiently small to guarantee converged results and the frequencies are cut off at some suitably chosen . In practice, the integral in (3) is replaced by a discrete sum over a logarithmic grid, and the realtime Green’s function must be measured separately for every time point. This is discussed in more detail in Appendix C. Alternative approaches to measure dynamical correlation functions are discussed in Ref. 18.
To measure the realtime Green’s functions (2), we relate them to expectation values of unitary operators , , which can be measured directly; for details, see Appendix B.3. Note that this formulation allows for straightforward determination of the superconducting components of the Green’s function, at the cost of twice as many measurements. In a naive approach, measurements are destructive and must be reprepared for each measurement.
We perform numerical simulations of our proposed quantum algorithms to establish a baseline of how many gates need to be executed to solve a simple impurity problem. Since these simulations on a classical computer scale exponentially in the size of the impurity system, we are limited to very small problems. We thus consider a single spinful impurity site () coupled to a bath of 5 spinful sites (). We run a selfconsistent DMFT calculation, based on a simulation of our quantum algorithm, for a Hubbard model on the Bethe lattice for two different strength of the onsite interaction corresponding to the Fermi liquid and the Mott insulating regime, respectively. Our results for the spectral function of the converged DMFT solution are shown in Fig. 3.
To evaluate the integral (3) in each iteration of the DMFT loop, we measure the realtime Green’s function on a grid of 1000 time points and take 400 measurements at each time point, where we need to measure both the particle and hole contributions, the imaginary and real part as well as spin components separately (unless the system is spindegenerate). These numbers have been chosen such that the error from the measurement of the Green’s function is small compared to the uncertainty in the bath fitting procedure, i.e. the limitations of DMFT with a small, discrete bath. With these choices, we need a total of measurements, each giving one bit of information. For each measurement, we prepare the ground state, which we found to require total gate operations in this instance.
For more complex problems, the preparation of the ground state will be much more costly, and should thus be avoided if possible. Since each measurement only extracts one bit of information, the state after the projective measurement may have significant overlap with the ground state. A relatively short quantum circuit, using a QPE, can then be used to project back into the ground state. This motivates the approach sketched in Fig. 2. For the simple test case considered here, the cost of performing a QPE is comparable to the adiabatic state preparation, and thus no advantage can be gained by attempting to project back into the ground state after measuring one bit. In general, however, adiabatic state preparation will scale with the square of the inverse gap, , while the cost of QPE scales only roughly linearly, , leading to quadratic advantage of QPE over repreparing the state. For large simulations, avoiding the preparation at the cost of performing QPE more often is therefore highly advantageous.
Further improvements can be gained with a fully coherent measurement procedure as described in Ref. 18 and also reviewed in Appendix D. This procedure quadratically speeds up sampling, reducing the scaling of the time required for a given accuracy from to . In our numerical example given above, the coherent approach will reduce the number of measurements needed to achieve the same accuracy from 400 to and thus yield roughly a tenfold improvement. In other applications, where a higher accuracy is required, the improvement will be more significant.
Having established a baseline for the number of gates that must be executed, we now address the question of how this number scales with the size of the impurity problem. The important contributions to the scaling are (i) the number of terms in the Hamiltonian, which determines the number of gates required to perform a single time step of the evolution, (ii) the number of measurements that must be taken, (iii) the time that is required to accurately prepare the ground state, and (iv) the time step required to reach the desired accuracy.
The number of terms in the Hamiltonian scales like (note that enters linearly while enters quartically; this is an example of exploiting the sparsity structure mentioned above). If gates can be executed in parallel, an even more favorable scaling is obtained by mapping the bath onto a set of chains (rather than the “star” topology used in (1)); this can be achieved by blocktridiagonalizing the quadratic bath terms using a Krylov approach [39]. Using this mapping, many terms can be executed in parallel, such that the scaling becomes independent of the size of the bath and scales as . The number of noncommuting terms in the Hamiltonian also modestly affects the required time step [16]. The selfconsistency condition requires measurements of the Green’s function matrix. However, in many cases orbital and spin symmetries of the system can be used to blockdiagonalize the Green’s function and thus reduce the number of independent measurements.
Since our algorithm spends significant time preparing the ground state, the requirements of the adiabatic state preparation play a crucial role in scaling. This generally depends on the minimal spectral gap along the chosen adiabatic path, which in turn will depend on the physical properties of the system in question. For gapped systems such as Mott insulators or superconductors where the order parameter fluctuations are gapped, the required preparation time is not expected to scale significantly with the size of the bath or the complexity of the impurity. Systems with very small gaps, such as Kondo systems, or in gapless regimes, such as Fermi liquids, likely pose greater challenges in particular for the accurate preparation of the ground state, where special care must be taken to find an optimal adiabatic path and choose a sufficiently long preparation time. However, it is important to note even for physical systems which are gapless, the finite size of the discrete impurity bath induces a nonzero gap in the problem we have to solve.
Taking the above scaling considerations into account, a relevant physical problem of 10 orbitals () with the corresponding number of 60100 bath sites seems within reach for a small quantum computer of about one hundred qubits. Such a problem would require on the order of measurements, which each can be achieved in a coherent run of about gates. While this leads to a large total number of gates of , it is important to keep in mind that in contrast to other approaches [13; 14], these gates need not be executed in a single coherent simulation, but are broken up into independent runs that can be executed sequentially or in parallel if several quantum computers are available. If the concrete quantum computing architecture allows the execution of more gates coherently, large improvements on the total gate count can be gained from exploiting the improved measurement schemes mentioned above.
Iv Outlook and Discussion
Our hybrid quantumclassical approach to materials simulations is similar in spirit to complete active space methods in quantum chemistry, which pick a subset of orbitals to be treated by a method with higher accuracy, but go beyond these methods by feeding back the solution of the quantum impurity problem into the DFT problem. The ideas put forward here are not restricted to the commonly used DFT+DMFT approach, but can be generalized to other quantum embedding approaches such as the recently proposed densitymatrix embedding theory (DMET) [40]. There, one strives to attain selfconsistency between an extended noninteracting lattice model and an interacting impurity problem. The parameters that must be determined selfconsistently are hopping parameters of the noninteracting model, which are chosen in such a way that the singleparticle equaltime Green’s functions of the two models match. This scheme requires knowledge only of equal time Green’s functions which are straightforward to measure on a quantum computer.
While it has been known for a long time [41; 8] that many quantum problems can be simulated on quantum computers with polynomial scaling, such scaling in the asymptotic regime is insufficient to make an algorithm practical, especially if the power of the polynomial is high and constants are large – as is the case for quantum chemistry solutions of molecules and materials [13; 14]. Recent improvements in algorithms and runtime estimates [15; 16; 17] make the solution of small but classically challenging molecules practical on small quantum computers.
With our hybrid quantumclassical algorithm small quantum computers will also be useful for the simulation of larger systems, and especially strongly correlated crystalline materials or complex molecules which exhibit a wide variety of interesting physical phenomena and have a broad range of applications. Materials simulations are today one of the major uses of supercomputing facilities, and will profit enormously from the availability of quantum computers as specialpurpose accelerators. Quantum algorithms like the one presented here have the potential to solve many of the problems that plague today’s simulation of correlated materials on classical computers, revolutionizing the field by opening new horizons for computational investigation of quantum materials.
A different class of algorithms that have been explored over the last few years are ”variational quantum eigensolvers” [42; 43]. These approaches, which are tailored for firstgeneration quantum computers that can only execute short circuits coherently, rely on Hamiltonians with only very few noncommuting terms to be practically feasible. As such, they are wellsuited to the approximate simulation of simple, local model Hamiltonians such as the singleband Hubbard model. However, simulating complex materials as discussed here will be prohibitive both in the number of required qubits to represent an extended system and in the number of relevant noncommuting interaction terms, which severely affect the scaling of the variational algorithms.
During the completion of this work, a related approach where a quantum computer is used as impurity solver within the variational cluster approach [44; 45] appeared in Ref. 46. Our approach differs in several crucial ways. Most importantly, the embedding method used in our paper is the more broadly applicable and widely used DMFT method. Furthermore, we have described a zerotemperature approach in this manuscript, while the method of [46] operates at finite temperature. Since the circuitbased quantum computers for which both approaches are designed perform unitary evolution of a quantum system, computing the required timedependent correlation functions for thermal (mixed) states incurs significant overhead both in the number of required qubits as well as the computation time. No attempt at estimating the scaling of the algorithm, or establishing a baseline for the gate counts, is made in Ref. 46.
Acknowledgements
The authors acknowledge useful discussions with Lei Wang and Ara Go. Part of this research was performed at the Aspen Center for Physics supported by NSF grant #PHY–1066293. MT was supported by Microsoft Research, the European Research Council through ERC Advanced Grant SIMCOFE and the Swiss National Science Foundation through NCCR QSIT.
Appendix A Introduction to Quantum circuits
(a)
(b)  (c)  (d) 

In this section, we show the quantum circuits necessary to implement the operations discussed in the main text. In these circuit diagrams, the horizontal lines indicate individual qubits and boxes indicate quantum gates. For example, in the following circuit, the Hadamard gate , which transforms between the and basis, is applied to one qubit:
(4) 
Other important ingredients are controlled gates, as shown in this circuit:
(5) 
Here, we first apply a controlled unitary, and then a controlledNOT (CNOT) gate. The operation ( or NOT) is applied to those components of the input state where the control qubit is in state , but not to those where the control is in state . Other gates we use are
(6) 
Fig. 4 shows how to perform the measurement of a unitary by applying the unitary controlled on an ancilla qubit that is in the state . This will entangle the ancilla qubit with the qubits on which the unitary is applied and make the expectation value accessible by measuring the ancilla.
Fig. 5 shows how a single time step of the Trotter evolution is implemented as a quantum circuit. In the present Hamiltonian, terms fall into three categories: (i) chemical potential terms of the form , (ii) hopping terms of the form , and (iii) an interaction term . Here, the subscript is a multiindex that contains both spin and orbital index. Fig. 5 shows the way these different terms are implemented as quantum circuits. Here, for sake of completeness, we always include the control qubit that is necessary to perform subsequent steps of the algorithm. In the case of the hopping circuit , in general a JordanWigner transformation has to be performed to correctly account for the fermionic sign structure. For an overview of how to achieve this, we refer to Ref. 15.
Appendix B Quantum Algorithms
b.1 Quantum simulation of time evolution
The most basic building block of our simulation algorithm is the ability to simulate time dynamics of quantum systems on a quantum computer. To this end, we first map the Hilbert space of the quantum system onto that of the quantum computer. The simplest method is to allocate one qubit per spinorbital and work in a secondquantized occupationnumber basis, i.e. use the state of the qubit to indicate whether the spinorbital is occupied or empty. Next, we need to apply the unitary to this state. While many approaches to approximating this unitary on a quantum computer are known [47; 48; 49; 50; 51; 52; 53; 54], we use here the simple approach of a TrotterSuzuki decomposition [47; 48]. We decompose the Hamiltonian as a sum of noncommuting terms , where the include both onebody and twobody terms, and make the approximation
(7) 
where is some integer. This approximation becomes exact as , and it may often be advantageous to use higherorder decompositions [55]. The simple quantum circuits that implement for the various terms are discussed in Appendix A. For the evolution under a timedependent Hamiltonian, as required for the adiabatic state preparation, we update the parameters of the Hamiltonian in each of the time steps.
b.2 Quantum phase estimation
Given the ability to apply to the state, we can measure the energy of a given state using an approach known as quantum phase estimation [37; 38]. The basic idea of quantum phase estimation is to implement an interference experiment: an additional ancillary qubit (labelled PE) is used to control the application of so that if the PEqubit is in the state , then is applied, while if it is , the identity operation is applied to the state. Effectively, this interferes two distinct trajectories with a phase difference where is the energy of the state. This phase difference rotates the angle of the qubit and by measuring this angle one can determine the phase difference. By taking large , one can make the phase difference sensitive to small energy differences, allowing precise measurement of the energy by combining the information from measurements at several different times.
b.3 Measurement of Green’s functions
Within the computational model assumed here, we can perform measurements on individual qubits in a fixed basis. In order to measure the realtime Green’s function (2), we first relate its expectation values to those of unitary operators [18]. We can then use a standard approach that allows the measurement of expectation values of unitary operators. This approach, shown in detail in Appendix A, uses a controlled unitary operation to entangle an ancilla qubit with the qubits on which the unitary should be measured, and thereby makes the expectation value accessible through a measurement of the ancilla in the computational basis.
We define the unitary operators
(8) 
where
(9) 
When applied to an eigenstate with energy , the unitary can be simplified using . In general, (8) contains terms of the form , , etc. However, assuming the absence of superconductivity, operators that do not conserve particle number will have vanishing expectation values in the ground state and we obtain
(10) 
We can thus reconstruct the desired expectation values as
(11) 
To get the real and imaginary part of both and for a given requires 4 measurements. In the presence of superconductivity, where fermion number conservation is broken down to fermion parity conservation, crossterms such as and do not vanish when evaluated on the ground state. Finding and thus requires measurement of the real and imaginary part of all of the four operators , thus increasing the total number of measurements at each time point from 4 to 8.
Appendix C Selfconsistency and bath fitting
The selfconsistency loop which is used to determine the free parameters of the bath, and , is most conveniently and reliably executed in imaginary (Matsubara) frequencies, and we must therefore extract the Green’s function in imaginary time from our quantum impurity solver. In the following, we will omit orbital indices for notational simplicity; in the general case, the Green’s function must be assumed to be a matrix. We define the particle and hole contribution to the realtime Green’s functions (cf Eqn. (2)) as
(12) 
The standard timeordered Green’s function , where is the timeordering operator, can be recovered as
(13) 
Performing a Fourier transform on the Green’s function, , we find using the above definitions:
(14) 
where the lower bound in the time integrals has been introduced for later convenience, and can be taken to be on the order of the floating point precision when numerically performing the integrals. is a numerical broadening factor that should be taken small compared to the relevant physical energy scales of the system for extracting the spectral function, but can be taken to be if only the imaginaryfrequency Green’s function is desired. Viewed as a function of complex frequencies , the manybody Green’s function is analytic in the lower and upper complex halfplane, with nonanalyticities only along the real axis , and asymptotically behaves as for . Following standard definitions, the spectral function is given by
(15) 
In this definition of the spectral function, positive frequencies encode the particle contribution, and negative frequencies encode the hole contributions. For equilibrium systems, .
Given the Green’s function on the real frequency axis, or the spectral function, we can rely on the analyticity properties mentioned above and use a Hilbert transformation of the spectral function to obtain the Green’s function in imaginary frequencies:
(16) 
Using the integrals (for , )
(17) 
one can for example compute the particle contribution to be
Here, denotes complex conjugation. The computation for the hole contribution proceeds analogously, to yield the final expression (for )
(18) 
We turn this into a discrete sum over a set of times at which the integrand is evaluated. For best convergence of the integral, we use an improved integrator scheme such as Simpson’s rule. We choose the times as
(19) 
We now describe the DMFT selfconsistency condition, which relates the free parameters of the impurity model to the parameters of the original Hubbard model. Here, we follow closely the notation of Section VI.A.1.d of Ref. 24. The selfconsistency condition can be stated as
(20) 
where is the chemical potential, denotes the selfenergy,
(21) 
is the noninteracting impurity Green’s function corresponding to a set of bath parameters, and is the impurity Green’s function measured from the solution of the interacting impurity problem. In general, the properties of the bath are encapsulated in the hybridization function . For the specific case of a discrete bath used in this paper, the hybridization function is related to the bath parameters by , and the noninteracting Green’s function is then given by [29]
(22) 
In most practical cases, the selfconsistency must be achieved iteratively by means of executing the selfconsistency loop. Given as input the noninteracting Green’s function (which in the case of a discrete bath as discussed here is parametrized by the bath parameters , through (22)) as well as the interacting impurity Green’s function obtained from the solution of the impurity problem for that set of bath parameters, one first computes the selfenergy using (21) and then evaluates
(23) 
to obtain a new estimate for the noninteracting Green’s function.
Having obtained the new noninteracting Green’s function in imaginary frequencies, we must obtain the bath parameters , (for a single spindegenerate impurity) such that the discrete version of the noninteracting Green’s function of Eqn. (22) best approximates the desired Green’s function. To this end, we optimize the cost function
(24) 
using a nonlinear optimization scheme in the parameters , . For more details on how to reliably perform optimization of the bath parameters, see Ref. 32.
In the case of the Bethe lattice [56] in the limit of infinite coordination number, the density of states follows the semicircular
form , , and 0 otherwise (here we set the hopping integral
to )
(25) 
Using (21), this can be simplified to yield the concise selfconsistency condition for the Bethe lattice [29]:
(26) 
In this case, a new noninteracting Green’s function can be obtained by simply evaluating the righthand side of (26) using the numerically obtained impurity Green’s function . This is then used as input for the same bath fitting procedure of (24).
Appendix D Gatecount estimates
d.1 Incoherent approach
The most naive possible workflow to measure an expectation value for a fixed time is as follows:

Prepare the ground state by adiabatic evolution from an easily prepared initial state, for example free fermions or the atomic limit.

Measure the expectation value of the unitary operator . This is achieved using the circuit illustrated in Fig. 4, which projects into a final state .
Since the most costly step of the above procedure is the preparation of , it would be very helpful to reduce the number of preparations. This seems feasible since the measurement only measures a single bit of information, and the final state will likely still have significant overlap with the ground state. To exploit this, we can apply a projective measurement that, if successful, projects the state back into the ground state. This can be achieved by performing QPE on the unitary to measure the energy. In the successful case, where the final measurement in the QPE yields the known ground state energy, the state has been successfully projected into the ground state and can be used as input to a new measurement; however, if the measurement yields an eigenstate of different energy, which is orthogonal to the ground state, the state would have to be reprepared. To avoid this, the measurement of the ancilla qubits in the final step of QPE can be replaced by the measurement of a single qubit which encodes whether the system is in the ground state or not. This yields the algorithm sketched in Fig. 2 of the main manuscript. The probability of returning to the ground state, which dictates how often the state must be reprepared from scratch, depends on the numerical details of the model and must therefore be estimated on a casebycase basis.
Numerical simulations of the example discussed in the main text show that following the above scheme, we need to prepare the ground state from scratch times, and have to perform QPE for a total of times.
d.2 Coherent approach
The coherent measurement procedure was first described in Ref. 18. For a detailed description of this approach, we refer to this reference. For the purpose of this paper, it suffices to note that this approach gains a quadratic speedup in the accuracy, i.e. the number of measurements required at each time step are reduced quadratically.
Footnotes
 See, e.g., Ref. 24
References
 P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864–B871 (1964).
 Norbert Schuch and Frank Verstraete, “Computational complexity of interacting electrons and fundamental limitations of density functional theory,” Nature Physics 5, 732–735 (2009).
 W. Kohn and L. J. Sham, “Selfconsistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
 Masatoshi Imada, Atsushi Fujimori, and Yoshinori Tokura, “Metalinsulator transitions,” Rev. Mod. Phys. 70, 1039–1263 (1998).
 J G Bednorz and K A Müller, “Possible high superconductivity in the BaLaCuO system,” Zeitschrift Fur Physik BCondensed Matter 64, 189–193 (1986).
 Elena G Kovaleva and John D Lipscomb, “Versatility of biological nonheme fe(ii) centers in oxygen activation reactions,” Nature Chemical Biology 4, 186 – 193 (2008).
 R.C. Albers, “Condensedmatter physics: An expanding view of plutonium,” Nature 410, 759–761 (2001).
 Seth Lloyd, “Universal quantum simulators,” Science 273, 1073–1078 (1996).
 Alán AspuruGuzik, Anthony D. Dutoi, Peter J. Love, and Martin HeadGordon, “Simulated quantum computation of molecular energies,” Science 309, 1704–1707 (2005).
 R. Barends et al., “Superconducting quantum circuits at the surface code threshold for fault tolerance,” Nature 508, 500–503 (2014).
 A. D. Corcoles et al., “Demonstration of a quantum error detection code using a square lattice of four superconducting qubits,” Nat Commun 6 (2015), 10.1038/ncomms7979.
 V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, “Signatures of majorana fermions in hybrid superconductorsemiconductor nanowire devices,” Science 336, 1003–1007 (2012).
 James D Whitfield, Jacob Biamonte, and Alán AspuruGuzik, “Simulation of electronic structure Hamiltonians using quantum computers,” Molecular Physics 109, 735–750 (2011).
 Dave Wecker, Bela Bauer, Bryan K. Clark, Matthew B. Hastings, and Matthias Troyer, “Gatecount estimates for performing quantum chemistry on small quantum computers,” Phys. Rev. A 90, 022305 (2014).
 M. B. Hastings, D. Wecker, B. Bauer, and M. Troyer, “Improving Quantum Algorithms for Quantum Chemistry,” Quant. Inf. & Comp. 15, 1–21 (2015), arXiv:1403.1539 .
 D. Poulin, M. B. Hastings, D. Wecker, N. Wiebe, A. C. Doherty, and M. Troyer, “The Trotter Step Size Required for Accurate Quantum Simulation of Quantum Chemistry,” Quant. Inf. & Comp. 15, 361 (2015), arXiv:1406.4920 .
 Ryan Babbush, Jarrod McClean, Dave Wecker, Alán AspuruGuzik, and Nathan Wiebe, “Chemical basis of trottersuzuki errors in quantum chemistry simulation,” Phys. Rev. A 91, 022311 (2015).
 D. Wecker et al., “Solving strongly correlated electron models on a quantum computer,” Phys. Rev. A 92, 062318 (2015a).
 J. Hubbard, “Electron Correlations in Narrow Energy Bands,” PNAS 276, 238–257 (1963).
 G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, ‘‘Electronic structure calculations with dynamical meanfield theory,” Rev. Mod. Phys. 78, 865–951 (2006).
 Walter Metzner and Dieter Vollhardt, “Correlated lattice fermions in dimensions,” Phys. Rev. Lett. 62, 324–327 (1989).
 Thomas Maier, Mark Jarrell, Thomas Pruschke, and Matthias H. Hettler, “Quantum cluster theories,” Rev. Mod. Phys. 77, 1027–1080 (2005).
 A Georges and G Kotliar, “Hubbard model in infinite dimensions,” Physical Review B 45, 6479–6483 (1992).
 Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg, “Dynamical meanfield theory of strongly correlated fermion systems and the limit of infinite dimensions,” Rev. Mod. Phys. 68, 13–125 (1996).
 V I Anisimov, A I Poteryaev, M A Korotin, A O Anokhin, and G Kotliar, “Firstprinciples calculations of the electronic structure and spectra of strongly correlated systems: dynamical meanfield theory,” Journal of Physics: Condensed Matter 9, 7359 (1997).
 Emanuel Gull, Andrew J. Millis, Alexander I. Lichtenstein, Alexey N. Rubtsov, Matthias Troyer, and Philipp Werner, “Continuoustime monte carlo methods for quantum impurity models,” Rev. Mod. Phys. 83, 349–404 (2011).
 Philipp Werner and Andrew J Millis, “Hybridization expansion impurity solver: General formulation and application to Kondo lattice and twoorbital models,” Physical Review B 74, 155107–13 (2006).
 Kristjan Haule, “Quantum monte carlo impurity solver for cluster dynamical meanfield theory and electronic structure calculations with adjustable cluster base,” Phys. Rev. B 75, 155113 (2007).
 Michel Caffarel and Werner Krauth, “Exact diagonalization approach to correlated fermions in infinite dimensions: Mott transition and superconductivity,” Phys. Rev. Lett. 72, 1545–1548 (1994).
 Dominika Zgid, Emanuel Gull, and Garnet KinLic Chan, “Truncated configuration interaction expansions as solvers for correlated quantum impurity models and dynamical meanfield theory,” Phys. Rev. B 86, 165128 (2012).
 Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort, “Efficient realfrequency solver for dynamical meanfield theory,” Phys. Rev. B 90, 085102 (2014).
 Ara Go and Andrew J. Millis, “Spatial correlations and the insulating phase of the high cuprates: Insights from a configurationinteractionbased solver for dynamical mean field theory,” Phys. Rev. Lett. 114, 016402 (2015).
 Daniel J. García, Karen Hallberg, and Marcelo J. Rozenberg, “Dynamical mean field theory with the density matrix renormalization group,” Phys. Rev. Lett. 93, 246403 (2004).
 F. Alexander Wolf, Ian P. McCulloch, Olivier Parcollet, and Ulrich Schollwöck, “Chebyshev matrix product state impurity solver for dynamical meanfield theory,” Phys. Rev. B 90, 115124 (2014).
 Martin Ganahl, Markus Aichhorn, Hans Gerd Evertz, Patrik Thunström, Karsten Held, and Frank Verstraete, “Efficient dmft impurity solver using realtime dynamics with matrix product states,” Phys. Rev. B 92, 155132 (2015).
 F. Alexander Wolf, Ara Go, Ian P. McCulloch, Andrew J. Millis, and Ulrich Schollwöck, “Imaginarytime matrix product state impurity solver for dynamical meanfield theory,” Phys. Rev. X 5, 041032 (2015).
 Alexei Y. Kitaev, “Quantum Measurements and the Abelian Stabilizer Problem,” (1995), arXiv:quantph/9511026 .
 Alexei Yu Kitaev, Alexander H Shen, and Mikhail N Vyalyi, Classsical and quantum computation, 47 (American Mathematical Soc., 2002).
 Qimiao Si, M. J. Rozenberg, G. Kotliar, and A. E. Ruckenstein, “Correlation induced insulator to metal transitions,” Phys. Rev. Lett. 72, 2761–2764 (1994).
 Gerald Knizia and Garnet KinLic Chan, “Density matrix embedding: A simple alternative to dynamical meanfield theory,” Phys. Rev. Lett. 109, 186404 (2012).
 Richard P Feynman, “Simulating physics with computers,” International Journal of Theoretical Physics 21, 467–488 (1982).
 Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, ManHong Yung, XiaoQi Zhou, Peter J. Love, Alán AspuruGuzik, and Jeremy L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nat Commun 5, 4213 (2014).
 Dave Wecker, Matthew B. Hastings, and Matthias Troyer, “Progress towards practical quantum variational algorithms,” Phys. Rev. A 92, 042303 (2015b).
 Michael Potthoff, “Selfenergyfunctional approach to systems of correlated electrons,” Eur. Phys. J. B 32, 429–436 (2003).
 M. Potthoff, M. Aichhorn, and C. Dahnken, “Variational cluster approach to correlated electron systems in low dimensions,” Phys. Rev. Lett. 91, 206402 (2003).
 PierreLuc DallaireDemers and Frank K. Wilhelm, “Method to efficiently simulate the thermodynamic properties of the FermiHubbard model on a quantum computer,” Phys. Rev. A 93, 032303 (2016).
 H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959).
 Masuo Suzuki, “Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to manybody problems,” Communications in Mathematical Physics 51, 183–190 (1976).
 Andrew M Childs, Richard Cleve, Enrico Deotto, Edward Farhi, Sam Gutmann, and Daniel A Spielman, “Exponential algorithmic speedup by a quantum walk,” in Proceedings of the thirtyfifth annual ACM symposium on Theory of computing (ACM, 2003) pp. 59–68.
 Anmer Daskin and Sabre Kais, “Decomposition of unitary matrices for finding quantum circuits: Application to molecular hamiltonians,” The Journal of Chemical Physics 134, 144112 (2011).
 Andrew M. Childs and Robin Kothari, “Simulating Sparse Hamiltonians with Star Decompositions,” in Theory of Quantum Computation, Communication, and Cryptography, Lecture Notes in Computer Science, Vol. 6519, edited by Wim Dam, Vivien M. Kendon, and Simone Severini (Springer Berlin Heidelberg, 2011) pp. 94–103.
 Anmer Daskin, Ananth Grama, Giorgos Kollias, and Sabre Kais, “Universal programmable quantum circuit schemes to emulate an operator,” The Journal of Chemical Physics 137, 234112 (2012).
 Dominic W Berry and Andrew M Childs, “Blackbox hamiltonian simulation and unitary implementation,” Quantum Information & Computation 12, 29–62 (2012), arXiv:0910.4157 .
 Andrew M Childs and Nathan Wiebe, “Hamiltonian simulation using linear combinations of unitary operations,” Quantum Information & Computation 12, 901–924 (2012), arXiv:1202.5822 .
 Andrew M Childs and Nathan Wiebe, “Product formulas for exponentials of commutators,” J. Math. Phys. 54, 062202 (2013), arXiv:1211.4945 .
 Hans A Bethe, “Statistical theory of superlattices,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 150, 552–575 (1935).
 See, e.g., Ref. \rev@citealpgeorges1996.