Calculation of the Green's function on near-term quantum computers

The Green's function plays a crucial role to study the nature of quantum many-body systems, especially strongly-correlated systems. Although the development of quantum computers in near term is expected to enable us to compute energy spectrum and energy eigenstates of such classically-intractable systems, the 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 the 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 Hubbard model and demonstrated the validity of our proposals.


I. INTRODUCTION
The advent of a primitive but still powerful form of quantum computers, called noisy intermediate-scale quantum (NISQ) devices, is approaching [1].NISQ devices have a few hundreds to thousands of qubits under highly precise control although they are not faulttolerant.It is believed that the behaviour of NISQ devices will soon reach a stage that one can not simulate its dynamics by classical computers due to exponentiallyincreasing size of the Hilbert space with the number of qubits [2][3][4][5][6].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 [1].
One of the most promising tasks which near-term quantum computers will 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, material science, pharmacy, etc [7,8].The most celebrated algorithm of the quantum simulations for near-term quantum computers is variational quantum eigensolver (VQE) [9][10][11][12] in which energy eigenstates and eigenenergies of the system are obtained based on the variational principle of the quantum mechanics.Although the VQE is originally proposed for finding only the ground state, its extension to excited states was discussed in several literatures [13][14][15][16][17].
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: the Green's function and the spectral function [18][19][20].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 excitation of a system, which gives crucial information on the high-T c superconductivity [21], magnetic materials [22], and topological insulators [23].While several methods based on the Suzuki-Trotter decomposition of the time evolution operator [24,25] or the quantum phase estimation [26,27] were already proposed to calculate the Green's function on quantum computers, they require a large number of qubits and gate operations which are hard to be realized 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 [28][29][30][31][32] for an efficient calculation of the Green's function in real time.We extend the original VQS algorithm in order to calculate 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 functions [18][19][20].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) [13] or the multistate contracted VQE (MCVQE) [16], one can compute the spectral function (hence the Green's function).We confirm the validity of our methods by numerical simulations for the Hubbard model, a prototype of strongly-correlated systems.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.II, we briefly review the definition of the Green's function and the spectral function at zero temperature.In Sec.III, 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 [13] and the MCVQE algorithm [16] in Sec.IV.Section V is dedicated to demonstration of our methods by performing numerical simulations calculating the spectral functions of the 2-site Hubbard model.We present the extension of our methods to the finite temperature Green's function in Sec.VI.We discuss implications and possible future directions of our results, and make a conclusion in Sec.VII.Appendix A gives details of the numerical simulations.

II. 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 [18][19][20].
Let us consider a fermionic system described by Hamiltonian H which is composed of fermionic creation (annihilation) operators {c a , c † a } a , where a is a label to specify the fermionic mode.For example, a can be (k, σ), where k denotes the momentum and σ =↑, ↓ does the spin of the fermion.The retarded Green's function at zero temperature is defined as 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 constant k B are = k B = 1.For simplicity, throughout this paper we consider the Green's function G R ab with a = b = (k, ↑), namely, the Green's function in the momentum space with identical spin.We simply write 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 function where η → +0 is a factor to assure the convergence of the integral.The spectral function and the Green's function GR k (ω) have a relation Finally, we introduce the Lehmann representation of the spectral function which utilizes the energy eigenvalue E n and the eigenstate of the Hamiltonian |E n , where we call the first (second) term as particle (hole) part of the spectrum function.
In Sec.III, we compute the spectral function by performing the Fourier transformation to G R k (t) calculated by the VQS-based method while we compute it by the Lehmann representation (6) with quantities calculated by the VQE-based method in Sec.IV.

III. COMPUTATION OF GREEN'S FUNCTION WITH VARIATIONAL QUANTUM SIMULATION
In this section, we first review the VQS algorithm [28][29][30][31] 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.

A. Review of variational quantum simulation
Here we review the variational quantum real time simulation algorithm introduced in Ref. [28].Let us consider an ansatz quantum state created by a parametrized quantum circuit, where U 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 that U i (θ i ) is composed of a set of Pauli rotation gates, e iα (j) θiP (j) with a coefficient α (j) and a Pauli matrix P (j) , and other non-parametrized gates.For a given initial state |ψ( θ(0)) and Hamiltonian H, the VQS algorithm finds the solution of the Schrödinger equation, within the Hilbert space spanned by the ansatz quantum state, {|ψ( θ) } θ .Specifically, the time evolution described by Eq. ( 8) is mapped to the time evolution of parameters θ(t).Although there are several variational principles to map Eq. ( 8) to the equations for θ(t) [33][34][35], we simply choose the McLachlan's variational principle [35] in this paper because it is the most stable and physically reasonable among them [29].The McLachlan's variational principle [35] maps Eq. ( 8) to the equation of motion of the parameters θ(t) by minimizing the distance between the exact evolution of the Schorödinger equation and the evolution of the parametrized ansatz state under infinitesimal variation of time δt [28], where |ϕ = ϕ|ϕ is the norm of |ϕ .One can explicitly write down the equation determining where 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.[28,36].Finally, from the solution of Eq. ( 10), 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 holds for a different initial state |ψ 1 in general.

B. 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, 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 [9][10][11][12] or the variational imaginary time simulation algorithm [30,37].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 [7,8,[38][39][40], where P n is a tensor product of Pauli matrices acting on N qubits that satisfies is a (complex) coefficient.For example, one can adopt the Jordan-Wigner transformation [40] for the decomposition.The first term of Eq. ( 12) can be rewritten as Therefore, the problem reduces to finding a way to compute G| e iHt P i e −iHt 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 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 Then, G| e iHt P i e −iHt P j |G can be evaluated as the transition amplitude of P i between e −iHt |G and e −iHt P j |G , This can be evaluated by using the quantum circuit shown in Fig. 1.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 FIG. 1.Quantum circuit to compute Eq. ( 17).The upper horizontal line represents the ancillary qubit and the lower line does the qubits for the system of interest.The initial state for the ancillary qubit is taken as 2 −1/2 (|0 + e iφ |1 ).The expectation value of Z-measurement on the ancillary qubit yields Re e iφ G| U (1) ( θ1) † PiU (2) ( θ2)Pj |G /2.Hence by choosing φ = 0, π/2, we can measure both the real and imaginary part of Eq. ( 17).
2. Quantum circuit to compute Eq. ( 18).Again, the upper horizontal line represents the ancillary qubit and the lower line does the qubits for the system of interest.The initial state for the ancillary qubit is taken as 2 −1/2 (|0 + e iφ |1 ).The expectation value of Z-measurement on the ancillary qubit yields Re e iφ G| U ( θ) † PiU ( θ)Pj |G /2, so we can measure both the real and imaginary part of Eq. ( 18) by choosing φ = 0, π/2.The number of required controlled operations is only two in this case.We can further eliminate controlled operations by using the method proposed in Ref. [36].
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 roots 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 In this case, the quantum circuit for evaluating the term is significantly simplified as depicted in Fig. 2. In the followings, we 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−1 l=0 .Let us consider a state with ancilla, 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−1 l=0 .To be more concrete, the ansatz state for the VQS algorithm is, where we define |ψ l ( θ) s = U s ( θ) |ψ l s .The M matrix and V vector in Eq. ( 11) become and 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. ( 21), ( 22) 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. 2. 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 (20) for the VQS algorithm can be constructed by a quantum circuit shown in Fig. 3.
To sum up, the calculation of each term in Eq. ( 14), 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 nearterm quantum computers, such as the VQE [9][10][11][12].
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. 3 will find such U ( θ).
2. Evaluate G| U ( θ) † P i U ( θ)P j |G by the quantum circuit shown in Fig. 2.

IV. 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 [13] and the MCVQE method [16], and take advantage of the Lehmann representation of the spectral function (6).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.[13,16], so our contribution will be application of it for calculation of the Green's function.
To compute the Lehmann representation of the spectral function (6), we also need the transition amplitude of the fermions c k,↑ , c † k,↑ , such as Ẽj |c k,↑ | Ẽ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 | Ẽj 's are created from the same unitary gate U ( θ * ).Specifically, if we can easily make superpositions of the input states, |ψ j and |ψ k , Ẽj |c k,↑ | Ẽ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. ( 14) and decompose it into real part and imaginary part, c k↑ = A k + iB k , where A k and B k are Hermite operators.Then we have (24) Each term can be evaluated by using |ψ 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 superposition of them on quantum computers is not so difficult.Therefore, by substituting eigenenergies of Eq. ( 6) with Ẽj and the transition amplitudes with Ẽ0 |c k,↑ | Ẽk , etc., we can evaluate the spectral function, and the Green's function accordingly.

B. 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 [17,41], and is essentially the same as the MCVQE algorithm [16].
The procedure of the algorithm is as follows.First, we prepare orthonormal input states {|ψ j K−1 j=0 }, which are simple and easy to realize a superposition of them on quantum computers.Then we minimise the cost function with respect to parameters θ, where U ( θ) is the ansatz quantum circuit.
After the optimisation, the subspace spanned by {U ( θ * ) |ψ j } K−1 j=0 will be close to that spanned by the true K eigenstates {|E j } K−1 j=0 , where θ * is the parameters after optimisation.At this stage, U ( θ * ) |ψ j is generally the superposition of the excited states {|E j } K−1 j=0 and ψ j | U ( θ * ) † HU ( θ * ) |ψ j is not a good approximation to the true eigenvalue E j .Therefore, to obtain the eigenstates and eigenvalues of H, we solve the eigenvalue problem within the subspace spanned by {U ( θ * ) |ψ j } K−1 j=0 , where and the approximate eigenenergies appear as Ẽ j = E j,j .The transition amplitude 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 matrix C (k) m,n , and evaluate the spectral function by Eq. ( 6) and the Green's function.

V. 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 Hubbard model of two sites, where t is a parameter characterizing hopping between the sites and U denotes the strength of the on-site Coulomb repulsion [42][43][44].We simulate two proposed protocols in Secs.III and IV by classical computers with the fast quantum circuit simulation library Qulacs [45].
We exploit the Jordan-Wigner transformation [40] to map the fermionic Hamiltonian (30) into the qubit one with four qubits by using the library OpenFermion [46].

A. Numerical simulation of the method based on variational quantum simulation
We calculate the real-time Green's function of the model (30) at zero temperature by using the method described in Sec.III.First, we prepare the ground state of the model ( 30) by the standard VQE algorithm with a hardware-efficient-type ansatz [10] depicted in Fig. 4.Then, we perform the VQS algorithm.As an ansatz quantum state, we adopt the so-called variational Hamiltonian ansatz [47,48] inspired by the Suzuki-Trotter decomposition of the time evolution operator e −iHt .The variational Hamiltonian ansatz is defined through the qubit representation of the Hamiltonian, where P m is a (multi-qubit) Pauli matrix and c m is a coefficient.An ansatz state for the variational Hamiltonian  30) of U = 3 (top row) and U = 6 (bottom row).The time step is taken as dt = 0.1(0.003)for U = 3 (6).The exact spectral function is calculated by the exact dynamics of the Green's function in real time from t = 0 to t = 100 with step dt = 0.1.We take η = 0.2 for the calculation of the spectral functions.
ansatz is given by |ψ(θ) = U V HA (θ) |ϕ 0 , where and n d denote the depth of the ansatz.If the qubit representation of the Hamiltonian (31) has N P terms, the number of parameters of the ansatz will be N P n d .In our simulation, N P = 8 and we choose n d = 8, so there is 64 parameters in the parametrized quantum circuit whereas general unitary operators on the system have (2 4 ) 2 = 256 parameters.Further details on numerical calculations are described in Appendix A. The result is shown in Fig. 5.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.Next, we numerically simulate the method described in Sec.IV.We adopt the symmetry-preserving ansatz [49] drawn in Fig. 6, 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 (30).As an input states, we simply choose the computational basis states with desired particle number:

VI. 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 where |E n and E n denote the eigenstates and eigenvalues of Hamiltonian H, respectively.The corresponding spectral function is A. Variational quantum simulation for Green's function at finite temperature Equation ( 33) 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 FIG. 8.The ansatz quantum circuit for obtaining the purified Gibbs state |Φ(β) for N = 2. finite-temperature version of the density matrix renormalization group method [50,51].We note that several methods based on the typically of the Hilbert space to evaluate physical quantities at finite temperature [52][53][54] 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 |i s,e is the computational basis of the system and environment.By using the method of variational imaginary time evolution introduced in Ref. [29], 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. 8 will produce |Φ(β) .We note that |Φ(β) satisfies 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 where superscript of the Pauli operator P (s) i implies that it only acts on the system s.In Fig. 9, we show the ansatz quantum circuit to construct the unitary gate U s ( θ(t)).By substiting the above quantity into Eq.( 33), we can evaluate (each term of) G k (t; β).
We finally remark on another way to evaluate Eq. ( 33) based on the VQS algorithm.It is possible to obtain several approximate eigenenergies { Ẽn } K n=1 and eigenstates {| Ẽn } K n=1 of the system by the SSVQE or MCVQE algorithm, and to perform the VQS algorithm in Sec.III for each obtained eigenstate | Ẽn .The result approximates G R kn (t) in Eq. ( 33), so substituting it as well as Ẽn into Eq.(33) 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 β.

B. Computing Green's function at finite temperature with excited-states search method
Extension of the algorithm in Sec.IV to the finite temperature Green's function is rather simple.Since we already introduce the way to evaluate the quantity | Ẽn |c k,↑ | Ẽm | 2 , putting them into Eq.(34) 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.

VII. 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 or imaginary time simulation algorithm 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 the method for constructing the variational quantum circuit which acts as the time evolution operator for multiple initial states simultaneously, and make 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 [55] expressed as φ BA (t − t ) = i B(t − t ), A(0) 0 , B(t − t ) = e iH(t−t ) Be −iH(t−t ) (37) 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 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 to 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, materials science, etc.

FIG. 3 .
FIG.3.The ansatz quantum circuit for the VQS algorithm to construct the variational unitary gate Us( θ(t)) which approximates the time evolution operator U (t) = e −iHt for |G and Pj |G .Here, we assume that the ground state |G is obtained as |G = UG |ϕ0 .

FIG. 4 .
FIG. 4. Hardware-efficient-type ansatz used to obtain the ground state of the model (30) for the VQS algorithm.Each rotational gate RY (θ) = e iθY /2 , RZ (θ ) = e iθ Z/2 has a parameter angle and D denotes the depth of the ansatz.

ImG
FIG.5.Numerical simulation of the VQS algorithm to compute the Green's function in real time G R k (t) (left and center columns) and the spectral function (right column) for the model (30) of U = 3 (top row) and U = 6 (bottom row).The time step is taken as dt = 0.1(0.003)for U = 3(6).The exact spectral function is calculated by the exact dynamics of the Green's function in real time from t = 0 to t = 100 with step dt = 0.1.We take η = 0.2 for the calculation of the spectral functions.

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

FIG. 7 .
FIG. 7. Result of numerical simulation of the method described in Sec.IV for calculating the spectral function of the model (30) at U = 3 (left) and U = 6 (right).

ACKNOWLEDGMENTSA
part of the numerical simulations in this work was done on Microsoft Azure Virtual Machines provided through the program Microsoft for Startups.YON and SE acknowledge valuable discussion with Kosuke Mitarai and Xiao Yuan.IK was supported by Qunasys Inc.This work was also supported by MEXT Q-LEAP.