Quantum subspace expansion algorithm for Green's functions

I. INTRODUCTION

• Green’s functions are difficult to compute for strongly correlated materials, where nonperturbative methods such as DMFT are needed.Conventional approaches face challenging computational costs, motivating quantum algorithms for Green’s functions.Source: §I · p. 1
• Real-time Trotter methods require many steps for fine energy resolution, while quantum phase estimation requires deep circuits, limiting near-term applicability.Variational alternatives reduce depth but depend on accurate circuit-parameter optimization.Source: §I · p. 1
• QSE represents a target wave function as a linear combination of individually quantum-prepared basis states, with coefficients obtained classically from a linear algebra problem.QSE circuits are deeper than traditional VQAs but shorter than quantum phase estimation and can be parallelized over quantum computers.Source: §I · p. 2, A. Ground state · p. 2
• The QSEG method extends KVQA by computing the Green’s function as a continued fraction and representing Krylov states within QSE.The basis states are prepared using real-time evolution with multiple Trotter step sizes.Source: §I · p. 2

A. Ground state

• QSE ground-state coefficients and energy are obtained by measuring Hamiltonian and overlap matrix elements on a quantum device, then solving the resulting linear system classically.Increasing the basis improves agreement with the true ground state, while keeping the basis small is central to efficiency.Source: A. Ground state · p. 2

B. Green’s function

• The Green’s function is obtained from operator expectation values, with the density of states given by the imaginary part of the retarded Green’s function.The formulation considers greater and lesser Green’s functions and derives the retarded function from them.Source: B. Green’s function · p. 3
• Lanczos constructs an orthogonalized Krylov basis in which the Hamiltonian is tridiagonal, yielding the Green’s function as a continued fraction.The Krylov basis spans powers of the Hamiltonian applied to the initial state, and iteration may stop when b_n^2 is approximately zero.Source: B. Green’s function · p. 3
• The conventional Krylov procedure is limited by exponential Hilbert-space growth, while variational quantum implementations introduce circuit-parameter optimization challenges as system size increases.QSEG replaces that optimization with a quantum-subspace construction.Source: B. Green’s function · p. 3
• QSEG projects each Krylov state onto a chosen quantum-subspace basis, becoming exact when that basis is complete in the Krylov space.The Krylov-state basis must represent all Krylov states, unlike the ground-state basis, which represents only one wave function.Source: B. Green’s function · p. 3
• Quantum devices measure the Hamiltonian, overlap, and transition matrix elements, while classical computation performs the Krylov recurrence and inserts the resulting coefficients into the continued fraction.This separates quantum matrix-element estimation from classical linear-algebra and recurrence calculations.Source: B. Green’s function · pp. 3–4

A. Ground state basis

• Real-time evolution from an easily prepared reference state constructs a basis spanning the Krylov space and converges exponentially toward the ground state when the reference overlaps it.The reference may be chosen from a Hartree–Fock or approximate VQE state.Source: A. Ground state basis · p. 4
• A single-step basis faces a depth-resolution trade-off: small time steps make neighboring states similar, while large steps can miss high-frequency features.For each system there is an optimal time step and basis size.Source: A. Ground state basis · p. 4
• The multigrid approach is motivated by the need to resolve long-time features without the long circuits required by a single fixed time step.It is proposed as an alternative to variational time-evolution approximations.Source: A. Ground state basis · pp. 4–5
• The two-level multigrid basis combines large steps of length (n_k+1)∆t with smaller intermediate steps to reach long total times while retaining fine time resolution.Figure 1 illustrates the large-step and smaller-step sequence used to construct one basis state.Source: Figure 1 · p. 4, A. Ground state basis · p. 5
• For fixed ∆t, increasing n_k enlarges the reached time without increasing maximum circuit depth, although it increases Trotter error and basis size.QSE is relatively tolerant of Trotter errors because approximate evolution can still add useful basis degrees of freedom.Source: A. Ground state basis · p. 5

B. Basis for the Green’s function

• Real-time evolution applied to |χ0⟩=c†|GS⟩ generates basis states that cover the space spanned by powers of the Hamiltonian acting on |χ0⟩.This mirrors the ground-state basis construction but targets the Krylov states entering the Green’s function.Source: B. Basis for the Green’s function · p. 5
• The Krylov-state basis has size n_ψ=2(˜n_l+1)(˜n_k+1)−1 and maximum circuit depth d_V(1+˜n_l).Tilded parameters distinguish this basis from the ground-state basis parameters.Source: B. Basis for the Green’s function · p. 5
• QSEG can obtain |χ0⟩ through QSE, a VQA, quantum phase estimation, the Rodeo algorithm, or direct preparation when the creation operator is decomposed into Pauli strings.The algorithm itself is independent of which of these approaches supplies |χ0⟩.Source: B. Basis for the Green’s function · p. 5
• For nonzero multigrid parameters, the total maximum number of Trotter steps is n_l+˜n_l+2, whereas zero parameters recover ordinary single-step Trotter evolution.This bounds the basis-construction step count for the combined ground-state and Green’s-function bases.Source: B. Basis for the Green’s function · p. 5
• The required matrix elements can be computed independently, enabling runtime reduction through parallel execution across multiple qubit sets or quantum computers.The number of matrix elements grows with basis size, so minimizing the basis for a target accuracy remains important.Source: B. Basis for the Green’s function · pp. 5–6

Model

• The AIM contains interacting impurity sites and noninteracting bath sites, with H_0 describing onsite and hopping terms and H_int local Coulomb interactions.The site ordering places impurities first and bath sites afterward, for N=n_imp+n_b sites.Source: §IV · p. 6

Fermion-to-qubit mapping

• The Jordan–Wigner transform maps each fermionic spin-orbital to a bosonic orbital for quantum-circuit implementation.The spin-orbital mapping uses separate indices for spin-up and spin-down modes.Source: §IV · p. 6

Time evolution

• Quadratic evolution under H_0 is exact with Givens rotations, while commuting interaction terms exponentiate exactly before Suzuki–Trotter splitting combines H_0 and H_int.The symmetrized expansion reduces the Trotter error to O(∆t^3).Source: §IV · p. 6

Circuit depth

• For a single-impurity AIM, combining Givens rotations and interaction evolution gives a circuit depth d_V=4N+4.The interaction contribution requires two CNOT gates, in addition to the 4N CNOT layers from the quadratic evolution.Source: §IV · p. 6

Reference-state preparation

• The reference state is prepared by filling particles in the diagonal orbital basis and rotating to the local basis with C†.Its preparation depth is d_0=2N+max(n_↑,n_↓)−1 when the multi-fidelity estimation state is included.Source: §IV · p. 6

Circuit and model schematics

• The circuit construction uses Givens-rotation circuitry for the operator C and an AIM layout containing one impurity coupled to bath sites.The supplied schematics identify these circuit and model structures.Source: Figure 2 · p. 7, Figure 3 · p. 7

DMFT mapping

• The Bethe-lattice Hubbard model is mapped by DMFT to an Anderson impurity model whose bath hybridization satisfies the self-consistent relation ∆(z)=γ^2G(z).The AIM Coulomb interaction equals that of the Hubbard model, and units are chosen with γ=1.Source: §V · p. 7

Finite-size implementation

• The bath is discretized to seven sites, giving one impurity plus seven bath sites and requiring 16 qubits for the emulator calculation.The simulations use the Qulacs quantum emulator and focus on half filling.Source: §V · p. 7

Simulation procedure

• The DMFT self-consistency loop starts from the noninteracting solution and is iterated to convergence.The reported results cover both metallic and insulating regimes.Source: §V · p. 7

Ground-state convergence

• Figure 4 tracks QSE ground-state energy error ∆E versus Trotter steps n_l for different Coulomb interactions U during the first DMFT iteration.The comparison uses ∆t=0.1 and n_k=0.Source: Figure 4 · p. 7, §V · p. 7

A. Ground state

• The U=0 Hartree–Fock reference becomes less suitable as U increases because its overlap with the ground state decreases.For larger U, a different reference such as the U → ∞ state may become advantageous.Source: A. Ground state · pp. 7–8
• Increasing nl systematically reduces the QSE ground-state error, reaching approximately 10^-5 to 10^-3 at nl = 30.The required accuracy depends on the desired Green’s-function accuracy and energy resolution.Source: A. Ground state · p. 7
• For U=8, ΔE < 10^-3 is reached at nl = 7 with nk = 3, whereas nk = 0 requires nl > 25.The two-level scheme reduces Trotter steps for a fixed target accuracy and basis size.Source: A. Ground state · p. 8
• The two-level multigrid approach gives smaller ΔE for a given nl across AIM parameters, including at DMFT self-consistency.Its larger Trotter error does not significantly impair ground-state accuracy because the evolved states serve as basis states.Source: A. Ground state · p. 8

B. Green’s function

• QSEG reproduces the DOS accurately in both metallic U=2 and insulating U=8 regimes, with exact agreement approached at ˜nl = 80.At lower ˜nl, the insulating case converges faster, while quantitative metallic differences emerge around ˜nl = 9.Source: B. Green’s function · p. 9
• Keeping the basis size approximately constant, the two-level grid preserves good DOS agreement while reducing ˜nl; the required Trotter steps fall by about a factor of 9 for U=2.Accurate results use four Trotter steps for the insulating system and six for the metallic system, with a maximum circuit depth of 254 CNOT layers.Source: B. Green’s function · p. 10
• For the Green’s-function calculation, a single ground-state Trotter step is sufficient, allowing nl = 0 to minimize circuit depth.Increasing nk and ˜nk enlarges the basis without increasing circuit depth but requires more matrix elements.Source: B. Green’s function · p. 10

• The method uses quantum computers to represent Krylov states and obtain Hamiltonian and overlap matrix elements, overcoming the conventional Lanczos method’s exponential scaling barrier in principle.Source: §VI · pp. 10–11
• QSEG avoids the quantum-circuit parameter optimization required by variational quantum algorithms.Source: §VI · p. 10
• Two-level multigrid Suzuki–Trotter evolution prepares basis states with small time steps while tolerating Trotter inaccuracies because the states need only span the subspace accurately.Source: §VI · p. 11
• QSEG accurately computes Green’s functions for a 16-qubit Anderson impurity model obtained from DMFT for the Hubbard model on an infinite-connectivity Bethe lattice.Source: §VI · p. 11
• 4–6 Trotter steps suffice for the demonstrated insulating and metallic states, with circuit depths significantly below time-evolution-based Green’s-function methods.Source: §VI · p. 11

Matrix-element estimation

• Multi-fidelity estimation computes unitary-overlap matrix elements without Hadamard tests when a reference state satisfies the required overlap and phase conditions.The approach obtains the needed overlaps from state fidelities measured on a quantum computer.Source: Appendix A · p. 11

Matrix elements of S

• For particle-number-conserving Hamiltonians, choosing a reference state with zero or all particles makes the relevant Trotter-evolved reference states equal up to a phase.The phase is −h0T and vanishes when the Hamiltonian has no constant term.Source: Appendix A · p. 11

Matrix elements of H

• For H matrix elements, Pauli strings can change particle number, so the reference state is selected from the all-zero or all-one computational states according to the ground-state particle count.The choices exploit Pauli operators connecting only adjacent particle-number sectors and are intended to be easy to prepare.Source: Appendix A · p. 12

Krylov matrix elements

• The Krylov-state construction generalizes the same estimation technique to distinct states generated by Pauli strings, provided both have zero overlap with the reference state.The corresponding Hψ matrix elements require additional care because Hamiltonian Pauli strings may break particle-number symmetry.Source: Appendix A · pp. 11–12
• Because H can contain up to four particle-number-changing Pauli operators, the reference-state prescription is valid only when the total number of sites exceeds four.For ground-state particle count above four, the all-zero state is used; for counts at most three, the all-one state is used.Source: Appendix A · p. 12

GHZ-state preparation

• The state |R⟩+|φ0⟩ is prepared by applying a Hadamard gate, spin-sector CNOT chains, and the C† rotation from Eq. 34.The CNOT chains use n↑−1 and n↓−1 gates for the spin-up and spin-down sectors, respectively.Source: Appendix B · p. 12

Circuit depth

• The circuit requires max(n↑, n↓) −1 CNOT layers for GHZ preparation plus 2N layers for the C† preparation.Thus, the total CNOT-layer depth is the sum of these two contributions.Source: Appendix B · p. 12

Bethe-lattice hybridization

• The noninteracting Bethe-lattice hybridization is Δ0(ω)=1/4−ω2 under the convention γ=1.Here ω is the real frequency.Source: Appendix C · p. 12

Bath fitting

• The bath parameters are optimized by minimizing the squared Matsubara-axis distance between the finite-bath and target hybridizations.The fit uses ωn=(2n+1)π/β with β=100 and nM=100, emphasizing the low-energy region.Source: Appendix C · p. 12

DMFT iteration

• After determining the Green’s function, the procedure updates the hybridization through the DMFT equation and repeats until Green’s-function convergence.Tables I and II report bath parameters from the first and last DMFT iterations.Source: Appendix C · p. 12, Table I · p. 13, Table II · p. 13

Model setting

• At half-filling, the impurity-site energy is fixed by the interaction strength as ε11=−U.This is stated as an assumption for the half-filled setup.Source: Appendix C · p. 13

• The Sec. IIB algorithm cannot be directly applied to off-diagonal elements of the greater Green’s function.Source: Appendix D · p. 13
• The limitation concerns the off-diagonal components of the greater Green’s function.Source: Appendix D · p. 13
• Appendix D addresses off-diagonal greater-GF elements beyond the algorithm’s direct applicability.Source: Appendix D · p. 13