One-particle Green's functions from the quantum equation of motion algorithm

Many-body Green's functions encode all the properties and excitations of interacting electrons. While these are challenging to be evaluated accurately on a classical computer, recent efforts have been directed towards finding quantum algorithms that may provide a quantum advantage for this task, exploiting architectures that will become available in the near future. In this work we introduce a novel near-term quantum algorithm for computing one-particle Green's functions via their Lehmann representation. The method is based on a generalization of the quantum equation of motion algorithm that gives access to the charged excitations of the system. We demonstrate the validity of the present proposal by computing the Green's function of a two-site Fermi-Hubbard model on a IBM quantum processor.


I. INTRODUCTION
The exponential scaling of time and memory resources required to accurately simulate an interacting quantum system of many particles poses fundamental limits on what can be achieved by classical algorithms [1,2]. This computational burden lies at the heart of many open problems in condensed matter physics, such as the description of high-temperature superconductors [3] and frustrated magnetic materials [4], for which an efficient and reliable computational treatment is still missing. In this context, quantum computing has recently emerged as a new promising paradigm [5] that can potentially outperform classical counterparts on selected tasks [6], including the simulation of quantum systems [7][8][9][10]. In the near term, quantum computers may reliably control only a relatively small amount of qubits and perform a limited set of operations on them before significantly loosing coherence [11]. Nevertheless, these devices are quickly approaching the stage at which it becomes impossible to simulate their behavior classically [12]. This motivates the search for quantum algorithms that can effectively exploit such limited but potentially useful computational power. To this aim, one should in general design methods which require a relatively low qubit and gate count, but that can still offer a quantum advantage by accessing otherwise classically intractable regimes [13,14].
The calculation of static and dynamical properties of many-body quantum systems is in fact one of the most promising research areas in which quantum computers could offer significant advantages. In this context, several near-term methods have been proposed, starting from the well known variational quantum eigensolver (VQE) for calculating the ground state of a given Hamiltonian [14,15]. Quite generally, these methods optimize the * Corresponding author. jrizzo@sissa.it use of quantum resources by splitting the target problems into parts, distinguishing between those that can be efficiently solved on a classical computer and the ones [16,17] that instead should be delegated to quantum processors. Promising quantum solutions have also been proposed and implemented for the evaluation of the excitation spectra [18][19][20][21][22][23], for variational time evolution [24] and, more recently, to compute dynamic correlation functions [25][26][27][28][29]. The latter are of crucial importance for the study of optical, magnetic, and transport properties of many-body quantum systems [30]. As an alternative to approaches based on the many-electron wavefunction, Green's function techniques can provide solutions for the calculation of interesting physical properties, which in the frequency domain can be directly accessed and validated in photoemission experiments.
In this work, we propose a novel near-term quantum algorithm for computing the Green's function via its Lehmann representation, following the general setup presented by Endo et al. in Ref. [28]. Our method relies on an original extension of the quantum equation of motion (qEOM) algorithm [23] that allows us to compute the charged excited states of the system. The latter are used to efficiently compute the spectral amplitudes needed to evaluate the Green's function in the frequency domain. We test our proposed quantum algorithm numerically and experimentally on IBM Quantum devices by computing the Green's function of a two-site Fermi-Hubbard model.
The rest of this paper is organized as follows. In Sec. II we fix the notation for the formalism of many-body Green's functions. In Sec. III we describe our proposed method for computing the Green's function on a quantum computer via the qEOM technique. In Sec. IV we demonstrate the validity of our proposal by performing numerical and hardware simulations. Finally, in Sec. V we conclude by commenting on possible future improvements and research directions.

II. MANY-BODY GREEN'S FUNCTIONS
We start by reviewing the formalism of many-body Green's functions at zero temperature for fermionic quantum systems [31]. We consider a system of N identical fermions defined by its HamiltonianĤ, projected onto a basis of M fermionic modes {|ψ α } M α=1 . In this context we may introduce the fermionic creation (destruction) operatorsĉ † α (ĉ α ) satisfying the usual anti-commuting algebra, and we can define the single-particle retarded Green's function (GF) at zero temperature as [30] Here θ(t) is the Heaviside step function,ĉ α (t) = e iĤtĉ α e −iĤt is the Heisenberg evolution of the operator c α and the expectation value is computed over the N particle ground state of the Hamiltonian |ψ 0 . In our definition we take = 1. The GF as defined in Eq. (1) may be expressed in Fourier space by introducing the regularization factor lim η→0 + e −η|t| as The real frequency GF then becomes [30] If we insert in Eq. (3) the representation of the identity I = n |E n E n |, where |E n are the (orthonormal) eigenstates of the system, Eq. (3) can be written as where ) are the eigenstates of the system with N + 1 (N − 1) particles, also referred as the particle (hole) spectrum. Eq. (4) constitutes the Lehmann representation of the retarded GF. In what follows, for simplicity, we restrict our attention to the diagonal terms of the GF In fact, the generalization to Eq. (4) of the method we propose can be easily accomplished. The single-particle GF in the real frequency domain is an interesting quantity, since its imaginary part, the spectral function, is directly accessible via photoemission spectroscopy [30] A α (ω) = −π −1 Im(G R α (ω)). (6) In general, the single-particle GF characterizes the dynamics of quasiparticles in correlated systems [31]. Furthermore, other types of GF can be used to characterize the optical, magnetic and transport properties of a system in the framework of linear response theory [30], where the effect of an external perturbationBV (t) on an observableÂ can be expressed as where The classical equation of motion (EOM) method was initially introduced by Rowe in Ref. [32] as a way to obtain excitation energies and excited states of a given HamiltonianĤ. In the EOM approach we search for an approximation to the excitation operatorÔ † n = |n 0|, which generates the n-th excited state |n of the Hamiltonian starting from a given generic state |0 with non-zero overlap with the true ground state |0 . The excitation operator satisfies where E 0n = E n −E 0 are the exact excitation energies. If we take |0 = |0 , where |0 is the true ground state of the system, the annihilation conditionÔ n |0 = 0 is satisfied, and by operating from the left-hand side of Eq. (9) witĥ O † n |0 we may write This equation holds if the excitation operator is exact and if |0 is the true ground state. Using again the annihilation condition, we can equivalently write if the annihilation condition is satisfied. This second formulation has the advantage of guaranteeing real-valued energy differences, since the operator [Ô n ,Ĥ,Ô † n ] in Eq. (11) is Hermitian. We can find approximate solutions to Eq. (11) by considering it as functional of the excitation operatorÔ n that should be minimized. In fact, by taking a variation δÔ n and imposing the stationary condition δE 0n = 0 we get again Eq. (9), namely the functional of Eq. (11) is stationary when evaluated at the exact excitation operators. We can then approximate the excitation operators by linearly expanding them in a basis of elementary excitations. For example, if one wants to target the charged excitations of the system one may consider the ansatẑ where µ α is a collective index of the single-particle orbitals entering in the charged excitation operatorÊ By imposing the stationary condition on Eq. (11) after inserting Eq. (12) we obtain the generalized eigenvalue problem (GEP) equation where In this case, if the ground state has N particles one targets the sector with N + 1 particles. We note that, since the functional of Eq. (11) can be intended both forÔ † n andÔ n , one also gets the excitation energies corresponding to the N − 1 particles sector. This is particularly useful for the calculation of the GF in the form of Eq. (5). A quantum version of the EOM algorithm, the quantum equation of motion method (qEOM), was introduced by Ollitrault et al. in Ref. [23], where the ground state of the system is prepared via the VQE and the matrix elements of the EOM entering in Eq. (13) are then estimated directly on a quantum computer. This last step can be efficiently done by using for example the operator averaging technique [14] for evaluating expectation values on a state prepared on quantum hardware. The qEOM method displays two significant advantages with respect to the classical EOM. First, the matrix elements entering in the secular problem, see Eq. (13), can be efficiently measured on a quantum computer under a suitable qubit encoding, while for large systems it is generally unfeasible to evaluate them with classical methods.
Second, the qEOM can also rely, in principle, on better representations of the ground state itself. Indeed, while for the EOM one should only consider a ground state which can be easily manipulated classically, the qEOM has access to a quantum register which can host highly entangled quantum states, prepared e.g. via the VQE. Both these observations are crucial in establishing a possible quantum advantage brought by the qEOM method with respect to its classical counterpart.
In Ref. [23] the qEOM algorithm was experimentally tested for the evaluation of neutral excitations of small molecules, and was shown to be robust to noise and more accurate than quantum subspace expansion [33] in numerical tests. In this work, we propose a generalized version of the qEOM algorithm which specifically targets charged excitations. This will be used to compute the spectroscopic elements entering in the GF, using an ansatz of the form of Eq. (12), which can be expressed as where |0 is the approximate and normalized ground state of the system, obtained via the VQE in this case. The normalization factor is included to ensure that the excited states obtained via qEOM are also normalized, as it is implied by the Lehmann representation. The formula above shows that, to estimate the GF, one should first obtain an approximation to the N particle ground state via the VQE method, and then apply the qEOM algorithm on a set of elementary excitation operators. Finally, one can efficiently access the spectroscopic elements entering in the GF by computing the expectation value over the ground state as written in Eq. (16). Compared to other quantum algorithms for computing excitations of the system such as subspace search [20] or overlap based methods [19], the qEOM approach has the advantage of allowing the evaluation of the spectroscopic amplitudes with no extra quantum resources compared to the VQE. This feature is shared, e.g., with the quantum subspace expansion method, although the latter was shown to be less suited than qEOM for molecular simulations [23]. On the flip side, qEOM requires a rather large of number of measurements to reconstruct the relevant matrices. This, in the true spirit of near-term quantum solutions, shifts part of the computational burden to classical resources, and could be partially alleviated via optimized measurement protocols. In the next section, we discuss the results of the experimental and numerical implementation of our algorithm applied to the two site Fermi-Hubbard model.

IV. NUMERICAL AND HARDWARE SIMULATIONS
We now consider the Fermi-Hubbard model with N = 2 sites, in the basis {|1, ↑ , |2, ↑ , |1, ↓ , |2, ↓ }, with the on-site repulsive interaction U and the hopping t Heren i,σ =ĉ † i,σĉ i,σ is the occupation number operator and U, t > 0. In our simulations we always take U = 3 and t = 1. We start by computing the ground state of Hamiltonian in Eq. (17) via the VQE procedure. To this aim, we map the Hamiltonian from the fermionic Fock space to the Hilbert space of qubits using the Jordan-Wigner (JW) transformation [34]. After this step, the Hamiltonian can be expressed as a linear combination of 11 Pauli strings. These Pauli strings can then be divided in 3 groups where the Pauli strings belonging to the same group commute with each other. We exploit this grouping when measuring the VQE cost function by applying the standard operator averaging technique and measuring simultaneously the Pauli strings entering in the same group. This allows us to run only 3 quantum circuits to estimate the expectation value of the energy for each VQE optimization step.
The encoding of the two-site Fermi-Hubbard model requires 4 qubits. We choose the heuristic ansatz state [35] shown in Fig. 1, with 4 parameters to be optimized. Here, the A entangling blocks are of the form and can be expressed using the circuit where R(θ, φ) = R y (θ + π/2)R z (φ + π). We notice that this ansatz preserves both the particle number and the spin projection S z . To enforce specific values for these observables, the two X gates at the beginning of the circuit initialize the state with the correct spin projection (S z = 0) and number of particles N = 2. The circuit consists of 8 CNOT gates (6 entering in the A operators, 2 external). The actual number of CNOT gates then amounts to 11 on ibm cairo, due to the limited connectivity of the hardware. We show a scheme of the quantum processor in Fig. 2, highlighting the 4 qubits we use to implement our algorithm. The two A gates perform a rotation in the subspace with fixed number of particles and S z for each of the two qubit sets having the same spin (up or down). Moreover, the two CNOT gates generate enough entanglement between the two sets to represent the correct ground state for the Fermi-Hubbard dimer, while preserving particle number  and spin projection. When performing the VQE, singlequbit rotations should be attached to the ansatz circuit in Fig. 1, before measuring in the computational basis to estimate the VQE cost function at each optimization step.
We execute four types of simulations for the VQE and for the overall algorithm. We start performing statevector and qasm simulations in Qiskit [36]. We then equip the qasm simulator with a noise model (via the qiskit NoiseModel.from backend functionality) taking into account the noise levels of the actual hardware. Finally we test our algorithm on the ibm cairo quantum computer. For what concerns the VQE energy minimization procedure, we use the COBYLA [37] classical optimizer for statevector simulations, while qasm, noisy qasm and hardware simulations are performed employing the SPSA algorithm [38], which performs better in the presence of noise and allows one to estimate the energy gradient with a small amount of quantum circuits runs. As we pointed  out before, each evaluation of the cost function requires the execution of only 3 quantum circuits in our case. We perform 8192 shots for each measurement of the cost function, finding it sufficient to successfully perform the optimization routine. We then employ standard measurement error mitigation to compensate for readout errors [39]. In all simulations we use the ansatz in Fig. 1 and we iterate the VQE scheme until convergence. We plot the results of the VQE optimization procedure for differ-ent types of simulations in Fig. 3-6(a-b). We initialize the state randomly and we pick up the best VQE experiment in terms of the optimal energy reached over 10 realizations. In Fig. 3-6(a), we plot the estimated energy of the parametrized state at each optimization step. For qasm, noisy qasm and hardware simulations we plot both the estimated energy (via projective measurements) and the theoretical energy, that is the energy computed via a statevector simulation of the ansatz circuit by setting all  parameters at their actual values at the corresponding optimization step.
In Fig. 3-6(b), we plot the fidelity between the approximate ground state and the true one at each optimization step F (|ψ , |φ ) = | ψ|φ | 2 ∈ [0, 1]. The fidelity for the ground state is experimentally estimated by first extracting the coefficients of the quantum state via quantum state tomography (QST). Then the fidelity is computed analytically with respect to the exact ground state of the system. As for the energy estimate, for qasm, noisy qasm and hardware simulations we plot both the estimated fidelity (via QST) and the theoretical fidelity, obtained via the corresponding statevector simulation. A summary of the results for the VQE is shown in Table 1.  We show the results in the format measured (theoretical) where measured is the direct output from the experiment, while theoretical is the value computed with statevector using the corresponding circuit parameters. All the simulations have been performed with 8192 shots and measurement error mitigation has been considered. The estimated fidelities have been obtained via QST.
Having obtained an estimate for the ground state of the system, we then implement the qEOM algorithm to reconstruct the charged excited states of the system and to compute the GF. For the qEOM, we consider a set of 4 excitation operators to retrieve the 4 (2 particle and 2 hole) excited states of the system. This is a novel extension to the originally proposed qEOM algorithm [23]. The operators are chosen as follows: We remark that other choices lead to equivalently good results. In our case, we get a 4 × 4 qEOM generalized eigenvalue problem, targeting simultaneously the two particle and two hole excited states. In particular, our selection of excitation operators targets only the particle states with S z = −1/2 and the hole states with S z = 1/2. This allows us to compute the GF selecting only the states with a specific spin symmetry, thereby reducing the size of the GEP to be solved. The matrix elements entering Eq. (13) are then estimated via projective measurements, after performing a JW mapping.
We find the generalized eigenvalue problem to be quite sensitive to noise levels: indeed, small perturbations in the matrix elements may cause high changes in the qEOM eigenvalues that approximate the excitation energies. We alleviate this problem in two ways, first by increasing the number of measurements, hence improving the overall statistical precision, and then by measuring the excitation energies with a direct evaluation of the functional defined in Eq. (11). In the second case, we find better estimates for the excitation energies with respect to the one obtained directly from solving the qEOM generalized eigenvalue problem. We further improve the results by averaging over 13 realizations of the matrix elements estimation performed with 8192 shots, giving an effective number of total shots of approximately 10 4 . We also notice that the general problem of minimizing the effect of noise in GEPs is well known in the literature [40], and we leave open for future studies the possibility of applying classical post-processing techniques to enhance the stability of the method.
Finally, we estimate the GF in the Lehmann representation by computing the spectroscopic elements in the form of Eq. (16). Since the excitation operatorsÔ n can be expressed as a sum of 4 fermionic operators, we can decompose the whole operatorÔ nĉ † α (orÔ nÔ † n ) into a sum of Pauli strings and compute the spectroscopic ele-ments by only applying single-qubit rotations at the end of the ansatz circuit (as for the operator averaging technique applied to VQE) and then performing projective measurements.
In Fig. 3-6(c-d), we report the results for the GF computed in k space, for k = 0, π, that is considering the GF computed withĉ k,σ = (ĉ 1,σ + e ikĉ 2,σ )/ √ 2. In our case, we are actually evaluating the element G k↑,k↓ of the GF, due to definition of the excitation operators. As shown in Fig. 3-6(c-d) we find a good agreement both between noisy simulations and the exact solution, and between the noisy qasm and the hardware simulations. This suggests that in this case the noise model does not severely underestimate the effect of noise in the real device. We also confirm the intuition that having a ground state which perfectly incorporates, despite the effect of noise, some conservation laws (i.e., spin and particle number) is crucial to ensure that the resulting excited states have the expected symmetries. However, while the eigenstates of the Hamiltonian corresponding to different energies should always be orthogonal, we find a considerable overlap between some of the evaluated excited states obtained via qEOM due to the approximate nature of our method. Since the Lehmann representation assumes an orthonormal basis of excited states, this has the effect of causing satellite peaks in the GF for noisy simulations, see Fig. 6(c-d). We leave deeper analysis and possible solution of this issue for future investigations.

V. CONCLUSIONS
In this work, we propose a novel method for calculating the many-body one-particle GF compatible with near-term quantum devices. This method, in alignment with the work of Endo et al. in Ref. [28], directly evaluates the transition amplitudes of the fermionic operators between the ground state and the charged excited states of the system. This is accomplished by an original extension of the qEOM algorithm introduced in Ref. [23], which allows to compute the charged excited states of the Hamiltonian on a quantum computer, and relies on the preparation of the ground state obtained from the VQE. The GF is evaluated in the frequency domain by classically combining the transition amplitudes and the excitation energies estimated via qEOM, with the use of the Lehmann representation of the GF.
To test the validity of the present proposal, we applied this method to the two-site Fermi-Hubbard model and we performed both classical and quantum simulations of the algorithm, the latter on the ibm cairo quantum processor. We demonstrate that the proposed algorithm is able to correctly reproduce the target GF in such proof-of-principle experiments, although a rather large (compared to the usual VQE requirements) number of measurements is needed to obtain reliable estimates of the excitation energies, due to statistical and hardware noise. Although we expect this measurement overhead to become more prominent for larger systems when the standard form of our algorithm is considered, more effective strategies for observable reconstruction may be applied to meet such demand in an efficient way [41][42][43][44], together with classical post-processing techniques aimed at enhancing the overall numerical stability and, possibly, the parallelization of the qEOM procedure.
All the steps involved in our proposal rely on estimating expectation values which can be efficiently evaluated on a quantum computer. While for the VQE the relevant observable to be estimated is the energy, in the qEOM algorithm the evaluation of the excited states relies on solving a GEP whose matrix elements cannot be efficiently obtained by classical means. At the same time, the dimension of the GEP problem scales only polynomially as the system size increases, thereby preserving the predicted quantum advantage.
Since the qEOM result is sensitive to the accuracy of ground state preparation, substantial improvements over classical EOM methods are expected as soon as quantum computing methods will give access to those regimes where classical approximate ground state techniques perform poorly. Moreover, a considerable advantage of using qEOM to compute excited states compared to other quantum algorithms is that it does not require additional quantum resources (e.g., deeper circuits), with respect to the ground state VQE calculation. However, this must be balanced with an increase of the overall number of required measurements.
It is important to stress that, under the proposed scheme, the quantum computation of the GF via the Lehmann representation is efficient only as long as a polynomial number of excited states is required for a reliable reconstruction of the relevant properties as the system size increases. Moreover, in practical implementations of the qEOM method described in this work one may also face a lack of orthogonality of the resulting excited states. While we have (efficiently) included the normalization factor in our simulations, adding methods to enforce orthogonality may deliver better results in the overall computation of the GF.
Future research prospects include extending the proposed method to compute other types of GFs [31] and response functions, possibly at non-zero temperature, as proposed in Ref [28], as well as scaling up the algorithm to larger quantum systems. The latter, other than being subject to general improvements of current quantum hardware, would require applying a full set of quantum error mitigation techniques [45][46][47][48][49][50] and considering more efficient and application-specific ansätze [51]. As already mentioned above, controlling the number of required measurements with effective observable reconstruction strategies [41][42][43][44] could also significantly extend the range of applicability of the qEOM approach.