Calculation of the Green's function on near-term quantum computers
Suguru Endo Iori Kurata Yuya O. Nakagawa



 endo@qunasys.com
 QunaSys Inc., Aqua Hakusan Building 9F, 1-13-7 Hakusan, Bunkyo, Tokyo 113-0001, Japan
 

 iori-kurata636@g.ecc.u-tokyo.ac.jp
 Department of Applied Physics, The University of Tokyo, Tokyo 113-8656, Japan
 

 nakagawa@qunasys.com
 QunaSys Inc., Aqua Hakusan Building 9F, 1-13-7 Hakusan, Bunkyo, Tokyo 113-0001, Japan
 




The Green's function plays a crucial role when studying the nature of quantum many-body systems, especially strongly-correlated systems.
Although the development of quantum computers in the near future may enable us to compute energy spectra of classically-intractable systems, methods to simulate the Green's function with near-term quantum algorithms have not been proposed yet.
Here, we propose two methods to calculate the Green's function of a given Hamiltonian on near-term quantum computers.
The first one makes use of a variational dynamics simulation of quantum systems and computes the dynamics of the Green's function in real time directly.
The second one utilizes the Lehmann representation of the Green's function and a method which calculates excited states of the Hamiltonian.
Both methods require shallow quantum circuits and are  compatible with near-term quantum computers. 
We numerically simulated the Green's function of the Fermi-Hubbard model and demonstrated the validity of our proposals.



Calculation of the Green's function on near-term quantum computers
    Yuya O. Nakagawa
    December 30, 2023
==================================================================






§ INTRODUCTION

The advent of a primitive but still powerful form of quantum computers, called noisy intermediate-scale quantum (NISQ) devices, is approaching <cit.>.
NISQ devices have a few hundreds to thousands of qubits under highly precise control but they are not fault-tolerant.
It is believed that the behaviour of NISQ devices will soon reach a stage that one cannot simulate its dynamics using classical computers, due to exponentially-increasing size of the Hilbert space with the number of qubits <cit.>.
Therefore, many researchers expect that NISQ devices will exhibit supremacy over classical computers for some specific tasks, even though they cannot execute complicated quantum algorithms requiring a huge number of qubits and gate operations due to the erroneous nature of them <cit.>.

One of the most promising tasks in which near-term quantum computers may outperform classical computers is quantum simulation, where one computes energy eigenvalues and/or eigenstates of a given quantum system.
It enables us to calculate and predict properties of quantum many-body systems, which is of great importance to many fields such as quantum chemistry, condensed matter physics, and material science <cit.>.
The most celebrated algorithm for quantum simulation on near-term quantum computers is the variational quantum eigensolver (VQE) <cit.>, in which energy eigenstates and eigenenergies of the system are obtained based on the variational principle of quantum mechanics.
Although the VQE was originally proposed for finding only the ground state, its extension to excited states was discussed in several recent works <cit.>.

However, other important quantities to investigate quantum many-body systems other than eigenenergies and eigenstates have remained relatively disregarded in the recent development of near-term quantum algorithms, i.e., the Green's function and the spectral function <cit.>.
They are fundamental to study quantum many-body systems, especially strongly-correlated systems;
for example, in condensed matter physics, the spectral function tells us that the dispersion relation of quasipaticle excitations of a system, which gives crucial information on  high-T_c superconductivity <cit.>, magnetic materials <cit.>, and topological insulators <cit.>. 
While several methods based on the Suzuki-Trotter decomposition of the time evolution operator <cit.> or quantum phase estimation <cit.> were already proposed to calculate the Green's function on quantum computers (including general multi-point time correlation functions <cit.> used in gauge theories and nuclear physics <cit.>),  they require a large number of qubits and gate operations which are hard to realize with near-term quantum computers.

In this paper, we propose two different algorithms for evaluating the Green's function on near-term quantum computers.
The first one takes advantage of the variational quantum simulation (VQS) algorithm <cit.> for an efficient calculation of the Green's function in real time.
We extend the original VQS algorithm in order to calculate the transition amplitude of general quantum operators between two different quantum states after time evolution.
The second one is based on the Lehmann representation of the spectral function <cit.>.
By calculating excited states of a given system and evaluating the transition amplitude of appropriate operators by using the subspace-search variational quantum eigensolver (SSVQE) <cit.> or the multistate contracted VQE (MCVQE) <cit.>, one can compute the spectral function (hence the Green's function).
We confirm the validity of our methods using numerical simulations of the Fermi-Hubbard model, a model of strongly-correlated system.
The extension of our methods to finite temperature and general correlation functions is also discussed.

The rest of the paper is organized as follows.
In Sec. <ref>, we briefly review the definition of the Green's function and the spectral function at zero temperature.
In Sec. <ref>, we propose a method to calculate them by using the VQS algorithm.
We describe another method to calculate the Green's function and the spectral function as a simple application of the SSVQE algorithm <cit.> and the MCVQE algorithm <cit.> in Sec. <ref>.
Sec. <ref> is dedicated to demonstration of our methods by performing numerical simulations calculating the spectral functions of the 2-site Fermi-Hubbard model. Sec. <ref> discusses how the depth of the ansatz for the VQS affects the accuracy of numerical simulations.
In Sec. <ref>, we discuss the feasibility of implementing our proposed algorithms on near-term quantum computers from the viewpoint of the number of gate operations and their error rate required. 
We present the extension of our methods to the finite temperature Green's function in Sec. <ref>.
We discuss implications and possible future directions of our results, and make a conclusion in Sec. <ref>.
Appendix <ref> provides detailed resource estimations of our algorithms. In Appendix <ref>, we describe details of the discussion in Sec. <ref>.  Appendix <ref> provides additional numerical results, e.g., four-site Hubbard model simulation.
Appendix <ref> gives details of the numerical simulations.



§ REVIEW OF THE GREEN'S FUNCTION AND THE SPECTRAL FUNCTION AT ZERO TEMPERATURE 

In this section, we briefly review definition of the Green's function at zero temperature for consistency <cit.>.

Let us consider a fermionic system described by Hamiltonian H which is composed of fermionic creation (annihilation) operators c_a, c_a^†, where a is a label to specify the fermionic mode.
For example, a can be (k,σ), where k denotes the momentum and σ=↑, ↓ denotes the spin of the fermion.
The retarded Green's function at zero temperature is defined as

    G^R_ab(t)  = -i Θ(t) ⟨c_a(t) c_b^†(0) + c_b^†(0)c_a(t) |_⟩0,

where Θ(t) is the Heaviside step function, c_a(t) = e^iHt c_a e^-iHt is the Heisenberg representation of the operator c_a, and ⟨…|_⟩0 = ⟨G|…|G|$⟩ denotes the expectation value by the ground state of the Hamiltonian,|G⟩.
We employ the natural unit where the Plank constantħand the Boltzmann constantk_Bareħ=k_B=1.
For simplicity, throughout this paper we consider the Green's functionG^R_abwitha=b=(k, ↑), namely, the Green's function in the momentum space with identical spin.
We simply write

    G^R_(k,↑),(k,↑)(t)  = G^R_k(t)

in all other parts of the paper.
We note that extension of proposed methods in this study to the Green's function with general indices is straightforward.

The Green's function is related to another important physical quantity to investigate quantum many-body systems, namely, the spectral functionA_k(ω).
It is defined through the Fourier transform ofG^R_k(t),

    G^R_k(ω)  = ∫_-∞^∞ dt   e^i(ω+iη)t G^R_k(t) 
    
    =: ∫_-∞^∞ dω' A_k(ω')/ω - ω' + iη,

whereη→+0is a factor to assure the convergence of the integral.
The spectral function and the Green's functionG_k^R(ω)have a relation

    A_k(ω) = - π^-1Im G_k^R(ω).

Finally, we introduce the Lehmann representation of the spectral function which utilizes the energy eigenvalueE_nand the eigenstate of the Hamiltonian|E_n⟩,

    A_k(ω)    =   ∑_n(  |⟨E_n|c_k^†|G||⟩^2 /ω+E_G-E_n+iη. 
          . +  |⟨E_n|c_k|G||⟩^2 /ω-E_G+E_n+iη),

where we call the first (second) term as particle (hole) part of the spectrum function.

In Sec. <ref>, we compute the spectral function by performing the Fourier transformation toG_k^R(t)calculated by the VQS-based method while we compute it by the Lehmann representation (<ref>) with quantities calculated by the VQE-based method in Sec. <ref>.




§ COMPUTATION OF GREEN'S FUNCTION WITH VARIATIONAL QUANTUM SIMULATION 

In this section, we first review the VQS algorithm <cit.> which calculates a quantum state after the time evolution by a given Hamiltonian within parametrized ansatz states.
Next we propose the method to compute the Green's function by extending the original VQS algorithm.



 §.§ Review of variational quantum simulation

Here we review the variational quantum real time simulation algorithm introduced in Ref. <cit.>. Let us consider an ansatz quantum state created by a parametrized quantum circuit,

    |ψ (θ⃗)⟩ = U(θ⃗) |φ_0⟩ = U_N_θ(θ_N_θ)… U_i(θ_i)… U_1(θ_1)|φ_0⟩,

whereU_i(θ_i)is some unitary gate with (real-valued) parameterθ_i,N_θis the number of parameters, and|φ_0⟩is a reference state to create the ansatz state.
We assume thatU_i(θ_i)is composed of a set of Pauli rotation gates,e^iα^(j)θ_i 𝒫^(j)with a coefficientα^(j)and a Pauli matrix𝒫^(j), and other non-parametrized gates.
For a given initial state|ψ(θ⃗(0))⟩and HamiltonianH, the VQS algorithm finds the solution of the Schrödinger equation, 

    d |ψ(θ⃗(t))⟩/dt=-i H |ψ(θ⃗(t))⟩,

within the Hilbert space spanned by the ansatz quantum state,{ |ψ(θ⃗)⟩ }_θ⃗.
Specifically, the time evolution described by Eq. (<ref>) is mapped to the time evolution of parametersθ(t). 
Although there are several variational principles to map Eq. (<ref>) to the equations forθ(t)<cit.>, we choose McLachlan's variational principle <cit.> in this paper because it is the most stable and physically reasonable among them <cit.>.

McLachlan's variational principle  <cit.> maps Eq. (<ref>) to the equation of motion of the parametersθ⃗(t)by minimizing the distance between the exact evolution of the Schrödinger equation and the evolution of the parametrized ansatz state under infinitesimal variation of timeδt<cit.>, 

    min δ‖( ∂/∂ t + iH ) |ψ (θ⃗(t))⟩‖,

where|φ⟩ =⟨φ|φ|$⟩ is the norm of |φ⟩.
One can explicitly write down the equation determining {θ̇_i}_i,

    ∑_j M_i,jθ̇_j=V_i,

where

    M_i,j   =Re(∂⟨ψ(θ⃗(t))|/∂θ_i∂|ψ(θ⃗(t))⟩/∂θ_j),
    
          V_i    =Im(⟨ψ(θ⃗(t))|H ∂|ψ(θ⃗(t))⟩/∂θ_i),

for i=1,…,N_θ.
We note that the matrix M and vector V can be efficiently obtained by measurements of quantum circuits as 
described in Refs. <cit.>.

Finally, from the solution of Eq. (<ref>), one can evaluate the parameter at time t+ δ t as θ⃗(t+δ t) ≈θ⃗(t)+ θ̇⃗̇(t) δ t.
By iterating this procedure from the initial parameters θ⃗(0), we can obtain the time evolution of the parameters θ⃗(t) and the quantum state |ψ (θ⃗(t))⟩.
We note that the VQS algorithm can be viewed as approximating the time evolution operator U(t) = e^-iHt by the variational quantum circuit U(θ⃗(t)), although the approximation is valid only when the operators are applied for the initial state |ψ(θ⃗(0))⟩: U(θ⃗(t)) is optimised to the chosen initial state |ψ(θ⃗(0))⟩ and U(θ⃗(t))|ψ_1⟩≠ e^-iHt|ψ_1⟩ for a different initial state |ψ_1⟩ in general.



 §.§ Computation of Green's function

Now, we discuss how we can apply the VQS algorithm for evaluating the Green's function.
What we want to evaluate is the retarded Green's function for t>0,

    G_k^R (t)     = -i ( ⟨G| e^i H t c_k ↑ e^-i H t c_k ↑^†|G⟩. 
       . + ⟨G|c_k ↑^† e^i H t c_k ↑ e^-i H t|G⟩).

The first and second term can be evaluated similarly, so we focus on the first term.

First, we prepare the (approximate) ground state of a given Hamiltonian H described by N qubits on near-term quantum computers, by using the conventional VQE method <cit.> or the variational imaginary time simulation algorithm <cit.>.
We denote it as |G⟩ = U_G |φ_0⟩ with a unitary U_G and a reference state |φ_0⟩.
Next, we decompose the fermion operator c_k,↑ into a sum of Pauli matrices <cit.>,

    c_k,↑→∑_n=1^N_kλ^(k)_n P_n,   
     c_k,↑^†→∑_n=1^N_kλ^(k)*_n P_n,

where P_n is a tensor product of Pauli matrices acting on N qubits that satisfies P_n^†=P_n,   P_n^2=I_N and λ^(k)_n is a (complex) coefficient.
For example, one can adopt the Jordan-Wigner transformation <cit.> for the decomposition.
The first term of Eq. (<ref>) can be rewritten as 

    ∑_i, jλ^(k)_i λ^(k)*_j ⟨G| e^i H t P_i e^- i H t P_j |G⟩.

Therefore, the problem reduces to finding a way to compute ⟨G| e^i H t P_i e^- i H t P_j |G⟩ on near-term quantum computer.



  §.§.§ Direct method to compute Green's function

A direct way to evaluate ⟨G| e^iHt P_i e^- iHt P_j |G⟩ by the VQS algorithm is as follows.
The time evolution of the states |G⟩ and P_j |G⟩ are approximated on quantum computers as

    e^-iHt|G⟩≈ U^(1)(θ⃗_1(t))  |G⟩,  
    
     e^-iHt P_j |G⟩≈ U^(2)(θ⃗_2(t)) P_j |G⟩,

where U^(1,2)(θ⃗) are parametrized unitary circuits and the initial parameters θ⃗_1(0) and θ⃗_2(0) of the VQS are taken so as to satisfy  U^(1,2)(θ⃗_1,2(0))=I_N. Note that, as U^(1)(θ⃗_1(t)) and U^(2)(θ⃗_2(t)) are optimised for initial states |G⟩ and P_j |G⟩, respectively, generally U^(1)(θ⃗_1(t))≠ U^(2)(θ⃗_2(t)) holds.
Then, ⟨G| e^i H t P_i e^- i H t P_j |G⟩ can be evaluated as the transition amplitude of P_i between e^-iHt|G⟩ and e^-iHt P_j |G⟩,

    ⟨G| e^i H t P_i e^- i H t P_j |G⟩
       ≈   ⟨G|( U^(1)(θ⃗_1(t)) )^† P_i U^(2) (θ⃗_2(t)) P_j |G⟩.

This can be evaluated by using the quantum circuit shown in Fig. <ref>. 

However, this quantum circuit necessitates a huge number of control operations from an ancillary qubit, because two controlled-unitary operations for U^(1)(θ⃗_1(t)) and  U^(2)(θ⃗_2(t)) are present.
Since control operations in near-term quantum computers are costly and have relatively low fidelity in general, the large number of them in the algorithm will deteriorate the performance of computations in real near-term quantum computers.
Therefore, we propose another more efficient method to evaluate ⟨G| e^iHt P_i e^- iHt P_j |G⟩ which will safely run on near-term quantum computers.






  §.§.§ Efficient method to compute Green's function

The problem in the previous method stems from the fact that the variational representations of the time evolution operator U(t)=e^-iHt are different between two initial states, |G⟩ and P_j|G⟩.
If we can construct the variational circuit U(θ⃗(t)) which simultaneously approximates the time evolution operator for the initial states |G⟩ and P_j |G⟩, we obtain

    ⟨G| e^i H t P_i e^- i H t P_j |G⟩
       ≈   ⟨G|( U(θ⃗(t)) )^† P_i U(θ⃗(t)) P_j |G⟩.

In this case, the quantum circuit for evaluating the term is significantly simplified as depicted in Fig. <ref>.
We will now describe how to construct such variational quantum circuit which simultaneously approximates the time evolution operator for general multiple states. 



Here, we consider the most general case where we have L multiple initial states for the time evolution, {|ψ_l⟩}_l=0^L-1.
Let us consider a state with ancilla, 

    |Ψ_0⟩= 1/√(L)∑_l=0^L-1|l⟩_a ⊗|ψ_l⟩_s,

where subscripts a and s denote the ancilla and the system of interest, respectively, and {|l⟩_a}_l is orthonormal state of the ancilla.
We note that when 2^k-1 < L ≤ 2^k, k ancilla qubits are needed.
By using the VQS algorithm to the state |Ψ_0⟩ with the variational quantum circuit I_a ⊗ U_s (θ⃗), we can construct the unitary operator U_s(θ⃗(t)) which approximately behaves as the time evolution operator for {|ψ_l⟩}_l=0^L-1.
To be more concrete, the ansatz state for the VQS algorithm is,

    |Ψ(θ⃗)⟩   := I_a ⊗ U_s  (θ⃗) |Ψ_0⟩
       = 1/√(L)∑_l=0^L-1|l⟩_a ⊗|ψ_l(θ⃗)⟩_s,

where we define |ψ_l(θ⃗)⟩_s= U_s (θ⃗)|ψ_l⟩_s.
The M matrix and V vector in Eq. (<ref>) become 

    M_i,j   = Re( ∂⟨Ψ(θ⃗(t))|/∂θ_i∂|Ψ(θ⃗(t))⟩/∂θ_j) 
       =1/L∑_l=0^L-1Re(∂⟨ψ_l(θ⃗(t))|_s/∂θ_i∂|ψ_l(θ⃗(t))⟩_s/∂θ_j),

and 

    V_i    = Im(⟨Ψ(θ⃗(t))| H ∂|Ψ(θ⃗(t))⟩/∂θ_i) 
       =1/L∑_l=0^L-1Im(⟨ψ_l(θ⃗(t))|_s H ∂|ψ_l(θ⃗(t))⟩_s/∂θ_i).

From this expression, one notice that this algorithm minimizes the average of δ‖( ∂/ ∂ t + iH ) |ψ_l (θ⃗(t))⟩‖ for l=0,…,L-1.
We also note that the algorithm itself can run without resorting to the ancilla because each summand in Eqs. (<ref>), (<ref>) can be computed by a distinct quantum circuit:
one can compute each term in different run of quantum computers and sum up them by classical computers.
The advantage of using the ancilla is that we can compute M and V for exponentially increasing number of input states in terms of the number of ancilla qubits. For example, when L=4 one should sum up results of four runs of quantum computers without the ancilla whereas one run is necessary with the ancilla (accompanying with the drawback of the complicated quantum circuit).
We remark that the evaluation of the transition amplitude between the time-evolved states |ψ_l(θ⃗)⟩_s= U_s (θ⃗)|ψ_l⟩_s requires some ancilla qubits in general such as Fig. <ref>.

In the case of calculation of the Green's function, the initial states are |ψ_0⟩ = |G⟩ and |ψ_1⟩ = P_j|G⟩.
The ansatz quantum state (<ref>) for the VQS algorithm can be constructed by a quantum circuit shown in Fig. <ref>.



To sum up, the calculation of each term in Eq. (<ref>), and consequently the Green's function, with the VQS algorithm proceeds as follows:

 
  0. Prepare the approximate ground state of a given Hamiltonian H by conventional methods on near-term quantum computers, such as the VQE <cit.>.
 We denote the ground state as |G⟩=U_G|φ_0⟩.
 
  1. Construct the variational quantum circuit U(θ⃗) which approximates the time evolution e^-iHt for two initial states |G⟩ and P_j|G⟩. The VQS algorithm with the circuit shown in Fig. <ref> will find such U(θ⃗).  
 
  2. Evaluate ⟨G|U(θ⃗)^† P_i U(θ⃗) P_j |G⟩ by the quantum circuit shown in Fig.<ref>.


Finally, we present a detailed resource estimation about the number of required distinct runs of the quantum circuits for this algorithm in Appendix <ref>. 




§ COMPUTATION OF GREEN'S FUNCTION WITH EXCITED-STATES SEARCH ALGORITHM 

In this section, we describe another method to compute the Green's function of a given quantum system.
We compute the energy eigenstates and transition amplitudes of fermion operators by the algorithm based on the SSVQE method <cit.> and the MCVQE method <cit.>, and take advantage of the Lehmann representation of the spectral function (<ref>).
We discuss two types of algorithms for calculating the excited states and the transition amplitudes.
The first one is based on the SSVQE algorithm with different weights where one obtains the excited states directly on quantum computers, while the second one is based on the SSVQE algorithm with identical weights where some classical post-processing after the use of quantum computers are required.
The computation of the first algorithm is simpler than that of the second one, but the convergence of the algorithm is better for the second one in general.
We note that the essential part of the algorithm described in this section is already discussed in Refs. <cit.>, so our contribution will be application of it for calculation of the Green's function.



 §.§ Computation by SSVQE with different weights

Let us consider finding K smallest eigenenergies and eigenstates of a given Hamiltonian.
The SSVQE algorithm with different weights finds the variational quantum circuit U(θ⃗^*) which makes input orthonormal states {|ψ_j⟩_j=0^K-1} into approximate eigenstates of H, |E_j⟩ := U(θ⃗^*)|ψ_j⟩.
The approximate eigenenergies are obtained as E_j := ⟨E_j| H |E_j⟩ = ⟨ψ_j| U(θ⃗^*)^† H U(θ⃗^*)|ψ_j⟩.
The algorithm performs this task by minimizing the cost function

    𝒞_0(θ⃗)= ∑_j=0^K-1 w_j ⟨ψ_j|U(θ⃗)^† H U(θ⃗)|ψ_j⟩,

with respect to the parameters θ⃗, where w_0>...>w_K-1>0 are weights which assure the approximate eigenstates {E_j }_j have the ascending order, E_0 ≤ ... ≤E_K-1.
After convergence of the classical minimization for 𝒞_0(θ⃗), one obtains optimal parameters θ⃗^* and can compute {E_j, |E_j⟩}_j.

To compute the Lehmann representation of the spectral function (<ref>), we also need the transition amplitude of the fermions c_k,↑, c_k,↑^†, such as ⟨E_j|c_k,↑|E_k|$⟩.
In general, the evaluation of the transition amplitude between different quantum states needs a complicated quantum circuit,
but in this case the evaluation will be done in a simple way due to the fact that|E_j⟩'s are created from the same unitary gateU(θ⃗^*).
Specifically, if we can easily make superpositions of the input states,|ψ_j⟩and|ψ_k⟩,⟨E_j|c_k,↑|E_k|$⟩ can be evaluated by simply taking the expectation value of c_k,↑ for several superpositions.
To see this, first we map c_k ↑ into qubit operators like Eq. (<ref>) and decompose it into real part and imaginary part, c_k ↑=A_k+iB_k, where A_k and B_k are hermitian operators.
Then we have

    ⟨ψ_j”| U(θ⃗^*)^† c_k ↑  U(θ⃗^*)|ψ_j'⟩ = ⟨ψ_j”| U(θ⃗^*)^† A_k  U(θ⃗^*)|ψ_j'⟩
     
    +i ⟨ψ_j”| U(θ⃗^*)^† B_k  U(θ⃗^*)|ψ_j'⟩.

Each term can be evaluated by using |ψ^±_j',j”⟩= U(θ⃗^*)(|ψ_j'⟩±|ψ_j”⟩)/√(2) and |ψ^i ±_j',j”⟩ = U(θ⃗^*)(|ψ_j'⟩± i |ψ_j”⟩)/√(2) as

    Re(⟨ψ_j”| U(θ⃗^*)^† A_k  U(θ⃗^*)|ψ_j'⟩)    = ⟨ψ^+_j',j”|A_k |ψ^+_j',j”⟩
       -⟨ψ^-_j',j”| A_k |ψ^-_j',j”⟩
    Im(⟨ψ_j”| U(θ⃗^*)^† A_k  U(θ⃗^*)|ψ_j'⟩)   =
    ⟨ψ^i+_j',j”|A_k |ψ^i+_j',j”⟩
       -⟨ψ^i-_j',j”|A_k |ψ^i-_j',j”⟩,

and similar equations for the B_k term.

In typical situations, the input states are taken as simple states, e.g., computational basis, so preparing superpositions of them on quantum computers is not so difficult.
Therefore, by substituting eigenenergies of Eq. (<ref>) with E_j and the transition amplitudes with ⟨E_0|c_k,↑|E_k|$⟩, we can evaluate the spectral function, and the Green's function accordingly. 



 §.§ Computation by SSVQE with identical weights

Next, we introduce another type of algorithms to obtain excited states and the transition amplitude between them.
This algorithm combines the SSVQE algorithm with identical weights and the quantum subspace expansion method <cit.>, and is essentially the same as the MCVQE algorithm <cit.>.

The procedure of the algorithm is as follows.
First, we prepare orthonormal  input states{ |ψ_j⟩_j=0^K-1 }, which are simple and easy to realize a superposition of them on quantum computers.
Then we minimise the cost function

    𝒞_1(θ⃗)= ∑_j=0^K-1⟨ψ_j|U^†(θ⃗)H U(θ⃗)|ψ_j⟩,

with respect to parametersθ⃗, whereU(θ⃗)is the ansatz quantum circuit.
After the optimisation, the subspace spanned by{U(θ⃗^*)|ψ_j⟩}_j=0^K-1will be close to that spanned by the trueKeigenstates {|E_j⟩ }_j=0^K-1, whereθ⃗^*is the parameters after optimisation.
At this stage,U(θ⃗^*)|ψ_j⟩is generally the superposition of the excited states{|E_j⟩ }_j=0^K-1and⟨ψ_j| U(θ⃗^*)^†H U(θ⃗^*)|ψ_j⟩is not a good approximation to the true eigenvalueE_j.
Therefore, to obtain the eigenstates and eigenvalues ofH, we solve the eigenvalue problem within the subspace spanned by{U(θ⃗^*)|ψ_j⟩}_j=0^K-1,

    ℋ𝒱= 𝒱ℰ,

whereℋ_i,j=⟨ψ_i|U^†(θ⃗^*)H U(θ⃗^*)|ψ_j⟩,𝒱isK×Kmatrix containing eigenvectors as its columns, andℰis a diagonal matrix whose diagonal elements are eigenvalues.
The approximatem-th excited state|E'_m⟩is expressed as 

    |E'_m⟩=∑_j 𝒱_j, m U(θ⃗^*)|ψ_j⟩,

and the approximate eigenenergies appear asE'_j = ℰ_j,j.

The transition amplitudeC_m,n^(k)=⟨E'_m| c_k ↑ |E'_n⟩can be computed as 

    C_m,n^(k)= ∑_j', j”𝒱_j”, m^* 𝒱_j', n⟨ψ_j”| U(θ⃗^*)^† c_k ↑  U(θ⃗^*)|ψ_j'⟩.

The quantity⟨ψ_j”| U(θ⃗^*)^†c_k ↑   U(θ⃗^*)|ψ_j'⟩can be evaluated in the way described in the previous subsection.
Thus, we can compute the transition matrixC^(k)_m,n, and evaluate the spectral function by Eq. (<ref>) and the Green's function.







§ NUMERICAL DEMONSTRATION 

In this section, we numerically demonstrate our proposed methods to calculate the Green's function and the spectral function at zero temperature.
We consider the Fermi-Hubbard model of two sites with the particle-hole symmetry,

    H    =    -t ∑_σ=↑,↓( c_1,σ^† c_2,σ + h.c.)
     + U ∑_i=1^2 c_i,↑^† c_i,↑ c_i,↓^† c_i,↓
           -U/2∑_i=1,2,σ=↑,↓ c_i,σ^† c_i,σ

wheretis a parameter characterizing hopping between the sites andUdenotes the strength of the on-site Coulomb repulsion <cit.>.
We set the hopping parametert=1throughout this paper.
We simulate two proposed protocols in Secs. <ref> and <ref> by classical computers with the fast quantum circuit simulation library Qulacs <cit.>.
We use the Jordan-Wigner transformation <cit.> to map the fermionic Hamiltonian (<ref>) into the qubit one with four qubits by using the library OpenFermion <cit.>.



 §.§ Numerical simulation of the method based on variational quantum simulation


We calculate the real-time Green's function of the model (<ref>) at zero temperature by using the method described in Sec. <ref>.
First, we prepare the ground state of the model (<ref>) by the standard VQE algorithm with a hardware-efficient-type ansatz <cit.> depicted in Fig. <ref>.
Then, we perform the VQS algorithm.
As an ansatz quantum state, we adopt the so-called variational Hamiltonian ansatz <cit.> inspired by the Suzuki-Trotter decomposition of the time evolution operatore^-iHt.
The variational Hamiltonian ansatz is defined through the qubit representation of the Hamiltonian,

    H_qubit = ∑_m c_m P_m,

whereP_mis a (multi-qubit) Pauli matrix andc_mis a coefficient.
An ansatz state for the variational Hamiltonian ansatz is given by|ψ(θ)⟩ = U_VHA(θ) |φ_0⟩, where

    U_VHA({θ_m^(d)}_m,d) = ∏_d=1^n_d( ∏_m exp( i θ_m^(d) P_m ) ),

andn_ddenotes the depth of the ansatz.
We note that we remove the identity operator from Eq. (<ref>) when constructing the ansatz.
If the qubit representation of the Hamiltonian (<ref>) hasN_Pterms (except for the identity operator), the number of parameters of the ansatz will beN_P n_d.
In our simulation,N_P=6and we choosen_d=8, so there are 48 parameters in the parametrized quantum circuit whereas general unitary operators on the system have(2^4)^2 = 256parameters. 
Further details on numerical calculations are described in Appendix <ref>.

The result is shown in Fig. <ref>.
The VQS algorithm nicely reproduces the exact dynamics (the panels in left and center columns), and the spectral function (right columns).
These figures illustrate the possibility of the VQS algorithm proposed in this study to calculate the Green's function.
In section <ref>, dependence of the results on the depth of the ansatz is analyzed.
Furthermore, numerical simulations forn_d=4and the four-site Fermi-Hubbard model are presented in Appendix <ref>. One can see that numerical results forn_d=4have almost the identical performance compared to the case ofn_d=8; therefore depending on strength of physical noises, the choice ofn_d=4may be recommended to avoid accumulation of physical errors.



 §.§ Numerical simulation of the method based on excited-state search


Next, we numerically simulate the method described in Sec. <ref>.
We adopt the symmetry-preserving ansatz <cit.> drawn in Fig. <ref>, which preserves the total number of particles in the system.
We use the SSVQE algorithm with identical weight and calculate five energy eigenstates of the model (<ref>).
As an input states, we simply choose the computational basis states with desired particle number: to calculate the particle (hole) part of the spectrum function, we choose|φ_0⟩=|0011⟩for the ground state,|φ_1⟩=|0001⟩, |φ_2⟩=|0010⟩, |φ_3⟩=|0100⟩, |φ_4⟩=|1000⟩(|φ_1⟩=|0111⟩, |φ_2⟩=|1011⟩, |φ_3⟩=|1101⟩, |φ_4⟩=|1110⟩) for the excited states. 

Figure <ref> shows the result of numerical simulation.
The SSVQE algorithm almost perfectly reproduces the exact result obtained by exact diagonalization.



§  DEPENDENCE OF ACCURACY OF THE VARIATIONAL QUANTUM SIMULATION ON DEPTH OF THE ANSATZ 

In this section, we provide a systematic analysis on how the accuracy of the numerical simulations of the variational quantum simulation (VQS) in Sec. <ref> depends on the depthn_dof the Hamiltonian ansatz (Eq. (<ref>)).
We run numerical simulations forn_d=4,5,⋯,10under the same conditions as in Fig. <ref>.

To see the accuracy of the simulations quantitatively, we calculate the mean absolute error (MAE) of the spectrum function atk=πin the region ofω∈[-5, 5],

    Δ_E (k) = 1/2N_ω + 1   ∑_n=-N_ω^N_ω| A^exact(k, ω = 5n/N_ω) 
       - A^VQS(k, ω = 5n/N_ω) |,

where2N_ω+ 1is the total number of data points andA^exact(VQS)(k,ω)is the spectrum function calculated by exact diagonalization (VQS).
We takeN_ω=5000.
Interestingly, as seen from Fig. <ref>, the MAE decreases with the inverse of the depth of the ansatz.
This dependence reminds us of the Suzuki-Trotter decomposition of the time evolution operator, i.e.,

    U(t) = e^-iH_qubitt≈ U^Trot_n_d(t) = (∏_m e^-i c_m P_m ·t/n_d)^n_d

with the error ofO(t^2/n_d), where we have used the notation in Eq. (<ref>).
In Fig. <ref>, we also show the MAE of the spectrum function calculated by the dynamics obtained from the approximate time evolution operator Eq. (<ref>) with the same time step and duration used for the VQS.
The MAE for this case exhibits1/n_ddependence as expected, but the values of the MAE are much larger than those for the VQS; the slope of the fit of the MAE with1/n_dis about six times smaller for the VQS  (1.820/0.285≈6.4).
The result shown in Fig. <ref> not only provides an estimation of errors in the VQS calculation when the ansatzes with various depths are used, but also illustrates the practical advantage of employing the VQS compared with the Suzuki-Trotter decomposition of the time evolution operator which has the same-depth quantum circuit.





§ FEASIBILITY OF QUANTUM ALGORITHMS ON NEAR-TERM QUANTUM COMPUTERS 

In this section, we discuss the feasibility of implementing our proposed algorithms on near-term quantum computers.
Let us consider 25 sites two-dimensional Fermi-Hubbard model on a square lattice, whose simulation requires 50 qubits and is almost intractable for classical computers.
The model is defined as

    H_2d = -t ∑_⟨i,j|,⟩σ( c_i,σ^† c_j,σ + h.c.)
     + U ∑_i c_i,↑^† c_i,↑ c_i,↓^† c_i,↓,

where⟨i,j|$⟩ runs nearest-neighbor sites on  a square lattice and σ=↑,↓ denotes the spin.
Based on the argument in Ref. <cit.>, using the Hamiltonian ansatz (Eq. (<ref>)) with the Jordan-Wigner transformation, we need N_two≈ 1000 two-qubit gates per depth of the ansatz when we employ the rotation Z gate and partial swap gates as elementary gates.
Although near-term quantum computers contain inevitable noise in gate operations, the technique of quantum error mitigation  <cit.> can suppress errors and recover noiseless results with a reasonable overhead when the error rate per gate ϵ_gate satisfies N_gateϵ_gate≲ 2, where N_gate is the number of gates. Therefore, when we adopt the Hamiltonian ansatz of two depths as an ansatz quantum circuit for the VQS or the SSVQE, we need at least 2/2000 = 0.1 % for the error rate of two-qubit gates, which has been achieved in current experimental setups <cit.>. 
We note that the error of single-qubit gate is ignored here because it is negligible compared with that of two-qubit gate.
Only a few additional controlled operations in the VQS-based algorithm introduced in Sec. <ref> is also be neglected.
In addition, for the VQS-based algorithm, gates for the unitary U_G which prepares the ground state of the system must be taken into account. The more detailed argument has been made in Appendix <ref>.




§ EXTENSION TO FINITE TEMPERATURE GREEN'S FUNCTION 

Here, we discuss the extension of our proposed methods to the Green's function at finite temperature.
For finite temperature T>0 or inverse temperature β = 1/T < ∞,
the retarded Green's function is defined as

    G^R_k(t;β)    = 1/Z(β)∑_n e^-β E_n G^R_k,n(t), 
    
    Z(β)    = Tr( e^-β H), 
    
    G^R_k,n(t)     = -i Θ(t) ( ⟨E_n| e^i H t c_k ↑ e^-i H t c_k ↑^†|E_n⟩. 
       . + ⟨E_n|c_k ↑^† e^i H t c_k ↑ e^-i H t|E_n⟩),

where |E_n⟩ and E_n denote the eigenstates and eigenvalues of Hamiltonian H, respectively.
The corresponding spectral function is

    G^R_k(ω;β)     = ∫_-∞^∞ dt   e^i(ω+iη)t G^R_k(t; β),
    
    A_k(ω; β)    = 1/Z(β)∑_n,m e^-β E_m
       ×(  |⟨E_n|c_k^†|E_m||⟩^2 /ω+E_m-E_n+iη.
      . +  |⟨E_n|c_k|E_m||⟩^2 /ω-E_m+E_n+iη).




 §.§ Variational quantum simulation for Green's function at finite temperature





Equation (<ref>) can be evaluated on quantum computers by combining the VQS method and the thermofield double technique which purifies the Gibbs state of the system.
The procedure is essentially the same as the finite-temperature version of the density matrix renormalization group method <cit.>.
We note that several methods based on the typically of the Hilbert space to evaluate physical quantities at finite temperature <cit.> might also be combined with the VQS algorithm.

Let us consider a N-qubit “environment" system (denoted by subscript e) in addition to the original N-qubit system of interest (denoted by subscript s).
First We prepare a state |Φ_0⟩ which satisfies Tr_e (|Φ_0⟩⟨Φ_0|) = I_s.
For example, we choose |Φ_0⟩ = 1/√(2^N)∑_i=0^2^N-1|i⟩_s ⊗|i⟩_e where |i⟩_s,e is the computational basis of the system and environment.
By using the method of variational imaginary time evolution introduced in Ref. <cit.>, one can obtain the state |Φ(β)⟩≈ Z(β)^-1/2 e^-β H/2|Φ_0⟩.
Namely, the variational imaginary time evolution for the total system with the “Hamiltonian" H ⊗ I_e and the variational quantum circuit drawn in Fig. <ref> will produce |Φ(β)⟩.
We note that |Φ(β)⟩ satisfies

    Tr_e ( |Φ(β)⟩⟨Φ(β)|) ≈1/Z(β) e^-β H.


Next, we perform the same VQS algorithm for the Green's function at zero temperature by replacing 
|G⟩ with |Φ(β)⟩.
It will bring out the variational quantum circuit on the original system U_s(θ⃗) satisfying

    ⟨Φ(β)| U_s(θ⃗)^† P_i^(s) U_s(θ⃗) P_j^(s)|Φ(β)⟩
       ≈⟨Φ(β)| e^iHt P_i^(s) e^-iHt P_j^(s)|Φ(β)⟩
       = 1/Z(β)∑_n e^-β E_n⟨E_n| e^iHt P_i^(s) e^-iHt P_j^(s)|E_n⟩,

where superscript of the Pauli operator P_i^(s) implies that it only acts on the system s.
In Fig. <ref>, we show the ansatz quantum circuit to construct the unitary gate U_s(θ⃗(t)).
By substiting the above quantity into Eq. (<ref>), we can evaluate (each term of) G_k(t;β).

We finally remark on another way to evaluate Eq. (<ref>) based on the VQS algorithm.
It is possible to obtain several approximate eigenenergies {E_n}_n=1^K and eigenstates {|E_n⟩}_n=1^K of the system by the SSVQE or MCVQE algorithm, and to perform the VQS algorithm in Sec. <ref> for each obtained eigenstate |E_n⟩.
The result approximates G^R _kn(t) in Eq. (<ref>), so substituting it as well as Ẽ_n into Eq. (<ref>) will give the Green's function at finite temperature (with truncating the summation up to n=K).
This type of the algorithm dose not need the environment qubits, but the number of energy eigenstates K to evaluate the Green's function with fixed accuracy would exponentially increase as the size of the system and the inverse temperature β.



 §.§ Computing Green's function at finite temperature with excited-states search method

Extension of the algorithm in Sec. <ref> to the finite temperature Green's function is rather simple.
Since we already introduce a way to evaluate the quantity |⟨E_n|c_k,↑|E_m||⟩^2, putting them into Eq. (<ref>) gives the spectral function at finite temperature with truncating the summation up to n=K, where K is the number of eigenstates obtained by the excited-states search algorithms.
Again, this method also has a potential problem of the exponentially increasing number of eigenstates required to compute the Green's function with fixed accuracy.





§ DISCUSSION AND CONCLUSION 

In this paper, we proposed two methods for calculating the Green's function compatible with NISQ hardwares.
One of the proposed methods uses the conventional VQE to prepare a ground state, and directly calculate the real time retarded Green's function for the obtained ground state by the VQS algorithm. We introduced a method for constructing a variational quantum circuit which acts as the time evolution operator for multiple initial states simultaneously, and makes the quantum circuit for computation of the Green's function significantly shallower.
Note that this method can be straightforwardly applied to evaluation of the linear response function <cit.> expressed as

    ϕ_BA(t-t')   =i ⟨B(t-t'), A(0)|_⟩0, 
    
    B(t-t')    = e^iH (t-t') B  e^-iH (t-t')

where A is an observable coupled to external field, for example, magnetic moment or charge density, B an observable to be measured. 
The other proposed method evaluates the transition amplitude of the fermion operators between energy eigenstates of the system by exploiting either the SSVQE or MCVQE method, and computes the spectral function and the Green's function with the use of the Lehmann representation.

In numerical simulation for the Green's function at zero temperature, both methods successfully reproduced the spectral function of the 2-site Fermi-Hubbard model.
We here discuss possible causes to hinder or deteriorate the performance of the methods for general large systems.
For the first method by the VQS algorithm, the choice of the ansatz for the real time evolution is crucial: once the variational quantum state is out of the correct trajectory in the Hilbert space, it is rare that the state returns to it.
Therefore the ansatz has to be chosen carefully so that the simulated quantum state has remained in a correct trajectory.
As for the second method by the excited-states search methods, it is in general unclear how many excited states are required to reach the desired precision of the Green's function.
It will possibly grow exponentially as the inverse temperature β and the size of the system N.
We leave the investigation of these problems and further comparison of two methods proposed in this study to future work. 

Since the Green's function is fundamental to study the nature of quantum systems, we believe our study will extend the possibility to utilize near-term quantum computers in condensed matter physics, quantum chemistry and materials science.




Note added:
Recently, a paper discussing calculation of the Green's function with NISQ devices appeared <cit.>.
The method used in that paper is based on the Lehmann representation of the Green's function, and the authors calculated excited states of the system by performing the VQE with penalty terms which are proportional to the overlaps between previously-found low excited states <cit.>.



A part of the numerical simulations in this work were done on Microsoft Azure Virtual Machines provided through the program Microsoft for Startups.
YON and SE acknowledge valuable discussion with Kosuke Mitarai, Xiao Yuan, Sam McArdle, Zhenyu Cai and Takamichi Umetsu. 
IK was supported by Qunasys Inc.
This work was also supported by MEXT Q-LEAP JPMXS0118068682.












§  RESOURCE ESTIMATION OF THE ALGORITHM BASED ON THE VARIATIONAL QUANTUM SIMULATION 



In this section, we discuss the number of distinct runs of quantum circuits required to implement our algorithm based on the VQS introduced in Sec. <ref>.
We firstly show that the error of the Green's function can be upper bounded by a product of the operator norm of observables involved in the Green's function and the trace distance between the ideal time-evolved state and the state computed by the VQS.
Next, based on the argument in Ref. <cit.>, we discuss the errors of the real time Green's function and clarify two sources of error: algorithmic error and implementation error.
Finally, we estimate the number of total distinct runs of quantum circuits needed to achieve a certain accuracy in computing the Green's function in the frequency domain, i.e., the spectral function. 




 §.§ Relationship between the error of the Green's function and the trace distance




  §.§.§ Difference of expectation values of an observable between two states

First, we prove the difference of expectation values of an observable M for two distinct quantum states (density matrices) ρ and σ can be upper bounded by the product of the operator norm of M and the trace distance of the two states.
It follows that

    Tr(M(ρ-σ))   = ∑_k m_k Tr( E_k (ρ-σ) ) 
       ≤∑_k |m_k| | Tr( E_k (ρ-σ) )|  
       ≤∑_k M Tr(E_k |ρ-σ |) = 2 M  D(ρ, σ),

where M is an operator norm of M (the largest singular value of M), and M=∑_k m_k E_k is the spectral decomposition of M (m_k is an eigenvalue of M and E_k is the corresponding projector).
We used |Tr(E_k (ρ-σ))| ≤Tr(E_k |ρ-σ| <cit.> and ∑_k E_k =I. 



  §.§.§ Error of the Green's function as the trace distance of the ideal state and the trial state
 
Next, we show the trace distance between the exact state after time evolution and the state obtained by the VQS gives the upper bound of the error of the calculated Green's function in our algorithm.
For simplicity, we only consider t>0.

Let us introduce two wavefunctions as

    |Ψ(θ⃗(t))⟩   = 1/√(2)(|0⟩_a ⊗ U(θ⃗(t)) |G⟩_s+|1⟩_a⊗ U(θ⃗(t)) P_j |G⟩_s ), 
    |Φ(t)⟩   = 1/√(2)(|0⟩_a ⊗ e^-iHt|G⟩_s+|1⟩_a⊗ e^-iHt P_j |G⟩_s ),

where a subscript a denotes an ancilla qubit and |G⟩_s is the ground state of the system.


In our algorithm based on the VQS, the Green's function is calculated based on the following decomposition (Eq. (<ref>) in the main text),

    G_k^R(t) = ∑_i,jλ_i^(k)λ_j^(k)*⟨G|e^-iHt P_i e^-iHtP_j |G⟩≡∑_i,jλ_i^(k)λ_j^(k)* G_ij^(0)(t).

The quantity G_ij^(0)(t) is approximated by 

    G_i,j(θ⃗(t)) ≡⟨Ψ(θ⃗(t))| X_a P_i |Ψ(θ⃗(t))⟩+i ⟨Ψ(θ⃗(t))| Y_a P_i |Ψ(θ⃗(t))⟩,

where the measurement of X_a P_i and Y_a P_i correspond to the cases we set ϕ=0, π/2 in Fig. <ref> in the main text, respectively.
By using Eq. (<ref>), the error of this approximation is upper bounded as

    ϵ_i,j≡| G̃_i,j(θ⃗(t)) - G_i,j^(0)(t) |≲ 2√(2) D(|Ψ(θ⃗(t))⟩,|Φ(t)⟩)+√(2)/√(N_m),

where G̃_i,j(θ⃗(t)) is an experimentally obtained value of G_i,j(θ⃗(t)) according to Eq. (<ref>) and we denote the trace distance between two pure states |φ⟩⟨φ| and |φ'⟩⟨φ'| as D(|φ⟩,|φ'⟩).
The second term in the right hand side indicates the shot noise, where N_m is the number of measurements for X_a P_i and Y_a P_i.
Finally, we can describe the error of the Green's function as follows:

    ϵ^R≡∑_i,j |λ_i^(k)| |λ_j^(k)| ϵ_i,j   ≤∑_i,j |λ_i^(k)| |λ_j^(k)| ( 2√(2) D(|Ψ(θ⃗)⟩,|Φ(t)⟩_0)+√(2)/√(N_m)) 
       = α( 2D(|Ψ(θ⃗(t))⟩,|Φ(t)⟩_0)+1/√(N_m)),

where α=√(2)∑_i,j |λ_i^(k)| |λ_j^(k)|.





 §.§ Algorithmic and implementation errors

We now analyze sources of errors of the real time Green's function in this subsection and clarify two of main contributions of them: algorithmic and implementation errors.
According to Ref. <cit.>, we can upper bound the trace distance between the ideal time-evolved state and the simulated state by the VQS as

    D(|Ψ(θ⃗(T))⟩,|Φ(T)⟩_0) ≤ D(|Ψ(θ⃗(0))⟩,|Φ(0)⟩_0 ) + ∑_n=1^N_step D(U|Ψ(θ⃗((n-1)δ t))⟩, |Ψ(θ⃗(n δ t))⟩),

where T is the time to be simulated, U is the exact time evolution operator of small time step δ t described as U=exp(-iH δ t), and N_step=T/δ t is the number of time steps.
The first term corresponds to the state preparation error of the ground state via the variational quantum eigensolver (VQE), which we denote as ϵ_s=D(|Ψ(θ⃗(0))⟩,|Φ(0)⟩_0 ).
The second term, which accumulates an error for each time step, can be decomposed into two types of errors: algorithmic and implementation errors.
Algorithmic errors contains error due to imperfection of the ansatz approximating the trial state and one due to a finite time step.
Implementation errors are caused by physical errors stemming from the noisy nature of near-term quantum computers and the shot noise.
We assume physical errors can be ignored because it can be suppressed by using quantum error mitigation techniques <cit.>, i.e., we consider only the shot noise as the implementation error.
By using the triangle inequality, it follows that

    D(U|Ψ(θ⃗((n-1)δ t))⟩, |Ψ(θ⃗(n δ t))⟩)  ≤ D(U|Ψ(θ⃗((n-1)δ t))⟩, |Ψ^(0)(θ⃗(n δ t))⟩) + D(|Ψ(θ⃗(n δ t))⟩, |Ψ^(0)(θ⃗(n δ t))⟩),

 where |Ψ^(0)⟩ denotes the state without implementation errors.
 The first term in the right hand side corresponds to the algorithmic error and the second term does to the implementation error for each time step.
 We denote these errors δϵ_A(n) and δϵ_I(n), respectively. 
 



  §.§.§ Algorithmic error

The algorithmic error can be written as <cit.> 

    δϵ_A(n) ≡ D(U|Ψ(θ⃗(n-1)δ t))⟩, |Ψ^(0)(θ⃗(nδ t))⟩) = √(Δ ^(2)_n δ t^2 +Δ ^(3)_n δ t^3 + O(δ t^4)),

where 

    Δ^(2)_n    = ⟨δΨ(θ⃗(n δ t))|δΨ(θ⃗(n δ t))|-⟩|⟨δΨ(θ⃗(n δ t))| Ψ(θ⃗((n-1)δ t))|⟩|^2, 
    |δΨ(θ⃗(n δ t))⟩   =-iH |Ψ(θ⃗((n-1)δ t))⟩-∑_k θ̇_k^(0)∂|Ψ(θ⃗((n-1)δ t))⟩/∂θ_k, 
    Δ^(3)_n    = HH^2 + 1/3H^3 + dR̃/dtd^2R̃/dt^2 + 1/3d^3R̃/dt^3
          + H^2(dR̃/dt^2 + d^2R̃/dt^2) + ( H + H^2) dR̃/dt,

the matrix norm …is induced by the vector norm, the operator d/dt is defined as d/dt = θ̇⃗̇(t) ·∂/∂θ⃗, and R̃≡ = I_a ⊗ U(θ⃗(t)) is the unitary circuit consisting of the ansatz (Eq. (<ref>)). 
The term proportional to Δ^(2)_n stems from limited representative capability of the ansatz state |Ψ^(0)(θ⃗)⟩ approximating the true state.
The term proportional to Δ^(3)_n is, on the other hand, due to a finite time step δ t. 
Therefore, the total accumulation of algorithmic errors from t=0 to t=T is 

    ϵ_A= ∑_n=1^N_stepδϵ_A(n) ≲√(Δ_max^(2))T +√(Δ_max^(3)δ t),

where Δ^(2)_max = max_n Δ^(2)_n and Δ^(3)_max = max_n Δ^(3)_n.
The first term in the right hand side does not depend on the time step δ t because this is the error caused by imperfections of the ansatz to represent the quantum state after the time evolution, whereas the second term can be suppressed by taking the small time step.


Here, we also show another important property of our algorithm based on the VQS.
By denoting |Ψ(θ⃗(t))⟩ as |Ψ(θ⃗(t))⟩ = 1/√(2)(|0⟩_a ⊗|ψ_1(θ⃗(t))⟩+|1⟩_a ⊗|ψ_2(θ⃗(t))⟩) for simplicity, we have

    Δ^(2)_n    =1/2(⟨δψ_1(θ⃗(n δ t))| δψ_1(θ⃗(n δ t))|+⟩⟨δψ_2(θ⃗(n δ t))| δψ_2(θ⃗(n δ t))|⟩) 
       -1/4( |⟨δψ_1(θ⃗(n δ t))| ψ_1(θ⃗(n δ t))|⟩|^2 + |⟨δψ_2(θ⃗(n δ t))| ψ_2(θ⃗(n δ t))|⟩|^2 )
       -1/2Re(⟨δψ_1(θ⃗(n δ t))| ψ_1(θ⃗(n δ t))|⟨%s|%s⟩⟩δψ_2(θ⃗(n δ t))| ψ_2(θ⃗(n δ t))) 
       =1/2(Δ^(2)_n, ψ_1+Δ^(2)_n, ψ_2) + 1/4|⟨δψ_1(θ⃗(n δ t))|ψ_1(θ⃗(n δ t))|-⟩⟨δψ_2(θ⃗(n δ t))|ψ_2(θ⃗(n δ t))|⟩|^2,

where 

    Δ^(2)_n, ψ_i   = ⟨δψ_i(θ⃗(n δ t))|δψ_i(θ⃗(n δ t))|-⟩|⟨δψ_i(θ⃗(n δ t))| ψ_i(θ⃗((n-1)δ t))|⟩|^2, 
    |δψ_i(θ⃗(n δ t))⟩   =-iH |ψ_i(θ⃗((n-1)δ t))⟩-∑_k θ̇_k^(0)∂|ψ_i(θ⃗((n-1)δ t))⟩/∂θ_k.

 From the expression above, one can see that Δ^(2)_n for |Ψ⟩ is not just a mean of Δ^(2)_ns for |ψ_1⟩ and |ψ_2⟩ but also includes the term δ_12 = 1/4|⟨δψ_1|ψ_1|-⟩⟨δψ_2|ψ_2||⟩^2.
 This term, δ_12, indicates that the algorithmic error due to the insufficient ansatz increases when we  want to find a unitary operator to evolve two input states simultaneously, compared with finding different unitaries for each state.
The term δ_12 clarifies the drawback of our algorithm to find a unitary which simultaneously approximates the time evolution operator for multiple states.


However, as seen in numerical simulation in Sec. <ref>, we successfully reproduced the real time Green's functions of the two site fermi-Hubbard model, which implies that the effect of δ_12 can be made small.
In general, the cost for employing different unitaries for multiple states and the circuit like in Fig. <ref>
is much larger than employing a single unitary and the circuit like in Fig. <ref>.




  §.§.§ Implementation error

The implementation error is caused by the shot noise and give errors in the solution of Eq. (<ref>) in the main text,

    M θ̇⃗̇=V⃗,

where 

    M_i,j    = . Re(∂⟨Ψ(θ⃗)|/∂θ_i∂|Ψ(θ⃗)⟩/∂θ_j) |_θ⃗=θ⃗(nδ t), 
    
    V_j   = . Im(⟨Ψ(θ⃗)|H ∂|Ψ(θ⃗)⟩/∂θ_j) |_θ⃗=θ⃗(nδ t).

At each time step, the matrix M and the vector V are obtained by measuring outputs of appropriate quantum circuits <cit.>, which include the shot noise.
We denote error-free values of the matrix M and the vector V⃗ as M_0 and V⃗_0, and observed values of them as M_0 + δ M and V⃗_0 + δV⃗, respectively.
The solution of Eq. (<ref>) calculated by M_0 and V⃗_0 is denoted by θ̇⃗̇_0 and the one calculated by M_0+δ M and V⃗_0 + δV⃗ as θ̇⃗̇_0 + δθ̇⃗̇.
In the first order of δ M_0 and δV⃗_0, it follows 

    δθ̇⃗̇≈ M_0^-1δV⃗ -M_0^-2δ M V⃗_0,

and we thus have 

    δθ̇⃗̇≤ M_0^-1δV⃗ + M_0^-1^2 V_0δ M.

Since we consider only the shot noise, δV⃗ and δ M can be described as

    δ M   ≈Δ_M/√(N_r), 
    δV⃗   ≈Δ_V⃗/√(N_r),

where N_r is the number of shots to evaluate M and V⃗.
When we denote the derivative of the state |Ψ(θ⃗)⟩ as

    ∂|Ψ(θ⃗)⟩/∂θ_k=∑_i=1^N_D,k g_k,i I_a ⊗(U_N⋯ U_k+1 P_k,i U_k-1⋯ U_1) |Ψ(θ⃗)⟩,

where g_k,i is a complex number and P_k,i is a Pauli operator (recalling that we have assumed the ansatz has the form of Eq. (<ref>) in the main text),
Δ_M can be written as

    Δ_M = 2 √(∑_k,q(∑_i,j |g_k,i^* g_q,j|^2)).
 
Similarly, Δ_V⃗ is written as

    Δ_V⃗=2√(∑_k (∑_i,j|g_k,i^* h_j|^2) ),
 
where the Hamiltonian is written as H=∑_j h_j P_j with P_j being a Pauli operator and h_j a coefficient.

In Ref. <cit.>, the implementation error is shown to be approximated as 

    δϵ_I(n)    ≡  D(|Ψ(θ⃗(n δ t))⟩, |Ψ^(0)(θ⃗(n δ t))⟩) 
       = √(δθ̇⃗̇^T B δθ̇⃗̇δ t^2 +O(δ t^3))
       ≲√(B )δθ̇⃗̇δ t  
       ≲√(B)Δ_I/√(N_r)δ t,

where Δ_I=M_0^-1Δ_V⃗+M_0^-1^2  V_0Δ_M and 

    B_i,j   = ∂⟨Ψ(θ⃗((n-1)δ t))|/∂θ_i∂|Ψ(θ⃗((n-1) δ t))⟩/∂θ_j
          - ∂⟨Ψ(θ⃗((n-1) δ t))|/∂θ_i|Ψ(θ⃗((n-1)δ t))⟩⟨Ψ(θ⃗((n-1) δ t))|∂|Ψ(θ⃗((n-1) δ t))⟩/∂θ_j.

Finally, the total implementation error will be

    ϵ_I ≡∑_n=1^N_stepδϵ_I(n) ≤√(B_max)Δ_I^(max)/√(N_r)T,

where the superscipt and subscript “max" on B and Δ_I  represents the maximum value for n=1,…,N_step.



 §.§ Resource estimation



  §.§.§ Evaluation of the error of the real time Green's function and the resource estimation

From Eqs.(<ref>)(<ref>)(<ref>), the error of the computed Green's function by our algorithm based on the VQS reads as

    ϵ^R    ≤α( 2(ϵ_s+ϵ_A+ϵ_I)+1/√(N_m)) 
       =α( 2 (ϵ_s+√(Δ^(2))T+√(Δ^(3)δ t)T +√(B)Δ_I/√(N_r)T ) + 1/√(N_m)),

where we omitted the subscript and superscript “max" for simplicity. The first and second terms in the right hand correspond to the state preparation error and the (part of) algorithmic error, respectively.
Both stems from imperfections of the ansatz to represent the ideal quantum state to be simulated and the analysis of them is not straightforward; e.g., it depends on details of the system (Hamiltonian) and/or the ansatz as well as the optimisation method used for the VQE.  We leave them for the future work.
On the other hand, the third, forth, and fifth terms in the right hand side can be suppressed by employing sufficiently small time step δ t and the large number of shots N_r, N_m for measuring the quantum circuits.
Those terms represent the algorithmic error due to a finite time step, the implementation error, and the shot noise in evaluating the Pauli matrices involved in the Green's function, respectively.

Let us evaluate the size of time step δ t and the number of shots N_r, N_m so as to bound each term in the right hand side of Eq. (<ref>).
To upper bound the third term by ε_A, we need the time step

    δ t= ε_A^2/4 α^2 Δ^(3)T^2.

It means that the number of time steps will be

    N_step=T/δ t=4α^2 Δ^(3) T^3/ε_A^2.

To upper bound the fourth term by ε_I, we need to set the number of shots to evaluate the matrix M as

    N_r= 4α^2 B T^2/ε_I^2.

Similarly, if we want to upper bound the fifth term by ε_m, it is required

    N_m = α^2/ε_m^2.


Finally, let us count the number of distinct runs of the quantum circuits to obtain the Green's function and discuss its dependence on the required error bound to the Green's function. 
The number of shots (distinct runs of quantum circuit) to populate M and V matrix is at most N_θ^2 N_D^2 + N_θ N_H N_D, where N_θ is the number of parameters (the number of elements of θ⃗), N_H is the number of terms of the Hamiltonian, and N_D is max_k N_D,k.
Hence the number of distinct runs of the quantum circuits to obtain the Green's function will be 

    N_tot   = N_step× N_r × (N_θ^2 N_D^2 + N_θ N_H N_D) + N_k × N_m 
       = 16 α^4 BΔ^(3)T^5/ε_A^2ε_I^2 (N_θ^2 N_D^2 + N_θ N_H N_D) + N_kα^2/ε_m^2,

where N_k is the number of the Pauli matrices in which we expand the electron creation operator c_k (see Eq. (<ref>) in the main text).
This is one of the main results for the resource estimation for calculating the Green's function based on the VQS, although this number is a very loose, or pessimistic estimation. 


We note that further simplification may be possible when we set ε_A=ε_I=ε_m=ε/3, i.e., we set the total error of the Green's function to ε, ignoring the errors originating from imperfections of the ansatz,

    N_tot   = 1296 α^4 BΔ^(3)T^5/ε^4 (N_θ^2 N_D^2 + N_θ N_H N_D) + 9N_kα^4/ε^2
       ≈1296 α^4 BΔ^(3)T^5/ε^4 (N_θ^2 N_D^2 + N_θ N_H N_D).

We will use this relationship in the next subsubsection.



  §.§.§ Evaluation of the error of the Green's function in frequency domain and the resource estimation

The Green's function in frequency domain is defined as

    G̃^R_k(ω)= ∫_0^∞ e^i (ω+iη) t G_k^R(t) dt  (η→ +0),

and its imaginary part Im G̃^R_k(ω) gives the spectrum function.
Suppose the Hamiltonian of interest has the energy spectrum { E_min, …, E_max}, where E_min(max) is the minimal (maximum) energy. We assume all energies are positive, E_min>0, without loss of generality.
According to the sampling theorem, to avoid aliasing in G̃^R_k(ω), we need to sample the real time Green's function G_k^R(t) with over twice the maximum frequency contained in it.
It means that we need to have the time step 

    δ t ≤1/2 f_max= 1/E_max/ π,

where f_max = E_max/2π is the maximum frequency of the real time Green's function, in the VQS algorithm.
It is important to note that this bound is constant with regard to the required accuracy to the real time Green's function while the time step shown in Eq. (<ref>) depends on it,
so we adopt Eq. (<ref>) as the time step to be considered. 
Meanwhile, the necessary evolution time T_0 to obtain information of the lowest frequency component is T_0=2π/E_min.

In our simple numerical integration scheme, the real time Green's function is sampled at t=0, δ t, …, T_0.
The Green's function in frequency domain is approximated as  

    G̃^R_k(ω)    ≈∫_0^T_0 e^i (ω+iη) t G_k^R(t) dt  (η→ +0) 
       ≈∑_n=0^N_step G_k^R(n δ t) e^i ω n δ tδ t + O(δ t)

where N_step = T_0/δ t is the total number of time steps for simulation.
We note that the first term in the last line is O(δ t^0).
When the errors of all data { G_k^R(n δ t) }_n=0^N_step are bounded by ε̃, the error of the imaginary part of G̃^R_k(ω) can be upper bounded as

    δ( Im[G̃_k^R(ω)] )
     ≲ N_step·ε̃·δ t
     = ε̃  T_0.

The same argument can be made for the real part.
Therefore, if we want to set total accuracy of G̃^(R)_k(ω) to ε, i.e., we should take ε̃ = ε / T_0.
In that case, the total number of measurements (Eq. (<ref>)) will be 

    N_tot   = 1296 α^4 BΔ^(3)T_0^9/ε^4(N_θ^2 N_D^2 + N_θ N_H N_D) 
       = (2π)^9 × 1296  α^4 BΔ^(3)/ε^4 E_min^9(N_θ^2 N_D^2 + N_θ N_H N_D).



Conversely, when the time step δ t is fixed and take only errors related to δ t into account, the error of the Green's function in frequency domain is

    δ( Im[G̃_k^R(ω)] )
     ∼ 2α T_0^3/2√(Δ^(3)δ t).




§ DETAILS OF DISCUSSION ON THE FEASIBILITY OF OUR ALGORITHMS IN SEC.<REF>



In this section, we provide details of the discussion on the feasibility of our proposed algorithms in Sec. <ref> based on Ref. <cit.>.

We take the the two-dimensional Hubbard model (Eq. (<ref>) in the main text) of N_site sites with the open boundary condition as an example.
By applying the Jordan-Wigner transformation, the Hamiltonian will have the form of H=∑_j h_j P_j, where h_j is real coefficient and P_j is Pauli matrices.
The Hamiltonian ansatz is defined as (Eq. (<ref>) in the main text)

    ∏_d=1^n_d( ∏_j exp( i θ_j^(d) P_j ) ),

where n_d is the depth of the ansatz and θ_j^(d) is a parameter to be optimized in the ansatz.
In Ref. <cit.>, it is shown that the number of single qubit gates and two-qubit gates to implement the single-depth (n_d=1) Hamiltonian ansatz of this model is

    N_single   =4 N_site^3/2+7  N_site-4 √(N_site), 
    
    N_two   =8 N_site^3/2+ N_site-4 √(N_site),

with assuming the Z rotation gate and the partial swap gate as elementary single and two-qubit operations.
Therefore, the total number of gates is n_d(N_single+N_two).
By putting n_d = 1 and N_site=25 as in the main text, we obtain N_single≈ 650, N_two≈ 1000.
Alternatively, when we set n_d=N_site, the total number of gates is N_single≈ 16400, N_two≈  2600.
In this case, the tolerable error rate for quantum error correction, N_gateϵ_gate≲ 2, gives ϵ_gate = 8 × 10^-2% even if we consider only the two-qubit gates.
This value indicates that the further  improvement of two-qubit gate fidelity is required.



§ ADDITIONAL SIMULATION RESULTS


Here, we show additional numerical results. The results for n_d=4 are presented in Fig. <ref>. We can see that the results are almost identical with the result of n_d=8 in Fig. <ref> in the main text although the peaks of the spectral function are slightly smeared.

Furthermore, we present results of the variational quantum simulation (VQS) for four-site Fermi-Hubbard model for further convincing the readers of the feasibility of our algorithm.
The Hamiltonian of the four-site Fermi-Hubbard model is defined as

    H = - ∑_i,=1,2,3,4,σ=↑,↓( c_i,σ^† c_i+1,σ + h.c.)
     + U ∑_i=1^4 c_i,↑^† c_i,↑ c_i,↓^† c_i,↓ -U/2∑_i=1,2,3,4,σ=↑,↓ c_i,σ^† c_i,σ.

Simulating this model requires eight qubits, so it is closer to the current experiments of the near-term quantum computers.
We employ the Hamiltonian ansatz (<ref>) of the depth n_d=16 which has 320 parameters (the number of Pauli terms in the Hamiltonian except for the identity operator is N_P=20).
The result of the numerical simulation of the VQS is displayed in Fig. <ref>, which well reproduces the exact result.
While general unitary operators on the system have (2^8)^2=65536 real parameters, our ansatz only uses real 320 parameters to simulate the dynamics for computing the Green's function by using the VQS algorithm in Sec. <ref>. 








§ DETAILS ON NUMERICAL SIMULATION  

For the numerical calculation in Sec. <ref>, we evaluate the derivatives appearing in Eq. (<ref>) by numerical differentiation by taking δθ_i=10^-4 for all parameters θ_i of the ansatz quantum state.
We solve Eq. (<ref>) by minimizing Mθ⃗̇⃗ - V⃗ by using  function implemented in Python library Numpy <cit.> with neglecting the singular values of M smaller than 10^-8.
The size of the time step dt is taken so small that there is no significant change in the final result when decreasing dt. 


