Comparative study on compact quantum circuits of hybrid quantum-classical algorithms for quantum impurity models

Predicting the properties of strongly correlated materials is a significant challenge in condensed matter theory. The widely used dynamical mean-field theory faces difficulty in solving quantum impurity models numerically. Hybrid quantum--classical algorithms such as variational quantum eigensolver emerge as a potential solution for quantum impurity models. A common challenge in these algorithms is the rapid growth of the number of variational parameters with the number of spin-orbitals in the impurity. In our approach to this problem, we develop compact ansatzes using a combination of two different strategies. First, we employ compact physics-inspired ansatz, $k$-unitary cluster Jastrow ansatz, developed in the field of quantum chemistry. Second, we eliminate largely redundant variational parameters of physics-inspired ansatzes associated with bath sites based on physical intuition. This is based on the fact that a quantum impurity model with a star-like geometry has no direct hopping between bath sites. We benchmark the accuracy of these ansatzes for both ground-state energy and dynamic quantities by solving typical quantum impurity models with/without shot noise. The results suggest that we can maintain the accuracy of ground-state energy while we drop the number of variational parameters associated with bath sites. Furthermore, we demonstrate that a moment expansion, when combined with the proposed ansatzes, can calculate the imaginary-time Green's functions under the influence of shot noise. This study demonstrates the potential for addressing complex impurity models in large-scale quantum simulations with fewer variational parameters without sacrificing accuracy.


I. INTRODUCTION
Accurately predicting the properties of strongly correlated materials poses a significant challenge in condensed matter theory, including long-standing challenges in the field, such as the mechanism of high-temperature superconductivity [1,2].Simulating these strongly correlated materials is difficult due to quantum superposition, which exponentially increases the accessible Hilbert space with the number of particles.Even with quantum computers with more than a hundred logical qubits, simulating solids with large numbers of degrees of freedom is still challenging.Quantum embedding theories, such as dynamical mean-field theory (DMFT) [3,4] or density matrix embedding theory (DMET) [5,6], aims to address this issue by limiting the correlated degrees of freedom in solid materials based on local approximation.
In DMFT, widely used in condensed matter physics, the original lattice system is divided into impurities with local interactions and a dynamical environment called a bath.This model is called a quantum impurity model.A self-consistent calculation is performed to update the parameters associated with the bath until the local Green's function defined on impurity matches that of the original lattice system with the dynamical mean-field.DMFT allows us to compute the single-particle excitation spectrum and successfully describes transitions from metallic * sakurairihito@gmail.com to Mott insulating behavior.The biggest numerical bottleneck in DMFT calculations is solving the correlated quantum impurity models, specifically computing local Green's functions for these interacting problems.While state-of-the-art classical algorithms have been adapted for use as impurity solvers, such as tensor networks [7][8][9] or quantum Monte Carlo methods [10], their applications are limited to models with only a few impurity and/or bath orbitals [7][8][9].This challenge stems from the exponential increase in quantum entanglement entropy and the notorious negative sign problem.
To exploit the growing potential for solving quantum impurity models on quantum devices, quantum algorithms based on quantum phase estimation [11,12] and adiabatic algorithms [13,14] have been proposed [15].Their practical implementation, however, may take decades because it requires large-scale error correction schemes.This led to a growing interest in variational quantum algorithms [16,17] for near-term quantum computers with limited hardware resources, often dubbed 'noisy intermediate-scale quantum' (NISQ) devices [18].A number of proof-of-principle demonstrations of solving quantum impurity models using NISQ devices have been conducted [19][20][21][22][23].
In near-term quantum algorithms, such as those in the NISQ era, it is crucial to utilize limited hardware resources effectively.Therefore, there is a need to discretize a continuous bath with fewer bath sites.This reduction can be achieved through the use of imaginary time formalism in DMFT [24][25][26][27][28].For example, a recent estimate for 20-orbital impurity models for iron-based supercon-ductors indicates that about 300 bath sites are sufficient for accurate discretization in the imaginary-time formalism [26].
Once this finite Hamiltonian representation of the quantum impurity model has been found, it is now in principle amenable to solution on a quantum device.For variational quantum algorithms, the first challenge is to define an appropriate ansatz which is flexible enough to span the solution to the problem, able to be efficiently evaluated via unitary quantum gates, and where the number of variational parameters, N P , does not grow prohibitively as the number of spin-orbitals N SO increases.Physics-inspired ansatzes based on unitary coupled cluster (UCC) methods [29][30][31] are widely used in previous studies for quantum impurity models [23,32,33].Among the family of UCC methods, for the unitary coupled cluster with generalized singles and doubles (UCCGSD) [34], N P grows as O(N 4 SO ).The computational times for computing imaginary-time Green's function grow even more rapidly, e.g., as O(N depth N 2 P ) [35] using the UCCGSD [34] and the variational quantum simulation (VQS) [17,36], where the depth of the circuit N depth ∝ N P .Thus, more compact ansatzes (circuits) are an important research direction for the success of simulating impurity models on quantum devices.
In this study, we develop compact ansatzes using a combination of two different strategies.First, we employ the k-unitary coupled Jastrow (k-uCJ) ansatz originally proposed for quantum chemistry, where N P scales only as O(N 2 SO ) [37].Second, we drop largely redundant variational parameters in both the UCCGSD and the k-uCJ ansatz based on physical intuition.This exploits structures in the Hamiltonian which are specific to quantum impurity models with a star-like bath geometry where the bath sites are connected via the Hamiltonian only through the impurity (see Fig. 1).In particular, we eliminate part of the two-particle excitations associated with direct excitations between bath sites, which does not change the scaling of N P but reduces the coefficient for a large number of bath sites.The scalings of the proposed ansatzes are summarized in Table I.We numerically demonstrate that the compact ansatzes describe ground-state energies and dynamic quantities, especially imaginary-time Green's functions, without compromising accuracy for typical quantum impurity models with/without shot noise, validating their potential in quantum impurity models.
The following outlines the contents of each section.Section II provides an overview of Green's functions and variational quantum algorithms for computing groundstate energy and dynamic quantities.This section also introduces the physics-inspired ansatzes used in this study.Section III introduces compact quantum circuits for quantum impurity models and compares the scaling of their variational parameters to those of the original ansatzes.Section IV compares the accuracies of groundstate energy and dynamic quantities such as spectral functions and imaginary-time Green's functions among ansatzes for typical quantum impurity models.Section V explores the effect of finite shot noise within the singlesite impurity model.Section VI reviews our results, compares them to existing methods, and highlights areas for future research.

II. REVIEW OF GREEN'S FUNCTIONS AND VARIATIONAL QUANTUM ALGORITHMS
A. Green's function We study a fermionic system in the grand-canonical ensemble, represented by the Hamiltonian H, with where c i /c † i are annihilation/creation operators for spinorbital i, and N represents the total number of spinorbitals.The hopping matrix, Coulomb interaction tensor, and chemical potential are denoted by t ij , U ijkl , and µ, respectively.The retarded (fermionic) Green's function is defined as where ĉa (t) = e iHt ĉa e −iHt and ĉ † b (t) = e iHt ĉ † b e −iHt represent the annihilation and creation operators for the spinorbitals a and b, respectively, in the Heisenberg representation.The θ(t) denotes the Heaviside step function.In this paper, we use ℏ = k B = 1.The thermal expectation, symbolized by ⟨• • •⟩, is evaluated in the grand canonical ensemble.
The retarded Green's function can be continued to the real frequency axis as where ω is a real frequency, while the imaginary-time Green's function is defined as where ĉa (τ ) = e τ H ĉa e −τ H .Note that the imaginarytime Green's function is anti-periodic as G ab (τ + β) = −G ab (τ ).The Fourier transform of the imaginary-time Green's function, known as the Matsubara Green's function, is given by where ω = (2n + 1)π/β and n ∈ N and β = 1/T .
The Matsubara 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 spectral representation with where z is a complex number and n, m runs over all eigenstates of the system with E m and E n being corresponding eigenvalues of H. On the real axis, these eigenvalues define individual poles for a finite system, or combine to form a branch cut for an infinite system.The retarded and advanced Green's functions are given by the value of G ab (z) just above/below the real axis.
Due to the branch cut on the real axis, G R ab (ω) ̸ = G A ab (ω) in general.The following relationship holds between the spectral function and the retarded and advanced Green's functions: where we used the formula 1/(x+i0 + ) = P(1/x)−iπδ(x), and P stands for the principal value.We now consider the limit of T → 0, where the ensemble average is restricted to the ground state(s) Ψ G .At sufficiently low temperatures, Eq. ( 4) can be rewritten as where A + = ĉa and B + = ĉ † b for 0 < τ < β/2, and A − = ĉ † b and B − = ĉa for β/2 < τ < 0. The signs ∓ are for τ > 0 and τ < 0, respectively, and In the presence of degenerate ground states, Eq. ( 11) should be averaged over all such states.In general, |G ab (τ )| decays exponentially in an insulating system, while algebraic in a metallic system.To ensure that G ab (τ ) is sufficiently small at the boundary, we need to increase β, which determines the upper limit of time evolution.

B. Variational quantum algorithms
In quantum computing, it is necessary to convert fermionic operators into qubit representations.There are several methods for this, such as the Jordan-Wigner transformation [38], and the Bravyi-Kitaev transformation [39,40].In this study, we use the Jordan-Wigner transformation given by ĉ 1. Ground-state calculation using VQE We use variational quantum eigensolver (VQE) [16,41].It begins by preparing an initial state |Ψ init ⟩ on a quantum computer.Then, a unitary operator described by a parameterized circuit with variational parameters θ, denoted as U (θ), is applied to the initial state, producing a quantum state, |Ψ(θ)⟩.Subsequently, the expectation value of each term in the Hamiltonian is measured using the quantum computer.This measured data is accumulated to compute the total expectation value of the Hamiltonian, ⟨H⟩, on a classical computer.The variational parameters are updated on the classical computer to minimize ⟨H⟩, and the process is iterated until the parameters are stably minimized.Provided the ansatz has sufficiently high expressive power and the optimization is carried out well using an appropriate initial state, the variational quantum state |Ψ(θ * )⟩ with optimized variational parameters θ * approximates the ground state |Ψ G ⟩ accurately.The success of the VQE therefore relies on finding an appropriate representation of the quantum state in terms of a sufficiently compact parameteric quantum circuit that can be optimized classically.

Recursive VQE for spectral moments
We detail here an approach to extend the scope of VQE to optimize the dynamics of the single-particle excitation spectrum via a compact moment expansion.This expansion allows access to a causal imaginary-time Green's function directly at zero temperature and in a fashion that allows for efficient quantum computation via a modified VQE [42][43][44].In a recent paper, direct measurements of the moment expansion expectation values via VQE have been proposed to compute the Green's functions [22].However, the proposed approach required measuring an increasing number of Pauli terms at higherorder moments and as systems increase in size, which we aim to mitigate via a recursive VQE approach to avoid this issue, as we will detail below.
The key physical quantities we aim to compute on the quantum device are the spectral moments of the Green's function.This quantity, which is classified as either hole or particle type at zero temperature, is defined in each case at the order m as follows: where These can be related to the matrix-valued spectral function, A(ω) rs defined in Eq. (10), as: The spectral moments defined in Eqs. ( 14) and (15) correspond to the Taylor expansions of the imaginarytime Green's function at the discontinuity points τ = 0 − and τ = 0 + , respectively.By increasing the number of moments, the imaginary-time Green's function can be systematically approximated over longer times τ .
Once the spectral moments for the particle and hole sectors are determined up to a maximum order N mom , we can appeal to the block Lanczos algorithm [45] to constructively build an effective single-particle Hamiltonian from these moments.This single-particle Hamiltonian spans the physical system and couples to it an auxiliary system whose dimensionality grows linearly with the number of system degrees of freedom and N mom .This auxiliary system acts as a zero-temperature dynamical self-energy, allowing correlation-driven changes to the original spectrum.These changes result from the projection of the eigenstates of this effective Hamiltonian back into the physical system.This auxiliary space is built in such a way that the resulting spectrum is causal, obeys required sum rules, and exactly preserves the initially provided moments, according to Eqs. ( 16) and (17).The resulting Green's function can be obtained directly in the Lehmann representation from the diagonalization of this effective Hamiltonian, providing the residues and energies of all the poles and allowing the Green's function to be easily transformed into any domain, including imaginary time.For more details of this procedure, see Refs.42-44, while similar approaches has also recently been applied in classical perturbative electronic structure methods to expand the self-energy [46,47].
We describe the procedure for calculating the moments defined by Eqs. ( 14) and (15) using a hybrid quantumclassical optimization algorithm, similar to VQE approach for the ground state.We assume that approximated |Ψ G ⟩ and E G are already computed using VQE.To simplify the exposition, we describe the construction of the particle sector moments, with the hole moments computed analogously.
First, we prepare a variational quantum state for the single-particle excited state ĉ † s |Ψ G ⟩.Because the operator is not unitary, we represent the resultant state as the action of a unitary multiplied by a scalar as where d 0 is a coefficient and the parametrized quantum state ϕ EX (θ 0 EX ) is defined by We choose to construct this state by defining an initial state ϕ 0 EX with N + 1 electrons and ensure that our parameterization for U (θ 0 EX ) conserves the electron number of the state.
The variational parameters θ 0 EX and coefficient d 0 can be computed as follows: After transforming ĉ † s into the qubit representation, we measure the cost function on the quantum computer via a circuit similar to a Hadamard test [48,49] (see Appendix A).The variational parameters are optimized to minimize the cost function C until convergence is achieved.After this optimization, the scaling coefficient measured on the quantum device.Finally, the zeroth order moment can be computed via the sampling of M p,(0) rs We can then subsequently compute the higher order moments up to N mom with (1 ≤ m ≤ N mom ) via a recursive approach, avoiding the need to measure over increasingly large numbers of Pauli strings for higher-order moments, as considered in Ref. 22. Using ϕ The variational parameters θ m EX and constant coefficient d m are determined by minimizing the cost function, . By performing m VQE steps optimizing these states, we can calculate the moments of order m as Similar ideas of hybrid quantum-classical variational optimization of alternative functionals for computing other (e.g.dynamical) properties have also been considered in other works [17,48,[50][51][52][53] As the ansatz used in optimizing all m states |ϕ m EX (θ m EX )⟩ becomes complete, it should enable the computation of the exact moments up to order m using the described approach.However, this optimization is also subject to various types of noises, including finite sampling errors of expectation values in a physical device, as well as optimization bottlenecks.This can result in numerical errors, which would likely accumulate exponentially at high orders of m.Nevertheless, as the magnitude of the moment also increases exponentially with respect to its order, we find that the numerical relative error in these moments compared to their exact benchmarks remains almost constant (see Appendix B).Finally, we note that while this approach has been presented for the computation of single-site Green's functions and moments, off-diagonal elements corresponding to matrix-valued Green's functions are possible, analogously to the approaches in Refs.22 and 48.

UCCGSD
The UCCGSD is a generalization of a unitary coupled cluster (UCC) [56][57][58][59][60][61] written as the exponential of an antisymmetric sum of excitation operators.The UCCGSD is formulated as follows: where |Ψ init ⟩ represents a product state, while Tn (n = 1, 2) and their respective conjugates T † n are excitation operators.The excitation operators Tn are T2 = 1 4 pqrs,αβγζ where T1 is a single-particle excitation operator, and T2 is a two-particle excitation operator.The indices p, q, r, s represent spatial orbitals, and α, β, γ, ζ represent spin.The composite indices pα, qβ, rγ, sζ span all spin-orbitals N SO .In this study, we removed one-particle and twoparticle excitations that change total S z .The t αβ pq and t αβγζ pqrs are complex-number variational parameters.The number of variational parameters 4 ), where N imp represents the number of spin-orbitals of the impurity and N bath the number in the bath.
Computing ⟨Ψ UCCGSD | H |Ψ UCCGSD ⟩ is exponentially expensive on classical computers because it results in a non-truncating Baker-Campbell-Hausdorff expansion.In contrast, quantum computers can compute this expectation value directly.We use a Trotter decomposition to implement Eq. ( 22) on a quantum computer.Classical optimization of variational quantum algorithms can partially mitigate the Trotterization error [62,63], but does result in a dependence of the final state on the ordering of the individual excitation operators.As commonly done, we set the Trotter step to 1, resulting in where qβ ĉpα } |Ψ init ⟩, demonstrating that the UCCGSD ansatz incorporates single-particle basis rotations into its definition [64].

k -uCJ
Let us first define the unitary cluster Jastrow (uCJ) ansatz and then the k-uCJ ansatz [55].The uCJ ansatz is defined as follows: The matrix K is complex and anti-Hermitian.The matrix J is symmetric, and its elements are purely imaginary.The |Ψ orb ⟩ is the single-particle basis rotated state defined in Eq. ( 25).This ansatz preserves the particle's number and total S z .The scaling with ).The uCJ ansatz is motivated via a tensor decomposition process that compresses the generalized two-particle excitation operators in the coupled cluster method.This compression results in a set of operators with only two indices.Similar approaches based on tensor decomposition have been proposed in Refs.65-69.Equation ( 28) can be implemented without Trotterization, as it involves only commuting number operators.By performing the Jordan-Wigner transformation on the equation, this term ĉ † pα ĉpα ĉ † qβ ĉqβ can be simplified to 1  4 (1 − Z pα )(1 − Z qβ ).The k-uCJ ansatz differs from the uCJ ansatz in that the operators J and K are applied multiple times, resulting in the k-uCJ ansatz, where variational parameters for different i are independently optimized.In a quantum embedding calculation, a continuous hybridization can be discretized with a finite number of bath sites.In particular, for a star-like geometry, the bath sites are connected only through the impurity.The number of bath sites, N bath , required for an accurate discretization scales linearly with N imp , albeit with a significant prefactor (on the order of ten [26]).Given the significant number of variational parameters associated with the bath sites, reducing the number of these parameters is critical for efficient quantum simulation of impurity models.
We propose compact ansatzes for quantum impurity models with a star-like bath geometries.We assume that two-particle excitation operators associated with two-body coupling between bath sites are not critical in the description of the ground states and spectral moments, given that two-particle interaction terms in the Hamiltonian are localized to the impurity space, and no Hamiltonian terms directly couple the bath sites.The ansatzes incorporating this assumption are referred to as "sparse ansatzes".In the present study, we construct sparse ansatzes based on the UCCGSD and the k-uCJ.We call them sparse UCCGSD and sparse k-uCJ, denoted UCCGSD(S) and k-uCJ(S), respectively.
For the UCCGSD, we remove two-particle excitation operators that involve more than three bath orbitals.Examples of such operators that involve three or four bath orbitals are ĉ † for the sparse variant (refer to Table I).Although N bath is proportional to N imp [26], ensuring that the scaling with respect to impurity size remains the same, the significant computational savings still result since N bath ≫ N imp .
For the k-uCJ, we apply a similar motivation to remove the operators acting between different bath sites while keeping the two-particle excitation operators between the impurity and the bath.For example, ĉ † 1 ĉ1 ĉ † 2 ĉ2 is dropped, as illustrated in Figs.1(c) and (d).As summarized in Table I, the scaling of N P in the k-uCJ ansatz scales as O (N imp + N bath ) 2 , while N P in the corresponding k-uCJ(S) sparse ansatz scales as O N 2 imp .Again, the prefactor is substantially reduced when TABLE I: Number of variational parameters for the UCCGSD, UCCGSD(S), k-uCJ, and k-uCJ(S).

IV. STATE VECTOR SIMULATION
In this section, we benchmark the k-uCJ and the proposed sparse ansaztes for typical quantum impurity models.We consider both single-site and two-site impurity models with N bath = 3 and N bath = 6, respectively.All calculations in this section are based on state vector simulations of quantum circuits.

A. Numerical details
The calculations were performed using the following libraries: QCMaterialNew [70] was used as a quantum circuit simulator, which is a Julia wrapper of Qulacs [71].We used Openfermion [72] for the Jordan-Wigner transformation and to calculate the exact eigenvalues of Hamiltonians.We performed DMFT calculations using DCore [73] to generate the single-site impurity models.We used dyson [74] library, in order to compute the Green's functions poles and residues from the spectral moments, as well as benchmark exact spectral moments via exact diagonalization.
For optimizing the variational parameters, we used the BFGS algorithm.We initialized variational parameters with random numbers.We observed that setting the initial guess to zero could lead the optimization to converge to a metastable solution.For ground-state calculations with VQE using the k-uCJ, we increased the number of terms k in the ansatz one by one, reusing the optimized variational parameters.In practice, at the beginning of the VQE calculations with k terms, we randomized the variational parameters in K1 and Ĵ1 but set those in Ki and Ĵi (2 ≤ i ≤ k) to the optimized variational parameters obtained in the previous calculation with k−1 terms.This procedure ensures that the optimized energy decreases or remains nearly stable with an increasing number of terms in the k-uCJ.
It is worth noting that the initial parameters significantly influence the accuracy of the optimized ground state and spectral moments.For ground-state calculations, we conducted VQE multiple times, each with a different set of initial parameters, to find the best variational state for the ground state.We used this best variational state for computing spectral moments.
Simulations were executed using an MPI parallelized program on a workstation with an AMD EPYC 7702P 64-core processor.Solving the largest model with 16 qubits and about 750 variational parameters in the k-uCJ took about five days on 55 cores using VQE and the recursive approach.

B. Single-site impurity model
We consider the single-site impurity model with particle-hole symmetry and N bath = 3 illustrated in Fig. 2

(a). The Hamiltonian is given by
where d † 1σ (ĉ † kσ ) are the impurity (bath) degrees of freedom of the fermionic creation operator with σ =↑, ↓, and k is an index for bath sites.The U represents the onsite Coulomb repulsion, V k is the hybridization term, µ (= U/2) is the chemical potential, and ϵ k denotes the bath energy.
We obtained the bath parameters using self-consistent DMFT calculations on a square lattice at zero temperature for U = 4 (metallic phase) and U =

(insulating phase).
The nearest neighbor hopping parameter was set to 1.
For U = 4, we obtained V k = {−1.26264,0.07702, −1.26264} and ϵ k = {1.11919,0.0, −1.11919}.For U = 9, we obtained V k = {1.31098,0.07658, −1.38519} and ϵ k = {−3.26141,0.0, 3.26141}.For the k-uCJ and the k-uCJ(S), we varied k from 1 to 5 to check convergence.The markers represent the best results obtained by varying the initial parameters 50 times for each ansatz.The lightly shaded areas indicate the variation in converged results depending on the choice of initial parameters for each ansatz.In all four ansatzes, the best ground-state energies are well reproduced.We also confirmed that the k-uCJ reproduces the ground-state energy with smaller N P than the UCCGSD.Also, the results for the sparse ansatzes in Figs.3(a) and (b) show that reducing the variational parameters associated with bath sites does not compromise the accuracy of the ground-state energies.It should be noted that the sparse ansatzes are efficient even for the metal-like system (U = 4), where the electronic structure is very much delocalized across the bath sites.
The following summarizes the reduction in N P for each ansatz by using the sparseness.In the UCCGSD, N P is reduced from 334 to 104.In the k-uCJ, N P is reduced from {64, 96, 128, 160, 192} to {58, 84, 110, 136, 162} for k = 1, 2, • • • , 5. The k-uCJ(S) has a small reduction in the number of parameters for this system, but this reduction becomes more significant with increasing system size and complexity (see Sec. IV C 1).
Figures 4(a) and (b) show the reconstructed spectral functions using the moment expansion for U = 4 and U = 9, respectively.For the k-uCJ, we set k = 5.We computed the exact values of the moments using exact diagonalization (ED).As shown in Fig. 4(a), for U = 4, all the ansatzes can reproduce the peaks around ω = 0.However, the quality of reproduction drops for ω ≥ 2. These discrepancies primarily arise from numerical errors during the moment computations via recursive VQE due to the limited representational ability of ansatzes and the optimization issues.This fitting error in the recursive approach grows exponentially with N mom , which prevents systematic improvement of reconstructed spectral functions with increasing N mom .Indeed, we observed no improvement for N mom > 7, although knowledge of the exact moments up to order N mom = 5 is largely sufficient to converge the spectral function over all frequencies.
As shown in Fig. 4(b), for U = 9, by increasing N mom up to 7, all the ansatzes accurately reproduced the positions of peaks for ω ≲ 6.In general, an insulating system has fewer spectral peaks than metallic cases, allowing the moment expansion by the recursive approach to reproduce the peak positions more accurately.Still, there is some variation among the ansatzes, likely due to the fitting error, especially around the small peak near ω = 4.The spectral function shows a tiny peak near ω = 0 as shown in the inset of Fig. 4(b).This is due to the k = 3 bath site nearly decoupled from the impurity, being physically irrelevant.
Here, we aim to quantify the difference between the spectral functions reconstructed from the exact moments and those calculated via the recursive approach.To this end, we utilize the Wasserstein metric, quantifying a difference between two distributions [75,76].Figures 6(a) and (b) show the imaginary-time Green's functions computed from the reconstructed spectral function by the moment expansion for U = 4 and U = 9, respectively.We use the reconstructed spectral function from the exact moments for each N mom as reference.In the k-uCJ, we set k = 5.In computing the reference data, we filtered out peaks below ω ≤ 10 −2 that are physically irrelevant.
In Figs. 6, for both U = 4 and U = 9, the differences among the ansatzes become less pronounced in the imaginary-time Green's functions compared to the differences in the spectral function.In Fig. 6(a), for U = 4, imaginary-time Green's functions exhibit a power-law decay.This necessitates a higher N mom in the moment expansion.However, for τ > 5, we observed that increasing N mom did not improve the accuracy due to the exponential growth in the fitting error with N mom in the recursive approach.Only the UCCGSD(S) result seems to diverge from the rest.Nonetheless, its deviation starting at τ = 5 aligns with the trends observed in other ansatzes, displaying a comparable pattern.In Fig. 6(b), for U = 9, imaginary-time Green's functions exhibit an exponential decay.The Green's functions computed by the recursive approach, even at N mom = 5, match the reference data, suggesting a smaller N mom achieves convergence compared to the metallic case.

C. Two-site impurity model
We consider the two-site impurity model with particlehole symmetry and N bath = 6, shown in Fig. 2

(b). The Hamiltonian is given by
where t represents the hopping between the two impurities.For V = 0.5 and V = 0.1, we use common bath parameters: U = 4, µ = U/2, t = 1, and ϵ k = {1, 0, −1, 1, 0, −1}.The case of V = 0.5 is expected to be more metallic than V = 0.1.Figures 7(a) and (b) show |δE G | for V = 0.5 and V = 0.1, respectively.For the k-uCJ and the k-uCJ(S), k was varied from 1 to 5 to check convergence.We omitted the VQE calculation with the UCCGSD because of its prohibitively large number of variational parameters.As before, the markers represent the optimal results obtained from 20 variations of the initial variational parameters for each ansatz.The lightly shaded areas highlight the dependency of each ansatz on initial guesses.
In the three ansatzes, the ground-state energies are reproduced with comparable accuracy.
Considering N P , both k-uCJ and k-uCJ(S) are more efficient than UCCGSD(S).The results for the sparse ansatzes in Figs.7(a) and (b) show that we can reduce the number of the variational parameters associated with bath sites without sacrificing ground-state accuracy in the cluster impurity model.The sparse ansatzes are also applicable for the case of V = 0.5, which exhibits more metallic characteristics.For the k-uCJ, N P is reduced from [256,384,512,640,768] to [226,324,422,520,618]

Spectral functions
Figures 8(a) and (b) show the reconstructed spectral functions using the moment expansion for V = 0.5 and V = 0.1, respectively.We computed the reference data from the exact moments for each N mom using exact diagonalization.In the k-uCJ, we set k = 5.Similar to the previous subsubsection, the UCCGSD and UCCGSD(S) calculations were omitted due to the prohibitive number of variational parameters.
In Fig. 8(a), for V = 0.5, increasing N mom tends to enhance the representation of several spectral peaks.Yet, it remains challenging to comprehensively capture the entire structure, mainly due to the fitting error, observing no improvement beyond N mom = 7.In Fig. 8(b), for V = 0.1, by increasing N mom up to 5, all the ansatzes accurately reproduced the positions of several peaks for ω ≲ 4.This indicates that an insulating system with fewer spectral peaks offers the advantage of accurately determining peak positions.However, variations around the small peak near ω = 1 among the ansatzes likely result from the fitting error.The spectral function shows a small peak around ω = 0 as shown in Fig. 8(b).This originates from bath sites weakly coupled with the impurity, being physically insignificant.
Figures 9(a) and (b) show the computed Wasserstein FIG. 8: Computed A 1↑,1↑ (ω) for each N mom .Panels (a) and (b) show the results for V = 0.5 and V = 0.1, respectively.In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments using exact diagonalization.The spectrum for V = 0.1 has a tiny peak around ω = 0, as shown in the inset.
metrics between the spectral functions reconstructed from the exact moments at N mom = 7 and those computed using the ansatzes at each N mom for V = 0.5 and V = 0.1, respectively.Due to the influence of noise, the distances, especially for N mom ≧ 5, stay at higher values than those without shot noise.Still, the Wasserstein metric tends to decrease as N mom increases, which is consistent with the improved reproducibility of the spectral functions reconstructed by the moment expansion.
Figures 10(a) and (b) show the imaginary-time Green's functions computed from the moment expansion for V = 0.5 and V = 0.1, respectively, with the k-uCJ ansatz with k = 5.We computed the reference data from the reconstructed spectral function using exact moments for each N mom .In this computation, we removed the physically irrelevant peaks below ω = 10 −2 in the spectrum (see the inset of 8).In Fig. 10, for both V = 0.5 and V = 0.1, the differences among the ansatzes become less pronounced in imaginary-time Green's functions compared to the cases of spectral functions.In Fig. 10 power-law decay.We observed no improvement by increasing N P , likely due to the fitting error in computing spectral moments.In Fig. 10(b), for V = 0.1, imaginarytime Green's functions exhibit an exponential decay.The results with N mom = 5 agree with the reference data.

V. FINITE SHOT SIMULATIONS
This section investigates the effects of shot noise for the single-site impurity model with N bath = 3.We first optimize variational parameters for the ground state and the intermediate states in the computation of the spectral moments [Eqs.(18), (20)] using state vector simulations as detailed in Sec.IV.Then, we measure the expectation values of the Hamiltonian and the transition amplitude (21) for each order of the moment m with a finite number of measurements.It should be noted that the effect of the shot noise was not considered during the optimization steps.This noise affects the measured scalar values, the ground-state energy E G and coefficients d 0 , d 1 , • • • , d mom in the recursive approach [Eq.(21)].We set the number of measurements to 30,000.In all four ansatzes, statistical errors with a finite number of measurements reduce the overall accuracy compared to the results without shot noise (see Fig. 3).Still, the ground-state energies can be reproduced with comparative accuracy among ansatzes.The results for the sparse ansatzes in Figs.11(a) and (b) show that reducing the variational parameters associated with bath sites does not compromise the accuracy of E G .The accuracy of the k-uCJ (S) is lower than that of the k-uCJ for the metallic system (U = 4), which may be attributed to statistical error.show the reconstructed spectral functions using the spectral moment computed with shot noise for U = 4 and U = 9, respectively.We set k = 5 in the k-uCJ.In Fig. 12(a), for U = 4, none of the ansatzes reconstruct the spectral peaks.These discrepancies primarily stem from numerical errors in the moment calculations.It should be noted that reconstructing a spectral function from the moments is not a well-conditioned problem (although more robust than traditional numerical analytic continuation from imaginary time due to the analytic procedure).Specifically, in the shot noise simulation, such errors are attributed to statistical error, the limited representational capability of the ansatzes, and optimization issues.The effect of statistical noise is dominant when comparing the calculation results to the case without shot noise Fig. 4. In Fig. 12(b), for U = 9, the shot noise induces small shifts in the positions of several peaks for ω ≲ 6 compared to the results computed without the shot noise.There are some variations among the ansatzes, likely due to the fitting error, but generally, the agreement is much improved compared to the more metallic U = 4 results.

C. Imaginary-time Green's functions
We now compute the imaginary-time Green's functions from the reconstructed spectral functions by the moment expansion with the shot noise.Figures 13(a) and (b) show the results for U = 4 and 9, respectively.In the k-uCJ, we set k = 5.
For both U = 4 and U = 9, despite the large deviations in the spectral functions due to the fitting error, these variations are suppressed in the reconstructed imaginary-time Green's functions.The results from all the ansatzes are consistent up to τ ≈ 1; then they start to deviate.This is because the imaginary-time Green's function is insensitive to changes in the associated spectral function.In Fig. 13(a), for U = 4 with N mom = 7, due to the shot noise, the black vertical line at τ = 1 marks the earlier start of deviation, while the gray vertical line at τ = 5 indicates the start without shot noise (see Fig. 6).In Fig. 13(b), for U = 4 with N mom = 5, the results with shot noise are in good agreement with the reference data.These results indicate the moment expansion can successfully calculate the imaginary-time Green's functions under the influence of shot noise.The imaginary-time Green's function, as calculated in this way, is sufficient for performing self-consistent calculations of DMFT.After convergence, some quantities com- puted from the imaginary-time Green's function (e.g., electron occupancy) are expected to be less sensitive to noise than real-frequency spectral functions.

VI. SUMMARY AND DISCUSSION
In this paper, we proposed compact quantum circuits for quantum impurity models with a star-like bath geometry by sparsifying the UCCGSD and k-uCJ ansatz.These forms have a significant parameter scaling of N 4 SO and N 2 SO respectively, which are reduced by removing insignificant variational parameters associated with two-body coupling between bath sites.This results in a reduced number of variational parameters scaling as O(N 4 imp ) and O(N 2 imp ) for the UCCGSD(S) and k-uCJ(S) ansatz respectively.We numerically demonstrated that the compact ansatzes can accurately reproduce the ground-state energies for typical quantum impurity models, with and without shot noise.In the moment calculations for dynamic quantities, to avoid measuring more Pauli-operator terms at higher orders, we proposed a recursive method similar to VQE.We also demonstrated that, when combined with the suggested ansatzes, the moment expansion effectively computes the imaginary-time Green's function, even in the presence of shot noise.
Before concluding this paper, we consider the proposed ansatzes and the spectral moments in the context of other approaches.A previous study utilized an adaptive variational quantum eigensolver (ADAPT-VQE) for impurity models [33,77,78].While ADAPT-VQE can provide near-exact solutions with a deep circuit, it demands more measurements for gradient computation than traditional VQE.Also, its success depends on the selected operator pool, which makes it hard to compare it to other approaches.In addition, it is instructive to compare the moment expansion to alternative approaches such as the VQS approach [35], with which the method bears many similarities.The moment expansion preserves the causal nature of the spectral functions; however, it encounters growing fitting errors in the recursive approach, most significantly in metallic systems.The VQS method might handle these systems more effectively via time evolution over a longer time span.Still, it could be costly since it requires computing all variational parameters at every time step.A more detailed comparison is left in future studies.
Finally, we discuss the potential future research directions.The initial parameter selection plays a crucial role in the accuracy of ground-state energies and moments.Specifically, the accuracy of the moment is closely tied to that of the ground state.The selection of optimal initial parameters to avoid local minima should be a critical area of future studies.Moreover, minimizing the number of measurements in VQE and the recursive approach is essential in the context of utilizing near-term quantum devices.One viable solution is the efficient grouping of observables for simultaneous measurement [79].It is also important to investigate how the noise in the measured Green's function and discretization errors of the bath propagate during self-consistent calculations in DMFT and affect quantities of interest, e.g., momentum-resolved spectrum.Developing methods for suppressing such errors is crucial.When optimizing in VQE under the presence of noise, it is crucial to employ optimization methods resilient to noise [80,81].At the same time, noise mitigation techniques [82,83] are essential.Furthermore, the potential applicability of sparse ansatzes to other impurity models with star-like geometry, such as multi-orbital systems, requires further investigation.Lastly, incorporating the concept of sparsity into classical variational algorithmic approaches, such as machine learning wave functions [84,85], may improve computational efficiency.In Fig. 15(a), for U = 4, no significant difference in relative error between ansatzes was observed due to the shot noise.Still, the relative error for each ansatz remains nearly constant.
FIG. 1: Schematic illustrations for the construction of sparse ansatzes.Panels (a) and (b) show the eliminated operators that involve more than three bath orbitals when constructing the UCCGSD(S) from the UCCGSD.Panels (c) and (d) show the eliminated operators acting between the different bath sites when constructing the k-uCJ(S) from the k-uCJ.

1 ĉ † 1 N 2 bath≃ O N 4 imp
ĉ2 ĉ2 and d † ĉ † 3 ĉ3 ĉ4 , where d † (ĉ † ) are fermionic creation operators for the impurity (bath) degrees of freedom, respectively.We illustrate these operators in Figs.1(a) and (b).For the UCCGSD, this reduces the number of variational parameters N P from O (N imp + N bath ) 4 to O N 4 imp +

FIG. 2 :
FIG. 2: Two quantum impurity models used in this study.Panels (a) and (b) show the single-site impurity model with N bath = 3, and the two-site impurity model with N bath = 6, respectively.

= 9 FIG. 3 :
FIG. 3: Computed |δE G | for the single-site impurity model.Panels (a) and (b) show the results for U = 4 and U = 9, respectively.In the k-uCJ and the k-uCJ(S), k increases from 1 to 5. The markers represent the smallest errors when the initial parameters are changed 50 times.The lightly shaded areas in the figure illustrate the dependency of the absolute errors on the initial parameters.

Figures 3 (
Figures 3(a) and (b) show the absolute errors in ground-state energies (|δE G |) for U = 4 and U = 9, respectively, compared to exact diagonalization results.For the k-uCJ and the k-uCJ(S), we varied k from 1 to 5 to check convergence.The markers represent the best results obtained by varying the initial parameters 50 times for each ansatz.The lightly shaded areas indicate the variation in converged results depending on the choice of initial parameters for each ansatz.In all four ansatzes, the best ground-state energies are well reproduced.We also confirmed that the k-uCJ reproduces the ground-state energy with smaller N P than the UCCGSD.Also, the results for the sparse ansatzes in Figs.3(a) and (b) show that reducing the variational parameters associated with bath sites does not compromise the accuracy of the ground-state energies.It should be noted that the Figures4(a) and (b) show the reconstructed spectral functions using the moment expansion for U = 4 and U = 9, respectively.For the k-uCJ, we set k = 5.We computed the exact values of the moments using exact diagonalization (ED).As shown in Fig.4(a), for U = 4, all the ansatzes can reproduce the peaks around ω = 0.However, the quality of reproduction drops for ω ≥ 2. These discrepancies primarily arise from numerical errors during the moment computations via recursive VQE due to the limited representational ability of ansatzes and the optimization issues.This fitting error in the recursive approach grows exponentially with N mom , which prevents systematic improvement of reconstructed spectral functions with increasing N mom .Indeed, we observed no improvement for N mom > 7, although knowledge of the exact moments up to order N mom = 5 is largely sufficient to converge the spectral function over all frequencies.As shown in Fig.4(b), for U = 9, by increasing N mom up to 7, all the ansatzes accurately reproduced the positions of peaks for ω ≲ 6.In general, an insulating system has fewer spectral peaks than metallic cases, allowing the moment expansion by the recursive approach to reproduce the peak positions more accurately.Still, there is some variation among the ansatzes, likely due to the fitting error, especially around the small peak near ω = 4.The spectral function shows a tiny peak near ω = 0 as shown in the inset of Fig.4(b).This is due to the k = 3 bath site nearly decoupled from the impurity, being physically irrelevant.Here, we aim to quantify the difference between the spectral functions reconstructed from the exact moments and those calculated via the recursive approach.To this end, we utilize the Wasserstein metric, quantifying a difference between two distributions[75,76]. Figures 5(a) and (b)show the computed Wasserstein metric between the spectral functions from the exact moments at N mom = 7 and those using ansatzes at each N mom for U = 4 and U = 9, respectively.As N mom increases, the distance between the two distributions decreases, consistently with the enhanced reproducibility of the spectrum at large N mom .

7 FIG. 4 := 9 FIG. 5 :
FIG.4: Computed A 1↑,1↑ (ω) for each N mom .Panels (a) and (b) show the results for U = 4 and U = 9, respectively.In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments using exact diagonalization.The spectrum of V = 0.1 around ω = 0 has a tiny magnitude of 10 −2 , as shown in the inset.

7 FIG. 6 :
FIG.6: Computed G (τ ) for each N mom .Panels (a) and (b) show the results for U = 4 and U = 9, respectively.In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments using exact diagonalization.The black vertical lines in panel (a) for N mom = 7 show where the deviation of the reconstructed spectral functions from the reference data starts.

1 FIG. 7 :
FIG. 7: Computed |δE G | for the two-site impurity model.The markers represent the smallest errors when the initial parameters are changed 20 times.In the k-uCJ and the k-uCJ(S), k increases from 1 to 5. The lightly shaded areas in the figure illustrate the dependency of the absolute errors on the initial parameters.Panels (a) and (b) show the results for U = 4 and U = 9, respectively.

1 FIG. 9 :
Figures10(a) and (b) show the imaginary-time Green's functions computed from the moment expansion for V = 0.5 and V = 0.1, respectively, with the k-uCJ ansatz with k = 5.We computed the reference data from the reconstructed spectral function using exact moments for each N mom .In this computation, we removed the physically irrelevant peaks below ω = 10 −2 in the spectrum (see the inset of 8).In Fig.10, for both V = 0.5 and V = 0.1, the differences among the ansatzes become less pronounced in imaginary-time Green's functions compared to the cases of spectral functions.In Fig.10(a), for V = 0.5, the imaginary-time Green's functions exhibit a

7 FIG. 10 :
FIG. 10: Computed G 1↑,1↑ (τ ) for each N mom .Panels (a) and (b) show the results for V = 0.5 and V = 0.1, respectively.In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments computed by exact diagonalization.The black vertical lines in panel (a) for N mom = 7 show where the deviation of the reconstructed spectral functions from the reference data starts.

= 9 FIG. 11 :
FIG. 11: Computed |δE G | with 30000 measurements for the single-site impurity model.In the k-uCJ and the k-uCJ(S), k was varied from 1 to 5. The markers represent the best result obtained by varying the initial parameters 50 times.The lightly shaded areas in the figure illustrate the dependency of the absolute errors on the initial parameters.Panels (a) and (b) show the results for U = 4 and U = 9, respectively.

7 FIG. 12 :
FIG. 12: Computed A 1↑,1↑ (ω) with measurements for each N mom .In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments using exact diagonalization.Panels (a) and (b) show the results for U = 4 and U = 9, respectively.

7 FIG. 13 :
FIG.13: Computed 1↑,1↑ (τ ) 30000 measurementsfor each N mom .Panels (a) and (b) show the results for U = 4 and U = 9, respectively.In the k-uCJ and the k-uCJ(S), we set k = 5.ED refers to the spectral functions constructed from exact moments using exact diagonalization.The black vertical lines in panel (a) for N mom = 7 indicate where the reconstructed spectral functions with shot noise begin to differ from those derived from exact moments.The gray line indicates the case without shot noise (see Fig.6).

Figures 16 (
Figures 16(a) and (b) show the relative errors of the spectral moments |δM p rs |/|M p rs | with a finite number of measurements, 30000 for U = 4 and U = 9, respectively.The markers in the figure denote the mean, and the lightly shaded areas indicate the standard deviation derived from the calculation repeated ten times with shot noise for each ansatz.The sparse ansatz is generally less accurate than the original ansatz due to the shot noise.In Fig.15(a), for U = 4, no significant difference in relative error between ansatzes was observed due to the shot noise.Still, the relative error for each ansatz remains nearly constant. qβ,rγ,sζ rγ ĉqβ ĉpα } |Ψ orb ⟩ ,