Effective calculation of the Green's function in the time domain on near-term quantum processors

We propose an improved quantum algorithm to calculate the Green's function through real-time propagation, and use it to compute the retarded Green's function for the 2-, 3- and 4-site Hubbard models. This novel protocol significantly reduces the number of controlled operations when compared to those previously suggested in literature. Such reduction is quite remarkable when considering the 2-site Hubbard model, for which we show that it is possible to obtain the exact time propagation of the $|N\pm 1\rangle$ states by exponentiating one single Pauli component of the Hamiltonian, allowing us to perform the calculations on an actual superconducting quantum processor.


I. INTRODUCTION
The many-body Green's function of an interacting quantum system is a central quantity in many stateof-the-art techniques for solid state physics and quantum chemistry. For instance, the GW method [1,2] and higher order diagrammatic expansions of the self-energy [3] use the Green's function to correct the DFT band structure of solids, while dynamical mean-field theory (DMFT) [4] tackles strongly correlated materials. Furthermore, non-equilibrium Green's functions [5] are used to study the ultra-fast dynamics of systems which are excited e.g. by a strong pulse.
The exact many-body Green's function possesses many important properties. First, it can easily be connected to observables (e.g. the density of the system n(x, t) is equal to its diagonal component), and the poles of its Fourier transform correspond to the excitation energies. Finally, the retarded Green' s function is connected to the exact linear response of the system.
The computation of the exact Green's function on a classical computer is known to become exponentially hard when increasing the size of the system under study, as it requires the knowledge of the exact ground state. It becomes thus attractive to devise algorithms for the calculation of Green's functions on quantum computers, since even though quantum algorithms make use of the many-electron wavefunction to evaluate the Green's functions, there is no actual need of measuring its components. At variance with equivalent classical approaches based on the Lehman representation, the wavefunction is first mapped to the quantum circuit, then optimized and finally used -without explicitly storing it full classical description obtained, for instance, through state tomography -to probe the Green's function. Several works presented in literature are based on the quantum phase estimation [6][7][8][9] or on the Suzuki-Trotter expansion of the propagator [10][11][12]. However, both of these * Corresponding author. francesco.libbi@epfl.ch techniques are known to be expensive in terms of circuit depth and number of controlled operations, and therefore are not suitable for near-term quantum computers. One of the first efforts to design an affordable quantum algorithm for the calculation of the Green's function goes back to Endo et al. [13]. In their work these authors suggest two alternative techniques, one in the frequency domain and the other in the time domain. The former exploits the Lehman representation and calculates the excitation energies of the system through the subspacesearch variational quantum eigensolver (SSVQE) [14], and an alternative approach has been recently presented in Ref. [15], where the SSVQE is replaced with a generalized quantum equation of motion (qEOM) technique [16] for charged excitations. The second approach of Ref. [13] is based on the real time evolution of a system through the variational quantum simulation (VQS) technique [17], which requires shallower quantum circuits compared to the Suzuki-Trotter expansion of the propagator. The algorithm of Ref. [13] requires the adoption of a variational form U which can accurately describe the time evolution of both the ground state of the system |ψ and a state of the form |P i ψ (where P i is a n-qubit Pauli operator) with the same choice of parameters. However, as shown in Sec. III, it is quite difficult to find an empiric variational form satisfying this requirement. In addition, such variational form is expected to require deeper circuits than the one needed to evolve one single quantum state. In this work, we propose a novel algorithm that make use of a variational form U evolving only the state |P i ψ , and show that it leads to much shallower circuits. The paper is organized as follows. In Sec. II we introduce the definition of the components of the Green's function, and show how these can be calculated through real-time evolution techniques. In Sec. III we discuss the choice of the ansatz for the VQS algorithm, listing the most important requirements and suggesting possible strategies to comply with these. Eventually we will see how the choice of a variational form falls naturally upon the variational Hamiltonian ansatz [18,19]. In Sec. IV we present our novel algorithm for the calculation of the Green's function, and in Sec. V we compare the results obtained from numerical simulations with those from the algorithm of Ref. [13] for the case study of a 2-site Hubbard model. We show that for this particular problem the exact evolution of the state P i |ψ can be obtained by exponentiating only one of the Pauli components of the qubit Hamiltonian, while the same simplification does not hold when evolving |ψ and |P i ψ at the same time. The quantum circuit corresponding to the algorithm suggested is shallow enough to be executed on current quantum computers. In Sec. VI we compare the results obtained with the two algorithms for the 3-site Hubbard model with open-boundary conditions and 4-site Hubbard model with periodic-boundary conditions, and discuss the scaling for larger problems. Finally, in Sec. VII we show the results for the 2-site Hubbard model obtained with state-of-the-art superconducting IBM Quantum processors.

II. GREEN'S FUNCTION THROUGH REAL-TIME EVOLUTION
Before presenting the algorithm for the calculation of the correlations, it is useful to define the components of the Green's function of relevance for this work. The lesser Green's function is defined as where |ψ is the ground state of the Hamiltonian and c † m , c l are fermionic creation and annihilation operators respectively. The greater component is instead The difference between these components is the retarded Green's function, which is the target of our calculations for its applications in linear-response theory: where θ(t) is the Heavyside step function. The real-time approach to the calculation of the Green's function consists in computing the states e −iHt |ψ and e −iHt c † m |ψ , which represent the time evolution of the states |ψ and c † m |ψ respectively. The choice of the time evolution technique, which is crucial for designing an efficient quantum circuit for the calculation of the Green's function, is discussed in the next paragraph.

III. CHOICE OF THE VARIATIONAL FORM FOR THE VQS
We briefly introduce the techniques suitable to follow the time evolution of a quantum system, and then focus on the choice of a proper ansatz for the VQS. The time evolution of a generic quantum state |φ 0 is obtained through the action of the time evolution operator: On a quantum computer, this can be achieved by applying the unitary gate expressing the exponential of the Hamiltonian to the qubit register encoding the state |φ 0 . The qubit Hamiltonian, which is built thanks to the mapping of the fermionic Hamiltonian to the qubit Hilbert space, is a weighted sum of multiple-qubits Pauli gates Therefore, the exact expression of the propagator is Since the operators P m does not necessarily commute with each other, it is not possible to express U as a product of exponentials in an exact way. Nevertheless, thanks to the Suzuki-Trotter expansion, the propagator can be approximated as a product of exponentials: Unfortunately, this approach requires deep circuits, especially when the time evolution is carried for long times, and therefore it is not optimal for near-term noisy quantum devices. A remarkable advantage in term of circuit depth can instead be obtained by adopting variational quantum simulation (VQS) algorithms [17,20]. These techniques consist in optimizing the parameters of a variational form U (θ = θ 1 , ..., θ M ) to satisfy the variational principles which hold for the exact time evolution. The variational principle most often adopted for quantum simulations is McLachlan's, due to its numerical stability, and reads δ||d/dt + iH |φ(θ(t)) || = 0 .
Eq. 8 leads to the following differential equation for the parametersθ: where and Here, we have adopted the shortcut notations |φ = |φ(θ(t)) and |∂ i φ = ∂ ∂θi |φ(θ(t)) . The McLachlan variational principle has a simple and instructive geometric interpretation, which is discussed in Ref. [21]: the variational form U (θ 1 , ...θ M ) represents a M-dimensional manifold embedded in the N dimension Hilbert space (N being equal to 2 n , with n number of qubits). In most cases the parameters θ i are real, and therefore it is more convenient to describe the Hilbert space as a 2N -dimensional real vector space, i.e. a space which is spanned by 2N complex vectors with real coefficients (see Ref. [21] for a detailed discussion). In a real vector space, the tangent space to the manifold M associated to the variational form U (θ) is spanned by the vectors: It is easy to check that T i |φ = 0, ∀i = 1, ...M , as |φ belongs to the subspace of vectors with norm 1. In Appendix A, we show that Eq. 9 can be elegantly written as with T i |T j being the Quantum Geometric Tensor [22,23]. Eq. 13 states that only the component of the vector H |φ which is parallel to the tangent subspace to the manifold at the point |φ drives the time evolution. Note that when H |φ is perpendicular to the tangent subspace at the manifold, the right hand side of Eq. 13 vanishes, so there is no time evolution (apart form a global phase). This happens correctly when |φ is eigenstate of the Hamiltonian, and T i |H|φ = E T i |φ = 0. However, this can happen also due to a poor choice of the variational form, leading to an unphysical evolution where the state remains constant. Unfortunately, the choice of a good variational form is in general quite problematic, as we will show in the following. A variational form is actually a circuit expressing the unitary operator U (θ) which is applied to the state |φ 0 to be evolved. As such, it must satisfy the following conditions. (I) It must be equal to the identity for the initial choice of the parametersθ 0 . It is not straightforward to identify the parameters which turn the variational form into the identity, and most often it is easier to design variational forms that are equal to the identity forθ 0 = 0. (II) The number of parameters must be much smaller than the size of the Hilbert space, otherwise the advantage of the variational algorithm would be lost. (III) It must be such that the the projection of the vector H |φ (which is proportional to the exact time derivative of the state |φ ) on the tangent subspace is as large as possible, as underscored by Eq. 13. (IV) In order to be realizable on near-term quantum computers, it must contain a limited number of controlled operations, and, more in general, lead to sufficiently shallow circuits. We can now distinguish two main kinds of variational forms. Empiric ansätze are expressive quantum circuits usually characterized by the alternation of CNOTs which create entaglement and rotations for the parametrization [24]. As opposite to empiric ansätze, the physically inspired ansätze are based on the physical properties of the system. The most important representative of this category is the variational Hamiltonian ansatz (VHA) [18,19]. This takes inspiration from the Suzuki-Trotter expansion of Eq. 7, and the parametrization is introduced by replacing the product c m t with a function θ m (t): here, n d is the depth of the circuit, and corresponds to the number of Trotter steps. We will now highlight the issues that can drive from the choice of an empiric ansatz, and show how those can be circumvented by adopting the VHA.
The initial value problem. As stated before, any ansatz must contain entangling operators which express an interacting many-body state. One of the most used entangling operator is the CNOT gate, since it is native on many superconducting devices. The position of these CNOTs must be carefully engineered in such a way that the ansatz is equal to the identity at the beginning of the evolution, as required by point I. This is not an easy task to achieve while keeping the circuit shallow at the same time. A strategy to do so can exploit the identity CNOT 2 = I and compose a circuit of units such as those represented in Fig. 1, where a layer of CNOTs is followed by a layer of rotations which are equal to the identity for θ = 0 and then another layer of CNOTs with the order inverted. We will see that this kind of structure looks very similar to the circuit to exponentiate the single Pauli operator P m of the qubit Hamiltonian. Another possibility is to use an ansatz of the following form: which guarantees by construction condition I. However, both strategies require to double the number of CNOTs, and are thus less suited for proof-of-principle demonstrations on near-term quantum processors (point IV).
The problem of the initialization of the time evolution. As anticipated above, the VQS can predict erroneously the initial state to remain constant in the time, if the vector H |φ 0 is orthogonal to the tangent subspace at the point |φ 0 of the manifold associated to the variational form. Unfortunately, this can happen often, and represents a major issue in the choice of an empiric ansatz. The reason is simple: the vector H |φ 0 points to one of the 2N possible directions of the real vector space associated to the Hilbert space. Choosing an empiric ansatz of M parameters, with M considerably smaller than N (in order to satisfy point II), consists in drawing a Mdimensional surface in a much larger space, while hoping that this surface be parallel to a specific direction, which is likely to fail. Things are made more difficult by the fact that different VQS parameters may span the same direction in the tangent subspace. An example of this is represented by the ansätze with tunable depth which are obtained by repeating the same blockŨ a number of times equal to the circuit depth n d . It can be easily proven that increasing the circuit depth does not increase the size of the tangent subspace: as expressed in Eq. 12, the tangent subspace is generated by the derivatives of the variational form with respect to its parameters; if we consider two corresponding parameters θ d1 i and θ d2 i belonging to the layers d 1 and d 2 , we can prove that the derivative of the ansatz U with respect to these two parameters are equal to each other at the origin (in the case in which this coincides with the pointθ = 0): in the case in whichŨ (0) = I and the blocksŨ are equal to each other. The orthogonality problem described above becomes more and more important when the number of qubits (i.e. the size of the Hilbert space) grows. Now, we know that |Hφ = m P m |φ , and we could try to engineer the empiric ansatz in a way such that the product of the vectors |T i and |P m φ is not vanishing. The most straightforward way to do this is to impose that ∂ m |φ = iP m |φ , but this is exactly the VHA ansatz introduced above. In fact thee VHA ansatz does not suffer from the problems analysed. First of all, it is equal by construction to the identity whenθ = 0. Furthermore, it can be proven (see Appendix B) that when adopting this variational form, the VQS solution coincides with the exact propagator for vanishing small times, with parameters and thus respects exactly the condition III. This happens for any state |φ 0 that is evolved. In many cases In Ref. [10][11][12] the unitary operator evolving the system is calculated through Suzuki-Trotter expansion.
Ref. [13], instead, express the propagator e −iHt through a variational form U (θ) whose parameters are determined through the VQS McLachlan algorithm. It is important to stress that such approach works only if U (θ) approximates accurately the time evolution of the states |ψ and P j |ψ ; otherwise it becomes necessary to adopt the quantum circuit shown in panel (b), which requires a controlled U (θ).
of interest, the combination of this ansatz with the VQS algorithm allows to reduce remarkably the circuit depth with respect to the Suzuki-Trotter expansion technique. The fact that the VHA is a natural choice for the time evolution of states especially when the size of the Hilbert space grows must be taken into account when designing algorithms for near-term quantum computers and in the following we leverage on its properties to propose a new and optimized algorithm for the calculation of Green's functions.

IV. OPTIMIZED ALGORITHM
The standard technique to calculate correlation functions [10][11][12] adopts the circuit shown in Fig 2a, where the propagator e −iHt is calculated using the Suzuki-Trotter expansion. Endo et al. [13] improved this algorithm by replacing the Trotter expansion, which requires very deep circuits, with a variational form U whose parameters are determined through the VQS algorithm. In order for this to work, the operator U has to accurately approximate the time evolution for the ground state |ψ and the state |P i ψ at the same time, i.e. with the same parameters θ i (the evolution of the ground state consists simply in adding the correct phase without modifying the state). This requirement allows one to eliminate the controlled variational form CU from the circuit (see Fig.  2b), which would be otherwise necessary when U evolves |ψ and |P i ψ with a different choice of the parameters for the two states. For this reason, we shall refer to this algorithm as control-free (CF). Requiring U to evolve two states at the same time seems to be a good price to pay to remove the CU operations, since a controlled unitary operation acting on n qubits is in general very complicated to realize and can require a large number of CNOTs for its implementation. However, it turns out that this strategy is not always advantageous: first of all, as discussed in the previous section, it is difficult to find a shallow empiric variational form that can propagate a desired state with a sufficiently small error. This task is even harder as the size of the Hilbert space grows, particularly if one requires the ansatz to evolve two states simultaneously. This leads naturally to the choice of the VHA ansatz. Second, we argue that it is possible to control the VHA by adding a minimal number of controlled operations, as shown in Fig. 3b. Therefore, the use of c-VHA does not necessarily imply a significant increase of controlled operations and, consequently, unaffordable circuit depths. Finally, evolving two states at the same time requires a deeper variational ansatz, with a larger cost in terms of controlled operations and circuit depth, as it will be proven numerically in the following sections. This is because the variational ansatz must approximate the ideal time evolution operator in at least two independent directions in the Hilbert space. It is easy to conclude that an algorithm for the Green's function which requires a controlled evolution operator evolving one single state at a time can lead to an improvement in term of circuit depth and number of CNOTs. This improved algorithm can be obtained under the assumption of a time independent Hamiltonian, which holds in many problems studied in the literature. In this case, the propagator is e −iHt , and its action on the ground state introduces simply a phase rotation e −iHt |ψ = e −iE0t |ψ , E 0 being the ground state energy which can be calculated through the variational quantum eigensolver (VQE) algorithm. Using this relation, the expression of the lesser Green's function becomes After mapping the fermionic Hamiltonian on a qubit register, the expression for the Green's function becomes where U is the ansatz of the VQS algorithm, and λ i are the coefficients of the expansion of the creation/destruction operators in term of Pauli gates. Since the operator P i U P † j is not Hermitian, the braket must be computed through the Hadamard test; the circuit is the one shown in Fig. 2. It can be easily observed that in this case the ansatz only implies the evolution of the state P i |ψ at the price of a controlled evolution operator. Furthermore, when P i = P j , it is possible to remove the control operation on these gates, as shown in Fig. 4b. Since the operator U † is removed from the expression,  [25]. Panel (a) shows an example of the circuit corresponding to the exponential of the four-qubit Pauli gate ZXZY (namely e −iθZXZY ). Panel (b) shows its controlled form. Notice that the controlled form is realized by simply controlling the central Rz gate alone. As a result, the controlled VHA simply requires the addition of a control operation for each appearing Pauli term, which is minimal considering that the exponential of a n-qubit Pauli gates requires at most 2(n − 1) CNOTS.
it is necessary to take correctly into account the global phase of the state U |P j ψ . This can be done in a very simple way, as suggested in Ref. [17], without the need to add a global phase to the ansatz used for the time evolution: indeed, the time dependent global phase θ 0 (t) ensuring that the state e iθ0(t) U (θ 1 (t), ..., θ m (t)) |φ 0 coincides with the time evolution of the state |φ 0 can be determined using the McLachlan principle. By manipulating Eq. 9 we obtain the following differential equation: where theθ i for i = 1, . . . , m are those calculated through the VQS method. In order to obtain the correct brakets, it is necessary to multiply the result of the circuit in Fig.  4 by e iθ0(t) .
The very same procedure can be followed to calculate the greater component of G, and thus the retarded correlation function. Since the proposed algorithm requires the evolution of one state at time, we shall refer to this algorithm as one-state (OS). In the next sections we will show the advantage of this algorithm for the 2-and 3-  . 4: Panel (a) shows the circuit to calculate the braket of Eq. 20. When φ = 0, the Z-measurement of the first qubit gives the real part of the braket, while the imaginary part can be obtained by setting φ = − π 2 . The controlled operations on the Pauli gates P i and P j can be further removed if P i = P j , as shown in panel (b).

V. 2-SITE HUBBARD MODEL
First we test our algorithm on the 2-site Hubbard model. It is possible to take advantage of the many symmetries of this model in order to reduce the depth of the ansatz used to evolve a generic state P i |ψ . In particular, here we prove that the exact evolution of the state P i |ψ is obtained by simply exponentiating one single Pauli component of the Hamiltonian (i.e. one among XZXI, YZYI, IXZX, IYZY). It is important to note that the present algorithm can exploit these peculiar symmetry properties, thanks to the fact that it only requires to implement the time-evolution of a single state P i |ψ . Conversely, the CF algorithm, where the VHA ansatz has to evolve the states P i |ψ and |ψ with the same parameters, requires not only the exponentiation of all Pauli terms in the Hamiltonian, but also the doubling of the circuit depth (depth=2), resulting in much deeper quantum circuits (see Fig. 5). We start introducing the Hubbard dimer Hamiltonian (23) Here we use the mapping (1 ↑, 1 ↓, 2 ↑, 2 ↓) → (q 1 , q 2 , q 3 , q 4 ). It is now possible to recognize two sets of operators, (iii) P 2 l = I .
These symmetries have very important consequences on the choices of the ansatz for the time evolution. We will derive now a series of properties which will lead to drastic simplification of the time evolution circuit for the propagation of the Green's functions.
Proposition 1. If the ground state |ψ of the Hamiltonian is non-degenerate, then ∀ P l , P m ∈ S i =⇒ P l P m |ψ = |ψ .
Proof: By definition, the ground state is the vector |ψ such that ψ|H|ψ is minimum. However, ψ|H|ψ = ψ|P † m P † l HP l P m |ψ from property (ii) above. Since we assume that the ground state is non-degenerate, then P l P m |ψ = α|ψ . Using the commutation of P l and P m (i) we have also that P m P l |ψ = α|ψ . Now, P l P m P m P l |ψ = α 2 |ψ . At the same time, since Pauli matrices square to unity, P l P m P m P l |ψ = |ψ ; it follows that α = ±1. Finally, ψ|P l |ψ + ψ|P m |ψ = ψ|P m P l P m |ψ + ψ|P l P m P l |ψ = α( ψ|P l |ψ + ψ|P m |ψ ); therefore, α = 1.
Following Proposition 1, it is possible to derive the following corollary. Proposition 2. If the ground state |ψ of the Hamiltonian is non-degenerate, then ∀ P l , P m ∈ S i =⇒ P l |ψ = P m |ψ . This result can be obtained by multiplying both the sides of the result of Proposition 1 by P l and using the property (iii) above.
Let us now proceed with the analysis of the time evolution of a generic state, which is obtained by applying the creation/destruction operators to the ground state |ψ .
Because of the Jordan-Wigner transformation, this corresponds to applying a Pauli operatorP (e.g. ZZXI) to the ground state; a tilde is used to remark thatP / ∈ S i . Let us consider P l , P m ∈ S i . For Proposition 2, we have that P l |ψ = P m |ψ = |χ . Therefore, P lP |ψ =P |χ + [P l ,P ] |ψ and P mP |ψ =P |χ + [P m ,P ] |ψ .
Using the simple rules described in Appendix C, we can then build the following table showing the sign of the term P lP |ψ (i.e. ±P |χ ) for each couple of P l (rows) andP (columns). It is now possible to relate these results XIII ZXII ZZXI ZZZX coeff. XZXI + + -+ -τ /2 YZYI - with the time evolution of the operatorP |ψ . At first order in the expansion in t, As the coefficients are the same for all the operators belonging to S i , then the contributions of the elements of S 2 always cancel each other, while that of the elements of S 1 is always 2P |χ , so where P 1 is one of the Pauli matrices belonging to S 1 . For the second order expansion in t, The second term in the r.h.s. of Eq. 27 vanishes as shown above, while the first term is equal to 4c 2 1 P 2 1 . It is straightforward to generalize the result at any order in perturbation theory, and conclude that the exact evolution ofP |ψ is obtained by applying the operator U = e iστ P l t (28) where P l ∈ S 1 and the sign σ can be read in Tab. I for each specific operatorP .
We make now use of this compact form of the propagator in the simulation of the two-site Hubbard model. Fig. 5 shows the results for the retarded Green's function calculated numerically with the Qiskit qasm simulator [26]. The Green's function obtained with the OS algorithm is represented by green squares. Thanks to the properties shown above, the ansatz adopted is the one of Eq. 28, which contains only 5 CNOTs (in its controlled form). The Green's function obtained with the CF algorithm is represented with light blue circles (d=1) and violet diamonds (d=2). The former contains 20 CNOTs and fails to reproduce the exact Green's function dynamics, while the latter is in good agreement with it, but requires 40 CNOTs, as shown in Tab. II. In summary, the novel OS algorithm manages to reach the same accuracy as the CF does; however, with a number of CNOTs and a total gate count of about one ninth.

VI. SCALING TO LARGER PROBLEM SIZES
In order to provide further evidence about the advantages of the OS algorithm, we also performed calculations of the Green's function for the 3-site Hubbard chain with  Table II: Number of 1-qubit and 2-qubit gates in the ansätze adopted to calculate the Green's function for the Hubbard model. The number reported for the OS algorithm corresponds to the controlled form of the ansatz, which is the one adopted to calculate the Green's function.
open-boundary conditions and the 4-site Hubbard chain with periodic-boundary conditions. In these cases, it is not possible to exploit the symmetries of the Hamiltonian to simplify the ansatz as much as for the 2-site Hubbard model. Nonetheless, we show that the suggested algorithm allows us to obtain an accurate Green's function with shallow circuits.
The 3-site Hubbard model is characterised by 2 degenerates ground states; here, we calculate the correlation function for one of these two ground states. The 4-site Hubbard model has instead a non-degenerate ground state. The retarded Green's function for both the models is computed at a k value such that c k = c 1↑ − c 2↑ . The result is shown in Fig. 6: we observe that the Green's function for the 3-site Hubbard model obtained with the new algorithm using a VHA of depth d = 3 shows the same accuracy as the one computed with the CF algorithm with a depth d = 5, which fails to reproduce the exact Green's function when the depth is reduced to d = 3. The same happens for the 4-site Hubbard model. The advantage in term of number of CNOTs and circuit depth is shown in Tab. II.
To further investigate the scaling properties of the OS algorithm, let us summarize the main differences with respect to the CF protocol. On one hand, our strategy generally requires smaller VHA depths, thus reducing significantly the number of quantum operations necessary for a faithful description of the dynamics. On the other hand, one additional ancilla-controlled operation is necessary per Pauli string appearing in the Hamiltonian. Under the decomposition scheme presented in Fig. 3a employing the native gate set of current IBM Quantum processors -the exponentiation of a single Pauli component of the Hamiltonian requires 2(w − 1) CNOTs, where w is the corresponding Pauli weight, i.e. the number of non-identity elements in the Pauli string. The total number of controlled operations required by the OS algorithm for a VHA of depth d 1 is thus: where n p is the number of Pauli components of the Hamiltonian, w i is the weight of the i-th component andw is the average Pauli weight, defined asw = 1/n p np i=1 w i . The CF algorithm requires instead, for a VHA of depth d 2 Combining Eqs. 29 and 30 it is easy to see that the suggested algorithm will show an advantage over the CF one when .
Clearly, the OS algorithm is favoured for larger average Pauli weights. In this respect, the n-site Hubbard chain with open-boundary conditions and first nearest neighbours coupling represents the least favourable case: indeed, the average Pauli weightw = (14n − 12)/(5n − 4) has an asymptotic value of 2.8, thus requiring d 2 /d 1 > 1.27 for the OS algorithm to be competitive. However, when replacing the open-boundary conditions with periodic-boundary conditions, or when assuming that the Hubbard sites are disposed on a 2d lattice instead of a chain, the appearance of Pauli components of maximum Pauli weight makes the condition of Eq. 31 much more favourable (see Appendix D). Finally, in a typical quantum chemistry problem, the interacting term h pqrs a † p a † q a r a s is mapped into It is worth reminding that, while all the examples above were treated by making use of the Jordan-Wigner fermion-to-qubit mapping, other strategies are known, such as the Bravyi-Kitaev mapping [27], which in general leads to lower average Pauli weight in the asymptotic limit. These alternative mappings also have an impact on the VHA implementation details and required depths, thus requiring a case-by-case assessment [28].
Finally, it is also important to remark that the controlled VHA requires in general a higher connectivity between the ancilla and the system qubits. As a result, native hardware coupling maps and circuit compilation protocols must also be taken into account when selecting the the most favorable strategy for specific applications.

VII. HARDWARE RESULTS
The algorithm presented in this paper allows one to reduce remarkably the size of the circuit used to compute Green's functions for several model Hamiltonians. In the case of the 2-site Hubbard model, the improvement is such that it allows the execution of the algorithm on a current quantum processor. More specifically, we demonstrate an implementation of our OS algorithm on IBM Quantum superconducting devices, accessible via the cloud [29].
We first obtain an approximate ground state using the VQE algorithm with the SPSA optimizer on ibmq manila. The adopted variational form is represented in Fig. 7 (within the red box), and consists of a concatenation of single R y qubit rotations and CNOTs.
For this experiment, we used a Hubbard Hamiltonian with t = 1 and U = 3. The theoretical ground state energy corresponds to E 0 = −4, while the experimental energy obtained with the VQE is -3.63. This result can be further improved on ibmq montreal, due to the lower noise level of the processor. By also applying readout error mitigation [30], we obtain a ground state energy of −3.91.
After the optimization, we calculate the diagonal component of the retarded Green's function, G R 1↑1↑ . The full circuit used for the calculation of the expectation values in Eq. 20 is given in Fig 7 for the case of P i = P j = XIII at t = 0; further details on the calculation setup are given in Appendix E. Due to the particle-hole symmetry, its real part is vanishing and only the imaginary part survives. The calculation is performed on ibmq montreal, and the results are shown in Fig. 8a. This figure, however is not meaningful for a visual comparison with the experiments, due to the frequency shift introduced by the error on E 0 . It is instead much more instructive to consider the Green's function obtained from Eq. 20, where all matrix elements are measured on the quantum processor while the value of the energy in the phase of e iE0t is taken to be equal to its theoretical value, E 0 = −4 (instead of taking the approximated variational energy computed on ibmq montreal ) in Fig. 8c. The reasons are threefold. First, the adoption of the energy E 0 eliminates the small frequency shift and allows to appreciate visually the accuracy of the algorithm for the evaluation of the brakets in Eq. 20. Second, the ground state energy is just a parameter in Eq. 20 that could be evaluated with higher accuracy using other techniques (such as imaginary time evolution). Finally and most importantly, a small error FIG. 7: Circuit required to calculate the braket in Eq. 20 when P i = P j = XIII at t = 0. The part of the circuit which is enclosed in the red rectangle corresponds to the variational form adopted for the VQE.
in the ground state energy only leads to a small energy shift of the poles of the Green's function in frequency domain, as it can be seen by comparing Fig. 8b with Fig.  8d, where the noisy energyẼ 0 and the exact energy E 0 are used respectively (the proof is given in Appendix E). This agrees with the general observation that, thorugh e.g. Fourier analysis, useful and accurate information in the frequency domain can be obtained even from noisy quantum data [12,31]. Overall, the agreement shown in Fig. 8c between the imaginary part of the exact Green's function and the one calculated on the hardware is very good, also considering that no further error mitigation scheme is adopted at this stage. The real part of the Green's function obtained on hardware shows small oscillations around the exact value 0. These oscillations are mainly caused by the errors of the hardware in reproducing the true ground state wavefunction. In the frequency domain, the agreement between the exact result and the quantum computation is overall very good, both for the shape and the position of the peaks, also when the noisy energyẼ 0 is adopted.

VIII. CONCLUSIONS
In this work, we proposed an optimized real time propagation algorithm for the calculation of the Green's function on a near-term quantum computer, which exploits the time independence of the Hamiltonian and the features of the VHA ansatz and requires the time evolution of the only state P i |ψ . We proved that this algorithm is more efficient in term of 2-qubit gates than the previously known ones, both for the 2-site 3-site and 4-site Hubbard models, and we estimates that the reported advantage is likely to persist when considering complex realistic quantum chemistry problems. In particular, this algorithm is able exploit the symmetries of the 2-site Hubbard model, where the exact propagator of the state P i |ψ is obtained by exponentiating only one of the Pauli components of the Hamiltonian. This leads to a very shallow circuit for the calculation of the Green's function which, can then be executed on a current quantum computer, with results in remarkable agreement with the exact solutions. We studied the scaling of the algorithm with the size of the problems, and proved that the algorithm proposed is favoured in problems where the Pauli weight of the components of the Hamiltonian is large. This is the case, for example, of quantum chemistry problems with interactions beyond first nearest neighbours. We believe that the present work will promote further research towards the development of hardware efficient quantum algorithms for calculation of electronic correlation functions, opening up new avenues in the study of complex many-body systems and their excitations on near-term quantum computers.  Similarly, since φ|H|φ is real, then i ∂ i φ|φ φ|H|φ = − ∂ i φ|φ φ|H|φ .

IX. ACKNOWLEDGMENTS
Substituting in Eq. 11 of the main text, In this section we show that, when using the VHA anzatz (Eq. 14), McLachlan-VQS solution coincide with the exact evolution for vanishingly small times. At t=0, the partial derivative of the VHA with respect to the parameter θ m is Let us consider two n-dimensional Pauli operators, P, P . We can write (without restriction) P = P 1 P 2 ...P m P m+1 ...P n P = P 1 P 2 ...P m P m+1 ...P n where the first corresponding m operators of P and P are Pauli operators which are different from each other and from identity (e.g P 1 = X, P 1 = Y ) and the corresponding operators from m + 1 to n are either equal to each other or one of them is equal to the identity (e.g P n = I, P n = Y or P n = Y , P n = Y ). We have that P P = (iε 11 1 P 1 )...(iε mm m P m )(P m+1 P m+1 ...P n P n ) (C1) = (iε 11 1 P 1 )...(iε mm m P m )(P m+1 P m+1 ...P n P n ) −(−iε 11 1 P 1 )...(−iε mm m P m )(P m+1 P m+1 ...P n P n ) with g < (t) = −ie iE0t ψ|c k e −iHt c † k |ψ periodic on T = 4π. We therefore calculate g(t) over this time range and then its Fourier transform g(ω). Thanks to the convolution theorem, therefore we just need to shift the function g(ω) of E 0 −Ẽ 0 in order to obtain the correct Fourier transform. We can repeat the very same procedure for the greater term and thus obtain the retarded component.