Quantum Computed Green’s Functions using a Cumulant Expansion of the Lanczos Method

.


Introduction
The Green's function (GF) is a quantity that allows access to in principle all single-particle properties, including electronic responses, and is therefore useful for calculating spectroscopic properties, such as photoemission, photoabsorption, and conductivity [1,2].It can also be used to capture electronic correlations in condensed matter and quantum chemical simulations, in addition to being the central quantity in many quantum embedding methods such as cluster perturbation theory and Dynamical Mean Field Theory (DMFT) [3,4,5,6,7].Despite its importance in these fields, the GF remains a difficult quantity Gabriel Greene-Diniz: gabriel.greene-diniz@quantinuum.com to calculate accurately from first principles on classical computers.The exact GF requires knowledge of all eigenvalues and eigenvectors of the Hamiltonian, so its calculation complexity is the same as for full Hamiltonian diagonalization.There is a long history of approximate methods for the many body GF in the classical domain of quantum chemistry, with the state-of-the art often exhibiting high-degree polynomial scaling [8,9,10].Quantum computation offers a novel and interesting route to obtain many-body GFs for strongly correlated systems.
In a previous study [11], an approach based on the Lanczos method was proposed to obtain the GF on quantum computers.This approach relies on a parameterised VQE method to construct the Krylov space, hence the orthonormal Lanczos basis is calculated iteratively in which a VQE-like algorithm is applied at each Lanczos iteration.This involves an overhead in the number of measurements and renders the approach more susceptible (relative to a typical ground state calculation) to the potential difficulties associated with variational optimization (such as barren plateaus) [12].This is also the case for recent proposal based on variational compilation [13].A method based on quantum subspace expansion (QSE) has also been proposed [14], in which the Hamiltonian and overlap matrix elements are obtained in a basis for the ground state and in another basis for the Krylov states involved in the GF.The basis sets are prepared using time evolution, and the matrix elements are subsequently measured on a quantum computer.The Krylov basis and Lanczos coefficients are then calculated iteratively on a classical computer.While an interesting approach, we note QSE requires more resources than VQE, and the reliance on trotterised time evolution may lead to a vulnerability to trotter errors as a function of the time step.A number of other recent proposals also utilise time evolution to measure the real-time GF [15,16], as well as variational methods to capture the imaginary-time GF [17,18].We also note the approach of Kosugi and Matsushita [19] in which (multi-)controlled operations between ancillas and state qubits are used to repre-sent the ladder operators involved in the (off-)diagonal elements of the GF matrix, with a subsequent application of quantum phase estimation (QPE).Other recently proposed methods obtain the GF using an algorithm to sample a Fourier series expansion [20], or using a 3-qubit iTofolli gate [21].Such approaches are prohibitive for the noisy intermediate-scale quantum (NISQ) era in terms of circuit depth scaling.Other algorithms aimed at greater suitability for near-term quantum computers have been proposed to obtain the GF in both the time and frequency domains [22], in addition to near-term algorithms for linear response properties of molecules [23] (which require ancillary qubits in some cases).
As mentioned above, the GF is a central quantity in quantum embedding methods such as DMFT [5].In recent years, a number of quantum computational approaches to DMFT have been published [24,25,26].These works involve simplifications which increase tractability, yet decrease generalisability, such as a restriction to two sites via the single Anderson impurity model (SIAM).In this picture the quasiparticle weight has an analytical form, to which the impurity-bath interaction is related in a simple way, via the limit of infinite dimensions [5].Previous works have also measured the GF in the time-domain to solve the impurity model in DMFT, which can lead to a large overhead in the number of measurements per time step [27], or a reliance on fault tolerant schemes such as QPE [28], or the use of a simplified sinusoidal form of the GF, or other techniques [25,26,29,30].A recent approach considers low energy expansions of the GF for an alternative embedding scheme based on the rotationally invariant slave boson method, which does not require calculating the full GF [31].We also note recent quantum computational approaches for the DMFT solution of the Hubbard-Holstein model [32], as well as novel spectral resolution techniques to measure the Anderson impurity GF [33].Finally, a very interesting recent approach adopted fast-forwarding circuit protocols to time-evolve the impurity state and converge the DMFT loop on noisy hardware [34].
We add to these previous works by providing a method to quantum compute the GF in the frequency domain that can be utilised in a hybrid quantumclassical algorithm for DMFT appropriate for NISQ.When considering which features of the GF are most important to capture in practice, the particular application of the GF calculation will clearly have an impact.Since the GF is a dynamical quantity, related to a wide range of real-world applications (such as the spectroscopic properties mentioned above), this motivates the calculation of frequency-dependent properties.GFs are also widely used in embedding applications.Therefore, we consider the following features of the GF in this work: accessing the frequency dependent spectral function of the GF (requiring an accurate representation of the imaginary part of the GF along the real frequency axis), and calculating the impurity GF in the DMFT embedding method (requiring an accurate calculation of the GF matrix elements along the imaginary axis to obtain convergence of the impurity problem).
In this work, a cumulant expansion of the Lanczos coefficients is utilised to provide measurable quantum circuits for each moment of the Hamiltonian [35,36].These measurable expectation values can be used to calculate the elements of the tridiagonalised Hamiltonian matrix in the Krylov basis, from which the well-known continued fraction expression of the GF can be evaluated [5,7].Thus, we extend the work of Vallury et al. [35] by utilising quantum computed moments to develop a quantum algorithm for calculating the full GF matrix.As a first example, we apply our method to the calculation of the GF of the fermionic Hubbard model.We note recent studies of the ground state properties of the Hubbard model using quantum algorithms [37,38], in addition to ground state preparation algorithms which rely on fault tolerant quantum computation [39,40,41].Rather than investigating efficient ways to find the ground state, we instead focus on a quantum approach to calculate GFs for finite chains of Hubbard sites.In addition, we show that this approach can be applied to DMFT in a straight-forward manner to iteratively calculate the quantum computed impurity GF to selfconsistency.Our methodology is implemented in the InQuanto quantum computational chemistry package [42,43].
Our results show that application of this approach to the H1-1 trapped ion quantum computer [44,45] results in spectral features (derived from the imaginary part of the GF) that compare well with ideal classical simulations, particularly in terms of peak positions, even without error mitigation techniques.Hence, expectation values of Hamiltonian moments are relatively robust to measurement noise for low lying Lanczos roots, which is consistent with previous applications of the method to infimum estimates of the ground state energy [46,35,36].The stability of this method also extends to iterative calculations of the impurity GF in the DMFT algorithm, even when (emulated) circuit based measurements are used to evaluate the Hamiltonian moments at every DMFT iteration.
This paper is organised as follows.In section 2 we describe the methods and show how the quantum computed moments can be used to obtain the GF.This is followed by a brief overview of DMFT applied to the Bethe lattice with multiple bath sites [6].In section 3, we present our results.Starting with the Hubbard model, we show that our method can be applied to multiple sites and benchmark our results against the classical Lanczos method.We then demonstrate this method on quantum hardware.Following this, we utilise the quantum computed GF in a DMFT algorithm.Section 4 provides the concluding remarks.

Quantum Computed Moments for Green's Functions
In this section, we give a brief overview of the cumulant expansion of the Lanczos method, and its application to computing the GF in a quantum computational setting.The Lanczos recursion method [47] can be viewed as an approximate diagonalisation scheme, resulting in a tridiagonal form of the Hamiltonian matrix which can be efficiently solved.In a typical classical Lanczos routine, the elements of the tridiagonalised matrix (referred to as Lanczos coefficients below) , where l = 1, 2, 3, β 0 = 0, and Ĥ is the fermionic Hamiltonian operator in second quantization.Vallury et al. [35] showed that this scheme can be re-expressed for a quantum computational context.The Lanczos coefficients can be expressed using the moments of the Hamiltonian operator Ĥ with respect to some initial Lanczos vector, ⟨ Ĥn ⟩ = ⟨v 1 | Ĥn |v 1 ⟩.Therefore, a quantum algorithm that can compute the moments can be used to generate the Lanczos coefficients.For example, for the upper 2 × 2 block of the tridiagonal matrix, i.e. α 1 , α 2 , β 1 , the quantities obtained from a quantum algorithm are the expectations values ⟨ Ĥ⟩, ⟨ Ĥ2 ⟩, and ⟨ Ĥ3 ⟩, and the coefficients are calculated as: ( In general, following the work of Hollenberg et al. [48,49], it can be shown that the Lanczos coefficients α l and β l can be expressed as functions of {⟨ Ĥn ⟩} n=1..2l−1 and explicit expressions can be obtained via the cumulant expansion of the Lanczos coefficients [35,48]. Once the Lanczos coefficients are obtained, they can be used to calculate the zero temperature GF matrix G(ω) in the continued fraction representation [5,50].In this representation, the diagonal elements of the single particle GF take the following form when de-composed into particle and hole parts [51] where E GS is the ground state energy, and the numerators n i,i (defined below) are related to normalised states obtained by operating on the ground state with the fermionic creation ( f † i ) and annihilation ( fi ) operators, applied to spin orbital (or qubit) i (Jordan-Wigner (JW) [52] encoding is assumed throughout this paper).These normalised states correspond precisely to the initial Lanczos vector for each diagonal GF element.Using the particle component g i,i (ω) as an example, we can write where n and |Ψ GS ⟩ is the ground state represented by a quantum circuit.In Eq. 2 the Lanczos coefficients α are labelled by their particle/hole contributions, according to whether the initial Lanczos vector is obtained by operating on |Ψ GS ⟩ with f † i (p) or with fi (h), respectively.
The off-diagonal elements are obtained in an analogous way using linear combinations of fermionic operators indexed by the matrix element.For example, this results in an initial Lanczos vector for g where n (p) . Following through as above, the resulting off-diagonal element G L i̸ =j (ω) is not the true off-diagonal element of the GF, but requires a modification which has a simple arithmetic form [7], shown in Eq. 5. Utilising the symmetry of the GF matrix, we finally obtain the GF off-diagonal element as . Thus, a method to quantum compute the matrix G(ω) is obtained, without the explicit reliance on multicontrolled ancilla qubits, time evolution, or phase estimation.On the other hand, relative to methods involving phase estimation, this approach in general requires a larger number of measurements due to the increase in the number of Hamiltonian moments with Lanczos root index l, hence we note trade-off between circuit depth and number of measurements when choosing between these approaches.We also note that this method does not require calling a VQE routine for each Lanczos root, nor a time-evolved subspace expansion to obtain the Lanczos basis vectors.
In practical terms, the cumulant expansion of the Lanczos coefficients combined with the continued fraction representation of G(ω) implies at least two strategies for implementation in a quantum algorithm: i) prepare the initial Lanczos vector explicitly with a state preparation circuit and measure the moments as expectations with respect to this circuit, or ii) measure each moment as the expectation value of an operator built by sandwiching Ĥn with ladder operators (or sums of ladder operators) indexed according to the initial Lanczos vector.
As an example, consider the n th moment contributing to a Lanczos coefficient required for g Using ii) we would prepare |Ψ GS ⟩ and measure Thus our methodology requires either |Ψ (p) i,i ⟩ or |Ψ GS ⟩ to be prepared on a quantum circuit, in addition to the measurement of Hamiltonian moments.While mathematically equivalent, i) and ii) can lead to considerable differences in the quantum circuit resources (depending on the state preparation methods), and catering for both approaches allows for a useful flexibility in the application of the quantum Lanczos approach to GFs.As discussed in Ref. [35], a naive counting of Hamiltonian terms results in an exponential growth in the number of Pauli strings with Hamiltonian power n.However, this can be significantly improved by the existence of large commuting sets in the Hamiltonian moment expansions [53].To mitigate the rapid growth we measured the commuting Pauli terms with a single measurement circuit.However, we note there are other efficient techniques that could be also applied based on the tensor product basis (TPB) and qubit-wise commutativity [35,54].
The following paragraphs provide further details on strategies i) and ii) for executing the quantum Lanczos approach to calculate a GF, followed by an outline of how the ground state can be prepared.

i) Initial Lanczos Vector Preparation
In the occupation number (ON) vector representation, the ground state is expressed as where x = {x i } represents a set of occupations of spin orbitals i for a given number of particles, and each basis configuration can be written as the ON vector of N i spin orbitals Given a circuit representing the ground state, the coefficients of the ON vector (Eq.8) can then be extracted by calculating the overlap (where ÛGS (θ opt ) is a unitary applied to a reference state |Ψ ref ⟩, as is typical of VQE) for all x, resulting in a set of c x values, and the corresponding {|x⟩}, hence the ON representation (Eq.8) is obtained.The ON representation of i,i ⟩ is then obtained by applying f † i or fi , respectively, to Eq. 8 and normalising (as in Eq. 3).The resulting expansion of |v 1 ⟩ as a linear combination of basis configurations with real coefficients can then be prepared on a quantum circuit using controlled Givens rotations as particle-conserving excitations [55].While the procedure defined by Eqs. 8 -10 in general exhibits an exponentially scaling overhead (due to the inclusion of all particle-number preserving configurations) for the representation of |v 1 ⟩ to be exact, it is used here to demonstrate the flexibility of our procedure to allow for arbitrary preparation of the initial Lanczos vector.Hence future preparation schemes with better scaling can be easily accommodated.

ii) Sandwiched Moment Expectation
As mentioned above, one can obtain each required term of the cumulant expansion from the measurement of an expectation, with respect to the Ĥ ground state, of a sandwiched Hamiltonian moment operator (see Eq. 7).This requires the preparation of a quantum circuit representing the ground state, which can be achieved by a VQE-optimized variational ansatz.Following the JW transformation of fi Ĥn f † i (and of f † i Ĥn fi for the hole part), the corresponding Pauli strings can then be measured with respect to the ground state circuit.In section 3.1.2,we report Pauli operators representing the measurable sandwiched Hamiltonian moments for GF elements of the Hubbard dimer for Lanczos coefficients corresponding to l = 2.

Ground State Preparation
To prepare the ground state, a number of approaches have been proposed in recent years based on hybrid quantum-classical algorithms such as imaginary time evolution [56,57], or VQE [12].These approaches generally involve the application of some ÛGS to a reference state, where ÛGS (θ) is expressed as an ansatz which depends on parameters θ that are variationally optimized to θ opt .For the latter, a wide range of chemically intuitive (e.g.unitary coupled cluster (UCC) [58]) and hardware efficient [59] ansätze have been proposed.In addition, adaptive methods also exist in which the ansatz is constructed iteratively for a given problem [60,61,62].We also note the continuing development of ground state preparation algorithms which rely on fault tolerant quantum computation [39,40,41].In this work we focus on a quantum approach to the GF and assume an accurate ground state has been provided from a separate calculation.In practice, we prepared the ground state circuit using a VQE algorithm with a parameterised ansatz, but our approach also supports the classical calculation of the ground state (which would exemplify a procedure that combines classical computing for ground state properties, with quantum computing for excited state properties).This can be done as long as the classical calculation yields an expansion of the state vector in a basis of occupation configurations, since the latter in principle can always be prepared using controlled excitation gates [55], as outlined for strategy i).

Fermionic Hubbard model
To demonstrate our approach, we quantum compute the GF of the fermionic Hubbard model, the Hamiltonian of which can be written as where µ, t, and U are the chemical potential, hopping amplitude, and on-site Coulomb interaction, respectively (t = 1 in these calculations, which sets the energy scale throughout this paper).s labels a site, <r, s> denotes nearest neighbor pairs of sites, and spins are labelled by σ = ↑, ↓.To establish consistency with the spin orbital index i of the GF matrix (see section 2.1), we use a linear mapping between the site-spin index and the spin orbital index s, ↑ (s, ↓) → i (i + 1) where s = 0, 1, 2, . . ., N s − 1 and i = 2s, hence G sσ, s ′ σ ′ ≡ G i,j .There are N s × 1 × 1 Hubbard sites arranged in a linear geometry.Since we use JW encoding throughout, the spin orbital index i also corresponds to a qubit index on a quantum circuit.Fermionic operators are mapped according to JW encoding as where i 2 = −1 (referring to the non-italic symbol i), and P i ∈ {I i , X i , Y i , Z i } refers to a Pauli rotation gate applied to qubit i (for i = 0 the product of Z d over d is omitted).This results in a qubit representation of the Hamiltonian as a sum of tensor products of Pauli operators.
Given the qubit representation of ĤHub , the corresponding GF matrix is obtained using the methodology outlined in section 2.1.We assume the halffilled regime (µ = U/2), set the initial reference state to the half-filled singlet configuration, and obtain the ground state using VQE with a parameterized ansatz.For the latter, we utilise qubit excitations [61] in the Hard-core Boson representation [63] to obtain the ground state using significantly fewer 2-qubit gates than UCC.Consider the Hubbard dimer, corresponding to 4 qubits.The initial reference state of a spin-up electron occupying site 0 and a spin-down on site 1 (|Ψ ref ⟩ = |1001⟩) can be prepared with Pauli-X gates.The ground state can then be obtained by applying two qubit excitation operators, i.e. unitary operators corresponding to a single-excitation and a doubleexcitation which act directly on the qubit Hilbert space (hence obviating the need for Pauli-Z gates to maintain fermionic exchange antisymmetry [64]) The doubles excitation portion of this circuit is shown in Fig. 1.By utilising particle-hole symmetry, as well as total spin symmetry (equal number of spin-↑ and spin-↓ electrons at half-filling), the number of 2-qubit gates for the double-excitation can be reduced as follows: step 1) treat the double-excitation initially as a singleexcitation involving Pauli operators acting on two Hubbard sites instead of four spin orbitals (equivalent to using the molecular spatial orbital index in a chemically aware strategy [63]), step 2) specify the qubit excitation between sites as the excitation between qubits representing the same spin on either site (e.g. for a spin-↑ to spin-↑ excitation in the Hubbard dimer, this corresponds to a qubit excitation involv-ing qubits 0 and 2), step 3) reintroduce spin orbital indexing by applying 2-qubit CNOTs which pair the relevant even-odd indexed qubits.The net result is a series of gates which perform the equivalent action of the double-excitation unitary e iθ2Y0X1X2X3 but with fewer CNOTs.This replaces the double-excitation circuit in Fig. 1.Due to spin symmetry and the lack of (s, ↑) → (r, ↓) cross-spin hopping between sites s and r ̸ = s, only one single-excitation unitary is needed.The resulting ansatz circuit is shown in Fig. 2. We also note that all measurable circuits in this work are compiled for hardware using the architecture agnostic quantum software compiler t|ket⟩ TM [65].
used to obtain the ground state of the Hubbard dimer, with an optimized 2-body qubit excitation.The X gate on qubit 0 corresponds to Hubbard site 0 occupied by an even-indexed electron.The S, S † , V , V † , and H gates are for rotation into the desired computational basis.The CNOTs with qubit 0 (2) as targets (controls) bracketing the Rz represent the exponentiated Pauli strings corresponding to the qubit excitations.The 3 rd and 4 th CNOTs translate the action of previous gates to odd-indexed qubits, resulting in a double excitation e iθ 2 Y 0 X 1 X 2 X 3 applied to the initial reference state |Ψ ref ⟩ = |1001⟩, followed by the single excitation e iθ 1 Y 0 X 2 .Hence in this circuit, 6 CNOTs are needed to represent the ground state.

Dynamical Mean Field Theory
We also apply our approach to quantum compute the impurity GF within the DMFT solution of the single band Hubbard model on a Bethe lattice [6].In the limit of infinite connectivity, the GF of this model can be interpreted as the impurity-site GF of an Anderson model Hamiltonian [5,66] in which Ĥimp = ĤHub (N s = 1), and baths sites are indexed by b. ϵ b,σ and V b,σ are variational parameters corresponding to the on-site bath energies and impurity-bath interactions, respectively.Here, we set ϵ b,↑ = ϵ b,↓ and V b,↑ = V b,↓ , and the Hamiltonian is constrained by particle-hole symmetry.We use the same topology for the Anderson impurity model as [27].
Solving for the eigenspectrum of ĤAnd yields the GF of the impurity+bath system, the upper 2 × 2 block of which corresponds to the impurity-site GF matrix G imp (with elements G imp i,j ).Since there is only 1 impurity site and no spin-flipping terms, impurityrelated quantities such as G imp are 2 × 2 diagonals and reduce to scalar functions of frequency in this case, however we write these as matrices in a basis of impurity spin orbitals i, j for consistency in notation and to emphasise generalisability.We solve for the G imp i,j elements using the methodology of section 2.1.
The dynamical interaction between the impurity and bath is governed by the U = 0 GF of the Anderson impurity model [6], the inverse of which is called hybridisation ∆, with elements where I i,j is an element of the identity matrix, and iω k is the k th Matsubara frequency.The Anderson model can also be related to G imp by a self-consistency con-dition [5,6] (with corresponding matrix ∆ sc ).We note that strict self-consistency can only be obtained for N bath = ∞.This limit can be numerically approximated by a finite bath by varying the bath parameters to minimise the cost function thereby fitting ∆ to ∆ sc , where | .| F denotes the Frobenius norm.Here, N ω is the number of fitting frequencies, and we note that all DMFT fitting is performed on a grid of imaginary Matsubara frequencies iω k = iπ(2k+1) where β is an inverse temperature which sets a frequency grid cutoff above iω k = 0.In our applications, we run the DMFT algorithm on 1, 2, and 3 bath sites, corresponding to 4, 6, and 8 qubits, respectively.For 1 bath site, bath fitting is performed on 26 imaginary frequencies and β = 8.While for 2 and 3 bath sites, bath fitting is performed on 64 imaginary frequencies with β = 16.Numerical minimisation of C results in a new set of bath parameters {ϵ b,σ }, {V b,σ }, which in turn define a new ĤAnd , which leads to a new G imp via the approach described in section 2.1.Eqs. 13 -16 therefore suggest an iterative algorithm, which can be terminated at a threshold value of the convergence error τ .For τ , we use the Frobenius norm of the change in ∆ sc between iterations m and m + 1, summed over N ω imaginary frequencies We assume the half-filled regime, hence µ is set to U/2 which corresponds to 1 electron on the impurity site in the ground state.To find the number of electrons in the bath, the following is performed: after finding the ground state of ĤAnd for a given set of bath parameters, the total number of electrons N e is obtained from the expectation of the total number operator.Bath sites are then filled according to the number of bath electrons = N e − 1.The resulting occupation state is then used to initialise the Lanczos procedure to solve for the impurity GF using the approach outlined in section 2.1.This procedure can be summarised by noting that we consider the Anderson impurity model to be in thermodynamic equilibrium with a reservoir of particles, hence it is treated as a grand canonical ensemble with the number of electrons in the Hubbard impurity set by particle-hole symmetry.In order to find the ground state of ĤAnd at each DMFT iteration, the UCCSD [58] ansatz is optimized using the VQE [12].

Green's Functions of Hubbard Chains
In Fig. 3 we compare a statevector simulated (noiseless) quantum GF to the classically calculated GF (throughout this work, results labelled 'classical ED' correspond to classical exact diagonalisation used to obtain the GF in the Lehmann representation at zero temperature [7]), for the non-interacting (U = 0), intermediate ( U t = 2), and strongly correlated ( U t = 8) regimes of the Hubbard dimer.In Fig. 3, the spectral function and real part of a diagonal element of the GF are plotted.The spectral function (which counts the number of states per energy normalised by π) is obtained from the imaginary part of the diagonal elements of the GF matrix where δ is a broadening term, and we set δ = 0.01.Both strategies i) and ii) (described in section 2.1) are used in Fig. 3, showing their equivalence in terms of physical results.In terms of the circuit resources used in either case, interesting differences arise due to the differences in state preparation and operator expectation measurement.In this case, the exact solution is obtained (i.e. the continued fraction GF matches the GF calculated from exact diagonalisation (ED)) for maximum l = 2.This corresponds to a maximum Hamiltonian moment of n = 3, as in Eq. 1.In Fig. 4 the spectral function is calculated for a 4-site Hubbard model, with an increasing number of Lanczos roots (i.e.increasing dimensionality of the Krylov space).This shows that at 8 Lanczos roots, reasonable accuracy is obtained in the position of spectral peaks (energies which exhibit poles of the GF), with most of the deviation relative to the exact result exhibited in the weights (rather than positions) of the peaks near ω ∼ ±2 (throughout this paper, ω is in units of t, and we set t = 1).
Figure 3: Noiseless simulations of quantum computed Green's function versus ω (in units of t) for the Hubbard dimer when U = 0 (top panels), U t = 2 (middle panels), and when U t = 8 (bottom panels), using two different initialisation strategies for the quantum Lanczos routine: blue '+' corresponds to strategy i), while red '×' represents strategy ii).In the legend, spin orbital indexes ij have been omitted and Ĥn,(p) = f Ĥn f † , Ĥn,(h) = f † Ĥn f .Left panels show the spectral function (number of states per unit energy normalised by π, defined in Eq. 18), and right panels show the real part of the G0,0 element (in units of 1  t where we use t = 1).We note that in this case the Lanczos routine reproduces the GF obtained from ED.

Strategy i) Initial Lanczos Vector Preparation: Hubbard Model
Following the noiseless VQE optimization of the ground state ansatz (Fig. 2), the ground state of the interacting Hubbard dimer exhibits the following expansion in terms of ON basis vectors where The quantum circuits representing Eqs.20 and 21 (in addition to the off-diagonal terms) can then be prepared using multicontrolled Given's rotations to obtain arbitrary (particle-conserving) linear combinations of basis states [55].In strategy i), in contrast to strategy ii), we do not measure the expectation of sandwiched moments with respect to the ground state, but rather take the expectation of the moments Ĥn with respect to the Lanczos vectors.Hence, the Pauli strings to be measured for each GF element do not change with the element index of the GF matrix.The Pauli strings take the following form after JW encoding the fermionic Hamiltonian moments (using Hubbard parameters U = 2, t = 1 as an example) (with I ≡ I 0 I 1 I 2 I 3 ).However, the circuits corresponding to the Lanczos vectors do change with GF matrix element index, and grow rapidly with the size of the system.Considering the scaling of the number of terms in Eq. 19 with respect to the number of qubits, the number of terms required to exactly represent the Lanczos vectors (Eqs.20 and 21) scales exponentially.Hence, the feasibility of strategy i) largely depends on how the Lanczos vectors are prepared on the quantum circuit.
A well known issue with the Lanczos method is the degradation of the Lanczos basis due to floating point precision [51].Typically, this degradation occurs for large values of l, where β l tends to become small and hence numerical noise is amplified when dividing by β l to normalise |v l ⟩.In classical schemes, this is typically treated by re-orthogonalising the current |v l ⟩ to the set of previously calculated vectors {|v k<l ⟩}.However, quantum noise and the statistical nature of circuit measurements required for the evaluation of Hamiltonian moments can also affect the Lanczos basis and associated Lanczos coefficients calculated from the measured moments.
To study this, the imaginary part of a diagonal element of G for the Hubbard dimer is obtained from measurements of expectations of Hamiltonian moments corresponding to the particle (g (p) 0,0 ) and hole (g (h) 0,0 ) contributions to G 0,0 .The results are obtained from the Quantinuum H1-1 emulator, a classical device emulator with a noise model corresponding to the noise profile of the H1-1 device [45].These emulated experiments correspond to one of the 14 measurable circuits to be described in section 3.2, representing the upper diagonal elements of the particle and hole GF matrices, both with 23 (7) total (2-qubit) gates and depth 13.These are plotted alongside the normalised absolute errors in α l and β l resulting from measurements, and the diagonal and off-diagonal components of the overlap matrix elements calculated classically from the Lanczos vectors which in turn are built from the measured α l and β l .Since the ED result for the Hubbard dimer GF is recovered for maximum l = 2, only α 1 , α 2 , and β 1 are used in this case, and we average the error in α l .Hence the errors in the Lanczos coefficients are calculated as where α l,exact are noiseless ideal values of the Lanczos coefficients.The next Lanczos vector |v 2 ⟩ is then constructed from these measured Lanczos coefficients, and its norm ⟨v 2 |v 2 ⟩ and overlap with v 1 are obtained classically (since the ground state parameters are obtained from an ideal noise free calculation, and the vector norms are calculated classically, the first Lanczos vector v ⟩ has unity norm by construction).The results of 4 separate emulator runs, corresponding to 4 separate sets of measurements, are shown in Fig. 5.It can be observed that, for a given set of measurements, the error in β l resulting from the quantum noise correlates to the deviation from unity norm of v 2 , which is not unexpected since β l governs the normalisation of the Lanczos vectors.This error in normalisation also translates to the deviation of spectral peak heights from the exact result, which gives a qualitative understanding of the effect of quantum noise in β l on the resulting GF: within each plot of Fig. 5, a larger |∆β (p/h) | results in a worse relative peak height (compared to ED) for g (p/h) 0,0 , where superscript "p" ("h") denotes contribution to the positive (negative) frequency spectral peaks.
Also, quantum noise at the circuit level can have a cumulative effect on the orthogonalisation of the Lanczos vectors at higher roots, the study of which requires measurements of corresponding higher moments of the Hamiltonian with respect to the ground state for models larger than the Hubbard dimer.Such measurements require circuits too deep for the current NISQ era, as errors due to gate infidelities and qubit decoherence would accumulate and dominate the measurements of associated Pauli operators.Hence, we leave the effect of quantum noise on the orthogonalisation of Lanczos vectors for higher lying roots as an open question for future studies, although we note the investigation of Vallury et.al. [35] into the effect of noise and device errors on the arithmeti-cal operations involved in the assembly of cumulants, in which it was found that infimum estimates of the ground state energy are improved in accuracy and more robust to device noise relative to variational ap-proaches to optimizing ⟨ Ĥ(θ)⟩ [36].The relation between numerical errors in the Lanczos basis and quantum noise could have interesting implications for error mitigation, which we briefly discuss in the conclusion.

Strategy ii) Sandwiched Moment Expectation: Hubbard Model
Using the particle GF with U = 2, t = 1 as an example, the first, second, and third powers of the Hamiltonian, sandwiched between ladder operators, and contributing to the diagonal element g (p) 0,0 , are mapped to the following sums of Pauli operator strings via JW encoding.
At variance to strategy i), in strategy ii) the measurable Pauli strings contributing to each GF element change with matrix element index.This is exemplified for 2, 3, and 4 Hubbard sites (4, 6, and 8 qubits, respectively) in Fig. 6, in which the number of Pauli strings in f0 Ĥn f † 0 , and in ( f0 + f2 ) Ĥn ( f † 0 + f † 2 ) (contributing to the off-diagonal element g (p) 0,2 ), are plotted as a function of n.
Interestingly, the total number of individual Pauli terms does not necessarily increase with the Hamiltonian power.This effect manifests in two ways: 1) The concatenation of Pauli strings that occurs when taking the power can result in the collapse of products into smaller equivalent strings; this is evident for smaller numbers of qubits where the number of Pauli strings can oscillate for even and odd powers of Ĥ, however the importance of this decreases rapidly as the system size grows (see Fig. 6).2) In general, the maximum number of Pauli strings for N qubits is 4 N , however the Hamiltonian will have less terms than this (depending on the interactions of the model) for a given N which will affect the number of individual strings that result from n powers of Ĥ; our results show that the number of Pauli strings for the Hubbard Hamiltonian saturates for values of n large enough for the Lanczos procedure to sufficiently cover the Hilbert (sub)space of the model (in general of lower dimen-sionality of the full Hilbert space of N qubits due to symmetry), which occurs at large enough values of the Lanczos index l defined in section 2.1.
The resulting manageable number of Pauli terms of small models, in addition to the optimized qubit excitation ansatz circuit for the ground state presented in section 2.2 (in this strategy, the circuits corresponding to the bra and ket of ⟨ Ĥn,(p/h) ii ⟩ do not change with GF matrix index, unlike strategy i)), facilitate the quantum calculation of the full GF matrix of the Hubbard dimer on Quantinuum's H1 ion-trap machine, the results of which are presented in the next section.3.2 Green's Functions Calculated from Hardware Experiments

Hubbard Model
In Fig. 7, the spectral function of the Hubbard dimer is obtained from the H1-1 trapped ion quantum computer.In addition, results from the Quantinuum H1-1 emulator (H1-1E) (with its noise model calibrated to the noise profile of H1-1 [45]) are also shown.Both the particle and hole GF matrices contain 16 elements, which immediately reduces to 10 elements each due to symmetry about the diagonal [7].This would naively result in 20 measurable circuits for the full GF (10 for both particle and hole GFs): one circuit per matrix element for particle and hole matrices (despite particle-hole symmetry, we explicitly calculate both the particle and hole matrices to more comprehensively test the performance on hardware).Due to certain elements within the particle and hole GF matrices being identical for the Hubbard dimer in this regime, the final total number of measurable circuits for the Hubbard dimer is 14.The total number of (2-qubit ZZ) gates per circuit ranges from 23 to 27 (7 to 9), with depth ranging from 13 to 18. Measurements of each circuit are performed with 8192 shots.Excellent agreement is observed in the spectral peak positions, indicating an accurate description of the excitation energies involved in the particle/hole transitions, when compared to the ideal classical simulation.Previous investigations have demonstrated the robustness of this method to quantum measurement noise and device errors for infimum estimates of the ground state energy [46,35,36], while other recent studies have also shown the robustness to noise of quantum subspace algorithms in general for accessing low lying eigenvalues [67,68].Our results are consistent with these previous works.The quality of the results is also indicative of low device errors on the H1-1 trapped ion quantum computer, since measurements of Hamiltonian moment expectations were performed without error mitigation.We also note the excellent agreement between H1-1 and its emulator.
The heights of the two central peaks in Fig. 7 are underestimated with respect to ED by an average (between positive and negative frequency peaks) of approximately 17% for H1-1.Repeated runs on the emulator H1-1E shows this is roughly consistent (Fig. 8), yielding a mean underestimation of 15% over the 5 runs for the two central peaks.With regard to the peak positions along the real frequency axis, the hardware results from H1-1 (Fig. 7) yield an average error relative to ED in the low lying peak frequency (at absolute frequency |ω| = 1.236 in the exact GF) of 0.5% (averaged over positive and negative frequency peaks), whereas the higher lying peaks (at |ω| = 3.238 in the exact GF) exhibit a frequency-averaged error of 1.3%.The repeated runs of the emulator (Fig. 8) yield a roughly similar average error of 0.6% in the low lying peak positions, and an average higher lying peak position error of 1.5%.
We note that the H1-1 device currently exhibits typical one-qubit and two-qubit gate error rates of 4×10 −5 and 2×10 −3 , respectively, and a typical state preparation and measurement error rate of 3 × 10 −3 [45].With this in mind, and to further investigate the impact of noise on our results, we performed bootstrapping of H1-1 hardware measurements to generate ensembles of GFs resampled from the original measured results.Statistical analysis of these bootstrapped samples provide effective error bars on the position and weight of peaks of the spectral function, thereby demonstrating the impact of random variations (e.g.due to device error) on the final output GF along the real frequency axis.
To perform such an analysis of the hardwarecalculated GF, the local spectral maxima A(ω peak ) were sampled at the transition frequencies (ω peak , of which there are 4 for the Hubbard dimer, located at frequencies ω peak = −3.238,−1.236, 1.236, 3.238 in the exact GF) for each bootstrapped GF.Each of these local maxima were then averaged over the ensemble, with their spread and standard deviations also obtained.A similar procedure was carried out for the positions of spectral peaks along the real frequency axis (the ω value for which A(ω) has a peak near ω = ω peak ), sampling the frequencies at which the spectrum has local maxima for each of the bootstrapped GFs.For the real part of the GF matrix element (right panel of Fig. 10), the value of the positive and negative extrema centred at the ω peak values were sampled over the ensemble.Fig. 9 shows the variation of spectral peak positions and heights with respect to the number of samples, showing reasonable convergence of the statistics at 500 bootstrapped samples.The ensemble of 500 bootstrapped GFs is plotted in Fig. 10.We observe that the low lying transition energies are well reproduced with respect to the exact result, with very small variation (standard deviation approximately 0.2% -0.4% normalised to the ensemble average at 500 samples as shown in Fig. 9) between bootstrapped samples.For the higher lying spectral peaks (ω peak ≈ ±3.238), standard deviations of the transition energy are larger (approximately 0.8% -1% at 500 samples in Fig. 9) which is unsurprising given the diminished accuracy of higher eigenspectra for truncated Krylov spaces, and is consistent with the results shown in Fig. 5.
In section 3.1.1we discussed how noise in the measured moments translate to errors in the normalisation of Lanczos vectors, which in turn translates to errors in the peak heights.Here we see that these errors exceed (in proportion) errors in the peak position frequencies (this is reflected in the larger average error relative to ED for the peak heights, in addition to the larger spread of spectral weights, compared to the transition frequencies).Hence variations due to device error seem to have a larger impact on spectral weights compared to errors in the values of transition frequencies.However spectral weights of low and high lying peaks are qualitatively in proportion with respect to exact diagonalisation: the dominance of low lying peaks over higher lying excitations is maintained even in the presence of quantum noise.This combined with the relatively accurate peak positions empirically demonstrates that our quantum Lanczos approach can reproduce important frequency dependent features of the GF in the presence of noise on the H1-1 quantum computer.
Whether the observed errors are large enough to prevent application of this quantum GF approach to more elaborate simulations, such as quantum embedding within DMFT, remains an interesting question.In the next section, we address this issue and demonstrate the usefulness of this technique for DMFT.

DMFT
The quantum computed moments approach is also used to obtain the GF of the impurity site in a DMFT algorithm.Running the DMFT procedure with the impurity GF (corresponding to 1 bath site) obtained from quantum computed moments at the ideal, noiseless statevector level reproduces the ED result (in which the impurity GF is calculated using classical exact diagonalisation at each DMFT iteration) exactly, corresponding to the black line in Fig. 11.Also in Fig. 11 the impurity GF corresponding to 1 bath site is quantum computed at the final DMFT step in which the Hamiltonian corresponds to converged bath parameters (where convergence was obtained using ideal noiseless simulations of the quantum GF at each DMFT iteration), using measurements of expectations of Hamiltonian moments performed on the H1-1 trapped ion quantum computer (green circles correspond to real hardware [45] in Fig. 11).In this case 2 Lanczos roots were used, which reproduces ED (consistent with the Hubbard dimer).For 1 bath site, 14 circuits were measured (corresponding to 14 impu-rity+bath GF matrix elements, similar to the Hubbard dimer), ranging from 58 (22) to 67 (25) gates (2-qubit gates), with depths ranging from 44 to 51.The larger number of gates compared to the Hubbard dimer is a result of the reduced symmetry of the impurity-bath Hamiltonian for DMFT.Each circuit is measured with 8192 shots.
Excellent agreement is observed between the hardware, emulator, and noiseless simulations, and we comment here on the observed robustness to error.At 8192 shots per circuit the absolute error in the expectation values of Hamiltonian moments due to shot noise (i.e.only from measurement sampling) is approximately 0.0062 for the first and second power, and 0.0127 for the third power (2 Lanczos roots requires up to the third power of the Hamiltonian).However, the GF matrix elements involve continued fractions of polynomials of Pauli strings which can lead to partial error cancellation when grouping Pauli terms into commuting sets: once the Pauli commuting sets are found, the measurement sampling algorithm [69] can associate certain Pauli operator expectations with both α l and β l , hence errors in those Pauli expectations can cancel when evaluating the α l and β l terms.We note that this is an interesting effect related to the combination of Krylov basis construction using Pauli-encoded Hamiltonian moments along with elements of the resulting tridiagonal matrix to build a Green's functions via continued fractions (the latter being the novelty of this work).
To investigate the impact of noise on DMFT convergence of the impurity GF, the final impurity GF obtained from the quantum computed moments method was calculated using a simple noise model and compared to the ideal noiseless result.To this end, bath parameters were first obtained by converging the im- purity GF at the ideal statevector level to τ = 0. Using these bath parameters, the impurity GF was then re-computed at the circuit level using a simulated noiseless backend [70] with 8192 shots per circuit, and with a range of values of a noise parameter λ.Focusing on two-qubit gate errors, our simplified noise model introduces a depolarising channel to all two-qubit gates in the circuit, where each two-qubit gate error can be expressed by the error channel E as which maps a corresponding state ρ to its linear combination with a maximally mixed state.This depolarising channel yields a corresponding error in the impurity GF.To quantify the latter, we re-use the threshold for DMFT convergence, τ , defined in Eq. 17.We calculate the additional error in the impurity problem due to depolarising noise as the value of τ resulting from comparing the ideal noiseless (statevector) converged ∆ sc (related to the impurity GF by the Anderson model self-consistency condition, see Eq. 15), ∆ sc ideal , to its value obtained in the presence of depolarising error channel E (with the same ground state and bath parameters), labelled as ∆ sc noise .
The error in ∆ sc as a function of λ is plotted in Fig. 12.A non-zero τ E at λ = 0 is due to finite sampling of the measurement distributions.We observe that in this range of λ, τ E increases monotonically with λ, and for small λ (≲ 0.01) reasonably small errors in the impurity GF matrix on the order of 10 −3 − 10 −4 are obtained.We also note that for this impurity GF errors of this magnitude do not change the number of GF poles and maintain the dominance of the two central peaks of the spectral function, as was observed in Fig. 11.While not directly equivalent to the two-qubit gate infidelity of the H1-1 device (2 × 10 −3 [45]), we expect two-qubit error rates on H1-1 to correspond to λ well below 0.01, which could translate to the ability to converge the DMFT procedure to correspondingly small errors in the impurity GF.As a step further, and to test this hypothesis, a hardware emulator in which the noise model is calibrated to H1-1 was also used to quantum compute the GF at each iteration of DMFT, for 20 iterations.The value of τ at each DMFT iteration is plotted in Fig. 13.Starting with random bath parameters, the initial error in ∆ sc is τ = 0.00865.After 20 iterations this error reduced to τ = 0.000237, a value sufficiently low as to obtain a reasonably accurate solution to the impurity problem.We also notice that this value of τ is slightly less than the λ = 0 value of τ E shown in Fig. 12 (0.000254), indicating that the impurity GF error between DMFT iterations near convergence performed on the H1-1 emulator does not exceed the error (relative to the ideal noiseless DMFT result) introduced by measurement sampling errors (without noise), for this case.The resulting spectral function is shown in Fig. 14.Hence the error due to quantum noise from the H1-1 emulator was sufficiently low to run the DMFT algorithm applied to 1 impurity site and 1 bath site for multiple iterations towards convergence.We also compute the impurity GF corresponding to converged bath parameters for 2 bath sites in the DMFT algorithm, in which 8 Lanczos roots are calculated.This is estimated to be reasonably accurate for 2 and 3 bath sites, in accordance with Fig. 4 in which it was shown that 8 Lanczos roots provide a good approximation to exact diagonalisation for 4 Hubbard sites.
Results for 2 bath sites are shown in Fig. 15.Initially, the DMFT algorithm was run in an ideal fashion, in  which the impurity GF was calculated using a noiseless classical (statevector) simulation of the quantum computed moments.In this setting, full DMFT convergence (τ = 0) was achieved after 14 DMFT iterations.This is shown by the green line in Fig. 15.
For comparison, a DMFT run was also performed in which classical ED was used to obtain the impurity GF at each DMFT iteration (again achieving convergence after 14 iterations), shown by the black line in Fig. 15.This compares well with the noiseless approach for 8 Lanczos roots.Following this, the H1-1 emulator (H1-1E) [45] was used to quantum compute the impurity GF using the statevector-converged bath parameters (blue squares in Fig. 15).For 2 bath sites, 65 circuits were measured (emulated hardware) at the final DMFT iteration to obtain the impurity+bath GF, ranging from 193 (81) to 204 (86) total gates (2qubit gates), with depths of 212 to 214.The spectral function is shown in Fig. 15.In this case 3 separate applications of the quantum method are applied to the emulator H1-1E, in order to simulate 3 separate hardware runs and assess the variation of results.We observe that the energy positions of low lying spectral peaks generally agree well with the exact result, and are also stable against statistical variations in measurements.Higher lying spectral peaks differ in position from the exact result more significantly, while peak heights (associated with normalisation of Lanczos vectors) are also less accurate and are more sensitive to variations due to measurement statistics.Finally, ideal simulations of this approach, in which the Hamiltonian moment expectations are evaluated in an noiseless statevector fashion at all DMFT iterations, are performed for 3 bath sites.As in the previous case, we also ran a variant of the DMFT algorithm in which the impurity GF is obtained from classical ED at each DMFT iteration.In both the quantum (statevector) and classical ED runs for 3 bath sites, the impurity GF was iterated to convergence (τ = 0) after 14 DMFT iterations.The resulting spectral functions are shown in Fig. 16.Disagreements arise in the total number of spectral peaks, and in the weight of low energy poles (consistent with observations made for the spectral function of the 4-site Hubbard model, see section 3.1 and Fig. 4).In this case, disagreements with classical ED are due to the continued fraction representation of the GF in Lanczos method; While the Lehmann representation of the GF used in the classical ED is exact for all temperatures T including the limit T → 0, this is not the case for the Lanczos calculated GF.In the latter, the exponential factor containing the inverse temperature β is discarded.Hence, despite the common practice, it is strictly an approximation to use a finite β in the fitting of the Lanczos-calculated impurity GF.However, this practice is justified by the common observation that the discrepancies between the Lanczos and ED impurity GFs vanish at sufficiently high β (sufficiently low T).We note that full DMFT convergence is achieved for the quantum Lanczos computed impurity GF, as well as for the classical ED GF, hence the DMFT algorithm finds slightly different solutions corresponding to the differences arising from approximations in the Lanczos computed GF.We choose a relatively low value of β for quicker convergence of the impurity GF, however we expect the differences between Lanczos and ED impurity GFs to vanish as β → ∞.Hence, discrepancies here are not due to errors from quantum noise or measurements, and these results show that this quantum approach in principle works for multiple bath sites in DMFT.

Conclusions
In this paper, a quantum computational approach has been presented to calculate the single particle Green's function in a spin orbital basis, using a cumulant expansion of the Lanczos procedure involving quantum computed moments.This is implemented in the In-Quanto package [42,43], and is an extension of a recent work which presented the quantum computed moments approach to obtain infimum estimates of the ground state energy [46,35,36].In our work, we utilise quantum computed moments to obtain the Lanczos coefficients, which in turn can be used to obtain the Green's function in the continued fraction representation.
Following a brief outline of the method, we show that our implementation allows for multiple strategies to initialise and run the Lanczos procedure; we present two separate strategies, one involving the explicit preparation of the first Lanczos vector on a quantum circuit (strategy i)), the other involving measurements of Pauli terms representing the Hamiltonian moments sandwiched between ladder operators indexed by the corresponding Green's function matrix element (strategy ii)).While strategy i) allows for a flexibility in the choice of representation of the Lanczos vector, we focus on strategy ii) for application of the method on trapped ion quantum computers due to the "NISQ-friendly" number of measurable Pauli terms for small sized Hubbard models.
Using our approach, we computed the GF matrix of the Hubbard dimer on the H1-1 trapped ion quantum computer, showing excellent agreement with the ideal noiseless result in terms of spectral peak positions.We also note the good agreement between H1-1 and H1-1E in Fig. 7, which also indicates the accurate approximation of hardware results when applying the emulator to larger models in section 2.3.Following the GF of the Hubbard Hamiltonian, we then apply our approach to obtain the impurity GF in a DMFT algorithm with up to 3 bath sites.Hardware results again show good agreement with the ideal noiseless result when applied to the final DMFT iteration, and emulated hardware results indicate that errors in the GF due to quantum noise do not prevent convergence of the DMFT algorithm when the GF is quantum computed at all DMFT iterations.
We note that our approach does not require ancillary qubits (only the ground state or Lanczos vector state circuits, and the Hamiltonian as a sum of Pauli operators).We also emphasise that no error mitigation has been applied to the hardware or emulator results in this paper.This is in contrast to a recent interesting work which utilised error mitigation to calculate linear response properties of molecules, and which would require ancillas for non-symmetric cases [23].In addition, another recent work also investigated a Lanczos recursion method to obtain the GF from quantum computers [71].In the latter, a state-preserving quantum counting algorithm is used to find the Lanczos coefficients.This quantum counting algorithm [71] relies on ancillary qubits and a QPE block, which likely render this approach difficult for near-term applications.This is in contrast to our approach which utilizes a cumulant expansion of the Lanczos coefficients in terms of the Hamiltonian moments, rather than a counting algorithm which requires QPE.
We also studied the scaling of the number of measurable Pauli terms with respect to the Hamiltonian moment index, for the sandwiched moment operators required for strategy ii), indicating polynomial scaling and a saturation for high-lying moments.Finally, we mention our investigation of errors in the Lanczos basis and corresponding coefficients arising from quantum noise, and how these can impact elements of the GF matrix.As mentioned in section 3.1.1,an accurate description of the relation between quantum noise and errors in the Lanczos basis could lead to very useful techniques to correct for these errors.For example, knowledge of this relation could be used to design error mitigation protocols in which noisy measurements of Hamiltonian moments, or the resulting errors in the Lanczos coefficients, could be corrected/mitigated by known calibration data of a particular device.This could potentially widen the application domain of quantum computed Green's functions by allowing for larger system sizes, represented by larger, more error prone circuits.We consider this a fruitful direction for future work.
Figure 1: Quantum circuit representation of the double qubit excitation operator e iθ 2 Y 0 X 1 X 2 X 3 .The V , V † , and H (Hadamard) gates are for rotation into the desired computational basis (such that HZH = X and V † ZV = Y ), while the CNOTs and Rz represent the exponentiated Pauli strings corresponding to the qubit excitation.Here 6 CNOTs are needed for the double excitation (2 extra CNOTs are also needed for the singles excitation not shown here, so 8 CNOTs altogether to express the ground state).

Figure 4 :
Figure 4: Spectral function of the the 4-site Hubbard model when U t = 2, using noiseless simulations of quantum computed Green's functions, using strategy ii) for the quantum Lanczos routine.Convergence towards the ED result of the quantum computed spectral function with respect to the number of Lanczos roots (maximum value of l) is observed.Peak positions are reasonably converged for ≥ 8 Lanczos roots.Dashed black line shows the result from exact diagonalisation.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

⟨
the coefficients depend on U and t Hubbard parameters, with c 1100 = c 0011 , c 1001 = −c 0110 due to symmetry and fermionic exchange, and for the noninteracting U = 0, t = 1 case, |c 1100 | = |c 0011 | = |c 1001 | = |c 0110 | = 1/2.Taking the diagonal element of the particle and hole GFs g Ĥn,(p) ii ⟩ and ⟨ Ĥn,(h) ii ⟩, respectively, the following expansions for these Lanczos vectors are obtained after applying the ladder operators to the ground state |Ψ

Figure 5 :
Figure 5: Imaginary part of G(ω)0,0 obtained from the Quantinuum H1-1 emulator, from 4 separate runs.For each run, the associated errors in Lanczos coefficients are shown, along with the norms of the Lanczos vector constructed from the (moment-)measured coefficients.

Figure 6 :
Figure 6: Number of Pauli strings in sandwiched moment operators (for diagonal and off-diagonal elements) versus Hamiltonian moment index n, using 4, 6, and 8 qubits, for non-interacting U = 0 (left panels) and interacting | U t | = 2 Hubbard models.For an interacting case of 4 qubits, oscillations in the number of Pauli strings are observed for even and odd n; the magnitude of these oscillations rapidly decrease as the system size increases.For larger numbers of qubits, the shape of the curve (semilog plot) is characteristic of polynomial scaling, and indicates a roughly convergent number of Pauli strings as a function of n as the Lanczos procedure approaches the full dimensionality of the Hilbert space (or symmetry-reduced subspace thereof).For 8 qubits, dotted lines are shown to indicate saturated numbers of Pauli strings.

Figure 7 :
Figure 7: Quantum computed GF obtained from the H1-1 trapped ion quantum computer.Solid black line shows the result from classical ED.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 8 :
Figure 8: Quantum computed GF obtained from the classical emulator of the H1-1 trapped ion quantum computer, with a noise model calibrated to current hardware data.Solid black line shows the result from classical ED.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 9 :
Figure 9: Analysis of bootstrapped GFs derived from measurements obtained on the H1-1 trapped ion quantum computer, showing the variation of results with respect to the number of bootstrap samples.Top left panel shows spectral peak heights averaged over each number of samples ( Ā(ω peak )), and bottom left panel shows the peak height standard deviations.Top right panel shows transition frequency (or energy) positions of the 4 peaks of the Hubbard dimer spectral function averaged over the number of samples ( ωpeak )), and bottom right panel shows the peak position standard deviations.All standard deviations σ here are expressed as a percentage normalised to the ensemble average for a given number of samples %σ(X, n sample ) = σ(X, n sample )/ X(n sample ) × 100, where X = A(ω peak ), ω peak .Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 10 :
Figure 10: Bootstrapped spectral functions derived from quantum computed GF corresponding to 500 samples, with original measurements obtained from the H1-1 trapped ion quantum computer.On the left panel, blue (red) error bars denote the standard deviations (minimum-maximum spreads) of Aω peak , centred on the mean, while green (orange) error bars denote the standard deviations (minimum-maximum spreads) of ω peak values.All error bars are centred on the ensemble mean.The right panel shows the corresponding information for the real part of the G0,0 matrix element (in units of 1 t where we use t = 1).Solid black line shows A(ω) from classical ED.In this figure, all standard deviations are in units of the corresponding quantity, i.e. blue bars represent units of A(ω) while green bars represent units of ω.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 11 :
Figure 11: Quantum computed impurity GF, in which the H1-1 trapped ion quantum computer is used to compute the impurity GF at the final DMFT iteration, following classical impurity GF computations for the previous iterations.Solid black line shows the final impurity GF obtained from classical ED.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 12 :
Figure 12: Error in ∆ sc , as defined in Eq.31, versus twoqubit gate depolarising error.At λ = 0 corresponding to zero noise, τE = 0.000254 which is close to the value of τ obtained from the H1-1 emulator after 20 DMFT iterations.

Figure 13 :
Figure13: Error in ∆ sc , as defined in Eq. 17 (which sets the threshold for DMFT convergence) at each DMFT iteration, in which the H1-1E emulator is used to quantum compute the impurity GF at each iteration.

Figure 14 :
Figure 14: Quantum computed impurity GF after 20 DMFT iterations, in which the emulator of the H1-1 trapped ion quantum computer is used to compute the impurity GF at each iteration.Black line shows the final impurity GF obtained from classical ED.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 15 :
Figure15: Quantum computed impurity GF of DMFT with 1 impurity and 2 bath sites, in which the emulator of the H1-1 trapped ion quantum computer is used to compute the impurity GF at the final DMFT iteration, following statevector computations for the previous iterations.Solid green line shows the converged impurity GF from ideal noiseless quantum computed moments.Solid black line shows the converged impurity GF when classical ED is used at each DMFT iteration.The emulator H1-1E is applied in 3 separate calculations of the quantum computed moments impurity GF (using bath parameters obtained from the noiseless run), to simulate 3 separate hardware runs and assess variation of results.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.

Figure 16 :
Figure 16: Quantum computed impurity GF for 3 bath sites, calculated to full convergence of the DMFT algorithm.Solid green line shows the converged impurity GF from ideal quantum computed moments.Solid black line shows the converged impurity GF when classical ED is used at each DMFT iteration.Dimension of spectral function (defined in Eq. 18) is number of states per unit energy normalised by π.