Hybrid quantum-classical algorithm for computing imaginary-time correlation functions

Quantitative descriptions of strongly correlated materials pose a considerable challenge in condensed matter physics and chemistry. A promising approach to address this problem is quantum embedding methods. In particular, the dynamical mean-field theory (DMFT) maps the original system to an effective quantum impurity model comprising correlated orbitals embedded in an electron bath. The biggest bottleneck in DMFT calculations is numerically solving the quantum impurity model, i.e., computing Green's function. Past studies have proposed theoretical methods to compute Green's function of a quantum impurity model in polynomial time using a quantum computer. So far, however, efficient methods for computing the imaginary-time Green's functions have not been established despite the advantages of the imaginary-time formulation. We propose a quantum-classical hybrid algorithm for computing imaginary-time Green's functions on quantum devices with limited hardware resources by applying the variational quantum simulation. Using a quantum circuit simulator, we verified this algorithm by computing Green's functions for a dimer model as well as a four-site impurity model obtained by DMFT calculations of the single-band Hubbard model, although our method can be applied to general imaginary-time correlation functions.


I. INTRODUCTION
The accurate and efficient computation of quantum many-body systems is critical in computational physics.As Feynman pointed out [1], quantum computers, which use quantum mechanics as their computational principle, are expected to efficiently simulate quantum systems.There are obvious advantages in this approach in the field of strongly correlated electron such as frustrated magnetic compounds and high-temperature superconductors, which are difficult to simulate using classical computers [2][3][4].However, even if a quantum computer with thousands or tens of thousands of logical qubits is realized, simulating macroscopic degrees of freedom in solids remains challenging.
We need a framework to reduce the degrees of freedom of solids to make them realistically computable by quantum computers.One promising approach is the combination of quantum computing and quantum embedding theories such as the dynamical mean-field theory (DMFT) [5,6].As illustrated in Fig. 1, DMFT maps the original system to an effective quantum impurity model comprising correlated orbitals embedded in an electron bath.The biggest bottleneck in DMFT calculations is solving the quantum impurity model numerically, that is, computing the Green's function.Although advanced classical algorithms such as the tensor network [7][8][9] and quantum Monte Carlo [10] have been developed, their applications are limited to a few correlated orbitals owing to the exponential growth of the computational cost with respect to the number of correlated orbitals.This difficulty * sakurairihito@gmail.com originates from the notorious negative sign problems and the rapid growth of quantum entanglement entropy.
It is highly desirable to solve the quantum impurity model with a quantum computer and overcome the above-mentioned difficulties.In recent years, theoretical proposals have been made to compute the Green's function in polynomial time using a quantum computer [11][12][13].However, these proposals assume a fault-tolerant quantum computer, which is beyond the capabilities of quantum devices with limited hardware resources such as noisy intermediate-scale quantum (NISQ) devices [14].Therefore, the development of methods to calculate the Green's function using quantum devices with limited hardware resources has been actively pursued [15][16][17][18].
In quantum embedding simulations, performing calculations in imaginary time rather than in real time is favorable for quantum computation with quantum devices in terms of the required quantum hardware resources.This is because the imaginary-time formalism allows us to discretize the environment with fewer auxiliary bath sites [19][20][21] and hence qubits.Although recently proposed algorithms has been used to compute the imaginary-frequency Green's function [17,18,22], it requires the evaluation of the expectation value of the square of the Hamiltonian or an effective Hamiltonian, which may be expensive for a large impurity model.Thus, efficient methods for computing the imaginarytime Green's function using quantum devices with limited hardware resources still remain to be explored.
In this study, we propose a quantum-classical hybrid algorithm to compute the imaginary-time Green's function on quantum devices with limited hardware resource devices.We apply a variational quantum simulation (VQS) [23], which is a variational algorithm for time evolution, to the computation of the imaginary-time FIG.1: Schematic of dynamical mean-field calculations and quantum impurity models.
Green's function.Our algorithm efficiently computes the imaginary-time Green's function by adopting a nonuniform mesh focusing on the shape of the imaginarytime Green's function.Using a quantum circuit simulator, we verified our method for typical impurity models, such as a dimer model and a four-site impurity model obtained by DMFT calculations of the single-band Hubbard model.
The remainder of the paper is structured as follows.In Sec.II, we review the finite-temperature and zerotemperature formalism of Green's functions.In Sec.III, we propose a variational quantum algorithm for computing the imaginary-time Green's function.In Sec.IV, we numerically demonstrate the proposed algorithm for typical impurity models using a quantum circuit simulator.Finally, we discuss the scalability and numerical stability of our algorithm, and compare it with other methods.

A. Finite-temperature formalism
We consider a fermionic system in the grand-canonical ensemble and denote its Hamiltonian as H.
where c i /c † i denote the creation and annihilation operators of the spin orbital i, and N is the number of spin orbitals.t ij , U ikjl , and µ denote the hopping matrix, Coulomb tensor, and chemical potential, respectively.The retarded Green's function is defined as where ĉa (t) = e iHt ĉa e −iHt and ĉ † b (t) = e iHt ĉ † b e −iHt denote the annihilation and creation operators for spin orbitals a and b in the Heisenberg representation, respectively.θ(t) is a step function.Throughout this paper, we take = k B = 1.The thermal expectation • • • is evaluated in the ground-canonical ensemble.The retarded Green's function can be transformed into the (real) frequency space as where ω is a real frequency.On the contrary, the imaginary-time Green's function is defined as where ĉa (τ ) = e τ H ĉa e −τ H . Notably, the imaginarytime Green's function is anti-periodic, as G ab (τ + β) = −G ab (τ ).The Fourier transform of the imaginary-time Green's function (Matsubara Green's function) is given by where ω = (2n + 1)π/β, where n ∈ N and β = 1/T .The imaginary-frequency Green's function G(iω) can be analytically continued from the imaginary axis to the full complex plane as G ab (z).The analytically continued G ab (z) has the spectrum representation with where z is a complex number and n, m runs over all eigenstates of the system (E m/n denotes an eigenvalue).
There are poles for a finite system or a branch cut for an infinite system on the real axis.The retarded Green's function is G ab (z)'s value just above the real axis: B. Zero-temperature formalism We now consider the limit of T → 0, where the ensemble average can be restricted to the ground state(s) Ψ G .At sufficiently low temperatures, for 0 < τ < β/2, Eq. ( 4) can be rewritten as where If the ground-state manifold is degenerate, then Eq. ( 9) must be averaged over all degenerate ground states.In general, |G ab (τ )| decays exponentially for an insulating system, whereas its decay is algebraic for a metallic system.We must increase β, which sets the upper limit of the time evolution such that G ab (τ ) is sufficiently small at the boundary.Similarly, for −β/2 < τ < 0, we obtain Equations ( 9) and ( 10) can be represented in a unified form as follows: where A + = ĉa and B + = ĉ † b for τ > 0, and A − = ĉ † b and B − = ĉa for τ < 0. The signs ∓ are for τ > 0 and τ < 0, respectively.Once G ab (τ ) is evaluated on a sufficiently fine mesh in [−β/2, β/2] for a sufficiently large fictitious inverse temperature β, we can transform the data to the imaginary-frequency space, where the DMFT calculations are usually implemented, by numerically evaluating the integral

A. Overview
We propose a variational algorithm for computing the imaginary-time Green's function.First, we choose β and introduce a fine mesh in [−β/2, β/2].Subsequently, we compute the Green's function for the positive and negative sides of τ with Eq. ( 11) using the following procedure: Stage 1: The ground state is computed by a variational quantum algorithm.

Stage 2:
The operator B± is applied to the ground state.
Stage 3: Imaginary-time evolution is performed by VQS.
Stage 4: The Green's function is evaluated by computing the transition amplitude for Â± .
Figure 2 illustrates the entire procedure.In the following subsections, we explain the stepwise details of the procedure.This method can be applied to general two-point correlation functions in the form of Ψ G | Â(τ ) B(0)|Ψ G , where Â/ B is an equal-time bosonic or fermionic operator.

B. Mapping to qubits
To calculate a fermionic system on a quantum computer, it is necessary to convert the fermionic operators from a fermionic second quantized representation to a qubit representation in advance.Typical methods include the Jordan-Wigner transformation [24] and Braviy-Kitaev transformation [25,26].
As an example, we consider the Hamiltonian in the form of Eq. (1).The Hamiltonian can be transformed to the qubit representation as where S k ∈ {X, Y, Z, I} ⊗m are the tensor products of Pauli operators with m qubits, which are transformed from the terms in Eq. ( 1); h p are the coefficients.In this study, we adopt the Jordan-Wigner transformation given by C. Stage 1: Ground-state calculations We use the variational quantum eigensolver (VQE) [27] to compute the ground state.VQE is a variational quantum algorithm for determining the ground state and its energy of the Hamiltonian in Eq. ( 13) using a quantum computer.
The flow of the algorithm is described as follows.
Step 1: Prepare an initial state |Ψ init (usually an unentangled product state) on a quantum computer.
Step 3: Measure S k in the Hamiltonian of Eq. ( 13) using a quantum computer.
Step 4: Calculate H on a classical computer.
Step 5: Update the parameters on the classical computer to reduce H .
By repeating steps 2-5, we obtain a set of converged parameters θ * .If the representation capability of the ansatz is sufficiently high and an appropriate initial guess is used, the optimized variational quantum state Ψ( θ * )) should approximate the ground state |Ψ G accurately.Throughout this paper, we assume that the quantum circuit conserves the number of electrons.The number of electrons in the Hilbert space to be searched can be FIG.2: Overview of the methods for computing the imaginary-time Green's function.
fixed by the number of electrons in the initial state and a number-conserving circuit.
In step 5, if we update the parameters using the gradient method, evaluating the derivative of H is necessary.This can be done on a quantum computer either by using numerical finite differentiation or parameter-shift rules [28].

D. Stage 2: Single-particle excitation
In stage 2, we compute a variational quantum state for the single-particle excited state B± |Ψ G .Because the operator is not unitary, we represent the resultant state as where c 1 is a coefficient and the parametrized quantum state φ EX ( θ EX ) is defined by where the initial state |φ EX must have N ± 1 electrons because the operator changes the number of electrons by ±1.
We compute the variational parameters θ EX and coefficient c 1 as follows: Step 1: Transform the creation operator B± to the qubit representation, and prepare an initial guess for θ EX .
Step 2: Prepare a variational quantum state φ EX ( θ EX ) and measure the cost function on a quantum computer.We repeatedly optimize the parameters to reduce the cost function until the parameters converge.
Step 3: For the converged parameters θ * EX , measure In steps 2 and 3, we evaluate the cost function and c 1 on a quantum circuit.We explain the method of evaluating c 1 in step 3 because the cost function in step 2 can be directly obtained by squaring the absolute value of c 1 [29].
First, we decompose B ± as the sum of its Hermitian part and its anti-Hermitian part using the Jordan-Wigner transformation, as in Eq. ( 14), and Eq. ( 15).We evaluate the transition amplitude on a quantum computer by measuring the Hermitian and anti-Hermitian parts of the following form: where P are Pauli operators with m qubits, and U 1 and U 2 are unitary operators with m qubits.Equation ( 18) can be measured using the quantum circuit in Fig.
3 [15,17], which requires one ancilla qubit.Let p 0 /p 1 be the probability of measuring 0/1 in the ancilla qubit.The real and imaginary parts of the transition amplitude can be measured separately by setting φ = 0 and π/2 in the R z gate, respectively, as As this method is based on a single ancilla qubit, we need complex quantum circuits owing to the use of the control unitary operators .
If we use a derivative-based optimization algorithm in step 2, we additionally need to measure the partial derivatives of C with respect to θ EX .The analytical form of the derivative of this cost function can be calculated using the quantum circuit of the same form to evaluate the transition amplitude.
In the following, we use the same quantum circuit to evaluate the quantities in the same form as Eq. ( 18) in stages 3 and 4.

E. Stage 3: Imaginary-time evolution
In stage 3, we calculate the imaginary-time evolution of φ EX ( θ * EX ) obtained in stage 2. As the imaginary-time FIG.3: Quantum circuit for calculating Eq. ( 18).This quantum circuit uses m qubits (bottom line) and one ancilla qubit (top line).The transition amplitude can be calculated by combining the Z measurement results for the ancilla qubit.
evolution is not unitary, we represent the imaginary-time evolved state [30] as where η(τ ) is a τ -dependent real number.A similar parametrization for nonunitary time evolution was proposed for the application of VQS to financial systems [31].Also, a different normalization factor calculation method was proposed in the context of calculation of the Gibbs partition function, which was the first paper to calculate the normalization factor [32].We perform the imaginary-time evolution on a (generally non-uniform) mesh, {τ 1 , τ 2 , • • • , τ N } (τ 1 = 0 and τ N = β/2).This can be achieved by determining the variational parameters θ EX (τ t ) and η(τ t ) on the mesh points sequentially from τ = 0 to τ = β/2.

Time-dependent variational principle
We first review the time-dependent variational principle (TDVP), on which VQS is based.The timedependent Schrödinger equation reads where τ is an imaginary time.The imaginary-time evolution of the normalized ket obeys where Ψ(τ ) We now parametrize these two kets as where φ|φ = 1, the vector θ(τ ) denotes the τ -dependent real variational parameters, and η(τ ) is a real parameter for the norm.In the TDVP, the time evolution of Eq. ( 22) is mapped to the time evolution of θ(τ ).In McLachlan's variational principle, we minimize the distance between the exact evolution and the evolution of the parametrized state under infinitesimal variation of the imaginary time δτ as where • • • denotes the Frobenius norm.There are several possible ways to solve this equation.We will explain them later .
Once the τ dependence of θ is determined, we can compute η(τ ) by solving the differential equation with an appropriate initial condition at τ = 0. Equation ( 26) can be derived by substituting Eq. ( 24) into Eq.( 21).

VQS
We explain how to perform the imaginary-time evolution of a quantum state on a discrete mesh in τ in VQS.Here, we define φ( θ(τ ) as U ( θ) |φ on the quantum circuit, where U ( θ(τ )) is a unitary operator with τ -dependent real parameters and |φ is an initial state.One can explicitly write the equation for determining the time derivative of the variational parameters at τ : where and N P denotes the number of variational parameters.
Equations.( 28) and ( 29) involve quantum circuits differentiated with respect to the variational parameters.As discussed in detail in Ref. 32, the matrix elements of M and C can be efficiently computed on a quantum computer using one ancilla qubit, for example, by differentiating the variational quantum circuits explicitly (see Appendix A for more details).
The linear system in Eq. ( 27) can be solved efficiently on a classical computer [33].Subsequently, we evolve the quantum state from τ to τ + ∆τ (∆τ > 0) by updating the variational parameters as Direct VQS Here, we propose an alternative way to perform the imaginary-time evolution, which is based on the direct minimization of Eq. ( 25) for a finite time step, ∆τ .We name this approach "direct VQS".The optimization problem in Eq. ( 25) is equivalent to the following optimization problem: where Ψ(τ ) = φ( θ(τ ) .The terms that do not depend on θ(τ + ∆τ ) are excluded from the cost function.Equation (31) becomes exact in the limit ∆τ → 0. This optimization problem can be efficiently solved using θ(τ ) as an initial guess.We can evaluate the quantities in Eq. ( 32) by using a quantum circuit as Eq. ( 18), after decomposing the Hamiltonian into a sum of terms in the form of Eq. ( 13).

Computational complexity of VQS and direct VQS
The computational complexity of VQS and direct VQS depends on the Hamiltonian and the ansatz.We discuss the computational complexity of VQS and direct VQS for a Hamiltonian with a general two-body interaction and a quantum impurity model.In this discussion, we assume the unitary coupled cluster ansatz with generalized singles and doubles (UCCGSD) [34,35] (see Appendix B for more details).For this ansatz, the number of variational parameters N P and the gate depth scale as O(n 4 so ) [n so is the total number of spin orbitals].After Jordan-Wigner transformation, the gate depth scales as O(n 5 so ).However, the gate depth reduces to O(n 4 so ) if we can maximally parallelize the terms in the unitary coupled cluster operator on a near-term quantum computer [36].

A Hamiltonian with a general two-body interaction
In the case of a Hamiltonian with a general two-body interaction, the number of terms N h in the Hamiltonian scales as O(n 4 so ).
[VQS] The bottleneck of imaginary-time evolution using VQS is the evaluation of the vector C and the matrix M on a quantum computer.The computational complexity for measuring the elements of C scales approximately as O(N depth N P N h ), where N depth is the depth of a quantum circuit.In particular, for the UCCGSD ansatz, O(N ]. On the other hand, approximately, the computational complexity for measuring the elements of M scales as O(N depth N 2 P ) ∝ O(n 12 so ).This is because we need to evaluate the quantum circuit(s) comprising O(N P ) parameters to evaluate each matrix element [32].
[Direct VQS] One does not have to evaluate the elements of the large matrix and vector, in contrast to the original VQS.The computational complexity scales as O(N depth N P N h N iter ) ∝ O(n 12 so N iter ), where N iter is the number of iterations required for optimization.If N iter does not depend strongly on N P , the direct VQS is as scalable as original VQS in terms of computational complexity.One has to evaluate the first and second terms of the last line of Eq. ( 32) and subtract the second term from the first one.This may be unstable under the influence of noise.
For comparing the computational complexity of VQE and that of direct VQS, both are roughly equivalent in terms of computation complexity because of the VQE's computational complexity O(N depth N P N h N iter ) ∝ O(n 12 so N iter ) with the UCCGSD ansatz.A more efficient and compact ansatz with a better scaling is desired.

Quantum impurity models
For a quantum impurity model with a starlike geometry (see Fig. 1), N h scales as O(N bath + N bath N imp + N 4 imp ), where N imp and N bath are the number of impurity orbitals and the number of bath orbitals, respectively.
For discretizing a continuous bath, N bath ∝ N imp are known to suffice [19,21] In Stage 4, we compute the imaginary-time Green's function where E G is the (approximate) ground-state energy obtained by VQE.We measure the quantity Ψ G | Â± |Ψ IM by using the same circuit as Eq. ( 18).

IV. NUMERICAL RESULTS
As an application of the two proposed algorithms using VQS and direct VQS, we solve a dimer model and a foursite impurity model obtained by DMFT calculations for the Hubbard model using a quantum circuit simulator.

A. Numerical details
In this study, we used Qulacs [37], pyed [38], Openfermion [39], and irbasis [40] to implement the proposed method.Qulacs is used as a quantum circuit simulator.We use the pyed library, which is based on TRIQS (Toolbox for Research on Interacting Quantum Systems) [41], for computing the reference data of the Green's function.Openfermion is used in the Jordan-Wigner transformation and to calculate the exact eigenvalues of models.We generate a sparse sampling mesh for a sufficiently large β = 1000 using irbasis (see Appendix C for more details).The sparse mesh covers −β/2 < τ < 0 and 0 < τ < β/2, when computing G(τ ) for τ < 0 and τ > 0, respectively.We perform DMFT calculations using DCore [42] to generate the four-site impurity model.
To deal with the numerical instability of the computing imaginary-time Green's function, we solve Eq. ( 27) using a truncated singular value decomposition as proposed in Ref. 23 and adopted additional tricks.See Appendix D for more details.
In the following results, we used the unitary coupled cluster ansatz with generalized singles and doubles (UC-CGSD) as the quantum circuit.We used the quasi-Newton method (the BFGS method) for optimizing the variational parameters.The gradients of cost functions are computed by a finite difference method.In Eq. ( 30), we set ∆τ = 10 −5 .In this study, we used random initial guesses as we observed that optimization was sometimes trapped in a metastable solution when starting from a zero initial guess.
In the following calculations, we used an MPIparallelized program.We used a workstation equipped with an AMD EPYC 7702P 64-core processor.To solve the largest model where the total number of circuit parameters is 1568, VQS took 10 hours with 22 cores, whereas direct VQS took 3 hours with 10 cores.

B. Dimer model
We first demonstrate the algorithms with a dimer model.Its Hamiltonian reads where niσ (≡ ĉ † iσ ĉiσ ) is the spin density operator at site i and spin σ.As illustrated in Fig. 4, the dimer comprises one interacting "impurity site" with an onsite repulsion U and one one-interacting "bath site."This corresponds to the case where there is one impurity site and one bath site (Fig. 1).We take U = 1, µ = U/2, b = 1 and V = 1.This model is not particle-hole symmetric because b = 0.
The exact ground-state energies for the number of particles n = 1, 2, 3, 4 are −1, −1.5, 0.2, 2, respectively.Thus, the global ground state has n = 2.In particular, there is a finite large energy gap between the ground and excited states.energy around the number of particles 2, which results in an exponential decay of the imaginary-time Green's function.This indicates that the system is insulating.
Figure 6 shows the off-diagonal component of the Green's function computed for τ > 0. This was done using Â+ = c 1,↑ and B+ = c † 2,↑ .The results clearly demonstrate that the off-diagonal component can be accurately measured using the proposed methods.
Figure 7 shows the Matsubara Green's function transformed from the τ domain.The Matsubara Green's function computed by VQS and direct VQS agree with the exact result from low to high frequencies.At high frequencies, the Green's function decays as 1/iω n , which is consistent with the fact that the Green's function has a discontinuity of 1 at τ = 0 owing to the noncommutativity of the creation and annihilation operators.Next, we consider the particle-hole symmetric four-site "impurity" model defined by the Hamiltonian where niσ = ĉ † iσ ĉiσ (σ =↑, ↓), and k is an index for "bath sites".This model corresponds to the case where there is one impurity site and three bath sites (Fig. 1).As shown in Fig. 4, the correlated impurity site is coupled to all three bath sites through the coupling terms V k = [0.0,−1.26264, 0.07702, −1.26264].There is no direct coupling between different bath sites.We take µ = U/2, k = [0.0,1.11919, 0.0, −1.11919].These parameters were determined by single-site DMFT calculations of the single-orbital Hubbard model on a square lattice with an onsite repulsion of U = 4 at half filling and zero temperature.The critical value of the Mott transition is U 12. The exact ground-state energies of the model are −3.16,−5.31, −5.49, and −5.51 for the number of particles n = 1, 2, 3, 4. Thus, there is only a small energy gap (0.02 between n = 4 and n = 3, 5).
Figure 8(a) shows the diagonal component of the Green's function, G 1↑,1↑ (τ ), computed by VQS or direct VQS as well as the exact one .The imaginarytime evolution was performed on 70 sampling frequencies in the same manner as for the dimer model.As G ii (τ ) = −G ii (−τ ) holds owing to the particle-hole symmetry, we only show the data for τ > 0. The exact Green's function decays slowly for 10 τ 10 2 owing to the small gap at first and then vanishes exponentially for τ 10 2 .
Figure 8(a) clearly demonstrates that the Green's functions can be accurately computed using our algorithms.At τ = 0, the computed Green's function agrees with the exact value within an accuracy of 10 −5 .This error comes from the fitting in stage 2. As τ increases, the absolute error increases owing to the discretization error in τ .The relative error however seems to stay constant at large τ , indicating the numerical stability of the present method.
Figure 8(b) shows the exponential component e η(τ )+τ EG of Eq. ( 33) and the remainder computed by VQS.For 10 τ 10 2 , the remaining part depends on τ , demonstrating that the slow decay of G 1↑,1↑ (τ ) is determined not only by the exponential part but also by the transition amplitude.For τ 10 2 , the remaining part is constant, where e η(τ )+τ EG dominates the exponential decay of G 1↑,1↑ (τ ) at large τ .
Figure 9(a) shows the Matsubara Green's function transformed from the τ domain.The Green's functions computed by VQS and direct VQS are in good agreement with the exact function.The Green's function decays as 1/iω n for the same reason as the dimer model.To confirm that this error originates from the discretization error in the imaginary-time evolution, we performed a similar simulation using a finer mesh comprising 139 sampling points, which were constructed by taking the midpoints of the original sampling points.Figure 9(b) shows the result, indicating that the effect is due to the discretization error as the error here is smaller than that shown in Fig. 9(a).Simultaneously, we observed that the error in the imaginary-time Green's function was also reduced (figure is not shown).
As indicated by the thin vertical lines in Fig. 8(b), to avoid numerical instability, we performed an imaginarytime evolution using an adaptively generated mesh in the τ domain.We decreased ∆τ when the imaginary-time evolution became unstable (for more details, refer to Appendix D).For the 70 and 139 sparse sampling points shown in Fig. 9, we used totally 121 and 159 mesh points for the imaginary-time evolution.Figure 10 shows how the imaginary-time evolution fails without the adaptive procedure.At τ 5, the imaginary-time evolution becomes unstable, which is signaled by the sudden increase in E τ .There are two possible reasons for this.First, the first-order approximation of the Taylor expansion in Eq. ( 30) is no longer a good approximation because ∆τ is extremely large.The second possibility is that the parametrization of the UCCGSD ansatz is redundant.This results in arbitrariness in the gradient θi (τ ), which makes Eq. ( 27) ill-posed.In other words, the condition number of M diverges.This may be because our test cases are still limited to small systems owing to the expensive computational cost of the imaginary-time on a quantum circuit simulator.To investigate this problem, solving a larger system with the same ansatz is necessary.Another interesting approach may be to construct a more compact ansatz with fewer parameters.Estimating the minimal number of shots for measurements assuming a fault-tolerant quantum computer is important.A previous study provides an analytic expression on the estimate of the total number of measurements for computing Green's functions by real-time evolution by VQS [23].
However, we cannot apply their expression to the computation of G(τ ) in a straightforward way due to the following reasons.The estimate of the number of shots depends on many factors such as the ansatz, the grouping of the terms of the Hamiltonian in measurements, and the condition number of the linear equation in Eq. (27) in VQS.In particular, the condition number of the matrix M in Eq. ( 28) also strongly depends on the system and the ansatz used.Actually, the condition number is divergent in the present cases.Thus, we somehow have to resort to numerical simulations.
Instead of estimating the total number of measurements for computing G(τ ) by imaginary-time evolution, we performed numerical simulations on the stability of our algorithm against shot noise assuming a fault-tolerant quantum computer.To be more specific, we studied how shot noise in the matrix and vector elements at an each time step [Eq.(28) and Eq. ( 29)] affects the imaginary-time evolution by VQS and the computed G(τ ).The shot noise was emulated by adding Gaussian noise with mean 0 to each element.The width of the distribution was set to the product of σ and the exact value, where σ (> 0) denotes the relative amplitude of the shot noise.
We performed simulations with various values of σ for the dimer model and the four-site model.In the following simulation, we removed the methods to deal with the numerical instability such as the "energy convergence condition" and "additional imaginary-time points" as discussed in Appendix D.
For reference, we demonstrate how the same shot noise as VQS affected the ground energy obtained by VQE.The shot noise was emulated by adding Gaussian noise with mean=0 to the expectation value of the Hamiltonian.When the width of the distribution σ = [10 −5 , 10 −3 , 10 −1 , 1, 10, 100], the relative errors of the expectation values are [6.545×10−6 , 0.001, 0.065, 0.460, 3.217, 169.851], respectively.To keep the relative error in energy low enough compared to the energy gap between the ground state and the first excited state, it will be necessary to reduce the shot noise to σ = 10 −3 .We confirmed that the imaginary-time evolution is stable up to as large as σ = 10 −3 .For σ 10 −1 , the imaginary-time evolution becomes more unstable than in the case of σ = 0.
We demonstrated the VQE adding the Gaussian noise in the measurement of the expectation value of the Hamiltonian.When the energy relative error for the width of the Gaussian distribution σ = [10 −5 , 10 −3 , 10 −1 , 1, 10, 100], the relative errors in energy are [5.144×10−6 , 0.0004, 0.032, 0.5091, 5.3307, 65.449], respectively.To keep the energy error low enough to compare to the energy gap between the ground state and the first excited state, it will be necessary to reduce the shot noise to σ = 10 −3 .

V. SUMMARY AND DISCUSSION
We proposed a quantum-classical hybrid algorithm to compute the imaginary-time Green's function on quantum devices with limited hardware resources by applying the VQS, which has been used to calculate the ground state.Using the quantum circuit simulator Qulacs, we verified this algorithm by computing Green's functions for typical impurity models such as the dimer model and four-site impurity model obtained by DMFT.The imaginary-time Green's function and Matsubara Green's function obtained using our algorithm agree well with the exact solution.Furthermore, we efficiently computed the imaginary-time Green's function by using a nonuniform mesh to reduce the number of imaginary-time τ points.For numerical instabilities occurring in regions where the mesh width ∆τ is large, we also computed the Green's function stably by applying a technique of adaptively generating mesh and imposing an energy convergence condition.
Quantum algorithms for computing the Green's function on quantum devices with limited hardware resources have been actively studied in recent years, and the complexity of quantum circuits needs to be discussed.First, we discuss the scalability of our algorithm and compare it with other similar methods.In calculating the excited states and VQS in our algorithm, the measurement must be repeated many times in the transition amplitude algorithm.This requires only a single ancillary qubit, which is approximately twice the depth of the quantum circuit used to calculate the ground state of the VQE.However, it would require many two-qubit unitary gates, which are challenging to implement in NISQ devices.
Recently, variational quantum algorithms to directly obtain the Green's function in the frequency domain was developed [17,18].These algorithms may be more general in that they directly compute the real-time Green's function as well as the imaginary-frequency Green's function.Finite-temperature static and dynamic correlation functions of spin systems have been calculated using the quantum imaginary time evolution (QITE) algorithm on a five-qubit IBM quantum device [22].These algorithms require measuring the expectation value of the square of the Hamiltonian or an effective Hamiltonian, which may be computationally demanding for a larger impurity model, whereas our method does not require measuring these observables.It is an interesting question which method is more efficient and stable in the presence of realistic noise.
Finally, we discuss future directions.It is interesting to perform simulations under realistic noise conditions with error mitigation techniques [43,44] for comparing the efficiency of the recently proposed various methods.For clarifying the causes of the numerical instability in the imaginary-time evolution, it is desired to apply the proposed methods to larger impurity models using a more compact/efficient ansatz.Possible directions are tensor decomposition methods [45][46][47][48][49][50] and an adaptive variational quantum imaginary-time evolution (AVQITE) approach [51,52].

FIG. 7 :
FIG. 7: Matsubara Green's function computed for the dimer model.Panels (a) and (b) show the results for the real part and imaginary part, respectively.

FIG. 8 :
FIG. 8: Computed G 1↑,1↑ (τ ) for the four-site impurity model.Panel (a) shows the result for τ > 0. The vertical thin lines denote the additional imaginary-time points added to avoid numerical instability (see the main text and Appendix D).Panel (b) shows different contributions to the Green's function (τ > 0).

FIG. 9 :
FIG. 9: Matsubara Green's function computed for the four-site impurity model using 70 (a) and 139 (b) sparse sampling points.The data in (a) correspond to the imaginary-time data shown in Fig. 8.

FIG. 10 :
FIG. 10: Numerical instability in computing the imaginary-time Green's function for a four-site impurity model by VQS (τ > 0).The panel shows Green's function computed without the adaptive construction of the mesh (see Appendix D).The inset shows E(τ ) and the exact ground-state energy for N + 1 particles (black horizontal line).
[21]te that N imp N bath holds for quantum impurity models describing real materials.For example, N bath and N imp required for a clustered DMFT calculation of the iron-based superconductor LaFeAsO were estimated to be N imp = 40, N bath = 332[21].For such a large impurity model, N h scales as