A Hybrid Quantum-Classical Method for Electron-Phonon Systems

Interactions between electrons and phonons play a crucial role in quantum materials. Yet, there is no universal method that would simultaneously accurately account for strong electron-phonon interactions and electronic correlations. By combining methods of the variational quantum eigensolver and the variational non-Gaussian solver, we develop a hybrid quantum-classical algorithm suitable for this type of correlated systems. This hybrid method tackles systems with arbitrarily strong electron-phonon coupling without increasing the number of required qubits and quantum gates, as compared to purely electronic models. We benchmark the new method by applying it to the paradigmatic Hubbard-Holstein model at half filling, and show that it correctly captures the competition between charge density wave and antiferromagnetic phases, quantitatively consistent with exact diagonalization.


I. INTRODUCTION
Understanding strongly correlated many-body systems is vital to many areas of science and technology, such as the development and analysis of functional quantum materials [1]. Due to the entanglement induced by correlations, macroscopic properties of quantum materials are often unpredictable from reductive single-particle models. Theoretical analysis of these systems with strongly entangled degrees of freedom has, however, been hindered by the exponential growth of their Hilbert space sizes with the number of particles. Understanding macroscopic properties of materials requires the analysis of sufficiently large model systems, which cannot be done accurately with classical computers. Quantum computing technologies, including hybrid quantum-classical algorithms [2][3][4] constitute an intriguing new direction for studying quantum many-body systems and especially quantum materials.
One of the promising representatives of hybrid algorithms is the variational quantum eigensolver (VQE) [2,[5][6][7][8]. This approach aims to accurately approximate ground states of quantum systems that can be naturally represented using qubits, such as spin and fermionic models. An example of such a protocol is shown in the upper panel of Fig. 1a: one uses a set of parameterized quantum gates to prepare a variational wavefunction and measures the expectation value of the Hamiltonian; then, one optimizes the parameters of these quantum gates using a classical computer. VQE has been implemented experimentally for small molecules [9][10][11][12][13], providing accu- * yao.wang@emory.edu rate solutions verifiable by exact methods. Hardwareefficient implementations of VQE have also been proposed for solid-state systems, including quantum magnets and Mott insulators [14][15][16][17]. Successful applications of the VQE approach to describe multi-orbital molecules and intermediate-size solid-state models pave the way for extending this technique to broader classes of materials.
However, realistic materials usually contain more complex interactions than simplified electronic models, such as the Hubbard model which only features local Coulomb interaction. The interaction between mobile electrons and the ionic lattice in solids, so-called electron-phonon coupling (EPC), underlies electric and mechanical properties of materials. Notably, it has been suggested that the interplay between the electron-phonon interaction and the electronic Coulomb repulsion is crucial for many novel quantum phases, such as unconventional superconductivity in cuprates [18][19][20][21][22][23] and twisted bilayer graphene [24][25][26][27][28][29]. Achieving predictive control of these quantum phases calls for developing reliable theoretical models describing materials with EPC [30][31][32], which has motivated studies based on small clusters [33][34][35] or perturbative couplings [36][37][38][39]. Quantum simulation of materials with strong EPCs, however, remains challenging due to the unbounded phonon Hilbert space [40][41][42][43]. The common spirit of quantum algorithms is to traverse quantum states encoded by the combination of available qubits. Therefore, even with a single electronic band and single phonon mode, the system has much higher computational complexity compared to electrons alone: inclusion of phonons in an L-site spinful system increases the Hilbert-space size from 4 L (for electrons) to 4 L (m + 1) L where m is the (truncated) maximal local phonon occupation. For materials with non-negligible EPCs, the required m ≫ 1 leads to an unreasonably large, and even unbounded Hilbert space. This issue prohibits not only classical simulations, but also an efficient encoding on a quantum machine.
To this end, we design a hybrid quantum algorithm which leverages the capability of VQE with quantum computers and the variational non-Gaussian description of non-perturbative polaronic dressing [44][45][46]. We prove the validity of our approach using the one-dimensional Hubbard-Holstein model and its variants, which is summarized together with the specifics of the algorithm in the "Methods" section. We then show that our hybrid quantum algorithm is able to reliably capture the groundstate properties of the paradigmatic Hubbard-Holstein model in all regions of the phase diagram, when compared to non-Gaussian exact diagonalization (NGSED) results. Our algorithm does not require any additional qubit overhead stemming from unbounded phononic degrees of freedom and the truncation to a low phonon occupation [40,42,43], as we implicitly sample the phonon Hilbert space. This makes it possible to analyze electronphonon systems over a broad range of parameters, including both adiabatic and anti-adiabatic regimes of the phonon frequencies. We conclude our analysis by investigating the scaling of the algorithm's performance with respect to the system size, indicating the reduction of exponentially increasing complexity. Moreover, unlike prior studies on EPC systems implemented in trapped ions [47,48], our approach is not limited to a particular hardware platform. This makes our algorithm a promising candidate for studying systems with strong EPC and electron-electron interactions beyond classically solvable problems.

A. Variational non-Gaussian VQE (NGS-VQE) method
We consider a prototypical correlated system, where electronic correlations are described by the local Coulomb interaction -the Hubbard model, while electron-phonon coupling follows the linear Fröhlich-type density-displacement interaction. While the latter is usually also simplified into local couplings -the Holstein model -recent experimental discoveries in cuprates have indicated the importance of nonlocal couplings [49][50][51][52]. Including all these interactions, we obtain the Hubbardextended-Holstein (HEH) model, whose Hamiltonian is Here, c i,σ (c † i,σ ) annihilates (creates) an electron at site i with spin σ, and a i (a † i ) annihilates (creates) a phonon at site i; n i,σ = c † i,σ c i,σ denotes the electron density operator for site i and spin σ. Among the model parameters, t sets the (nearest-neighbor ⟨i, j⟩) hopping integral, U sets the on-site repulsive interaction, ω 0 sets the Einstein phonon energy. g ij is the coupling strength between the phonon displacement at site i and electron density at site j. While our method can tackle any distribution of EPCs, as discussed later, we restrict ourselves to local g = g ii and nearest-neighbor coupling g ′ = g i,i±1 in one dimension.
To handle the strongly entangled electronic wavefunction and unbounded phonon Hilbert space simultaneously, we employ a variational non-Gaussian construction of the many-body wavefunction [44,69]. A universal electron-phonon wavefunction can always be written in the form of where the right-hand side is a direct product of electron and phonon states (denoted as |ψ e ⟩ and |ψ ph ⟩, respectively), with the variational non-Gaussian transformation U NGS = e iS and the (Hermitian) operator S being a polynomial formed by c, c † , a, and a † operators (with any sub-indices). The functional class of Eq. (2) is a complete representation of the full electron-phonon Hilbert space. The variational solution based on this approach is guaranteed to be accurate, provided that S can assume arbitrary polynomials. The accuracy usually converges at a relatively low order in the exponent [70]. When determining the order of polynomials in S, one should balance theoretical needs and computational feasibility: including high-order powers of electronic and phonon operators in S improve the accuracy and expand the applicability to complex models; however, these powers also lead to a large variational parameter space and a complex energy representation form. As benchmarked by exact diagonalization (ED) and determinant quantum Monte Carlo (DQMC) simulations of small clusters [45,46], it is sufficient to truncate the S operator to the lowest-order terms for the Holstein-type coupling [see Eq. (7) in the "Methods" section]. Denoting these lowest-order coeffi- cients as {f q } (q is the quantum number, naturally chosen as momentum for periodic systems), we have the variational non-Gaussian transformation U NGS ({f q }) fully determined by these parameters. Using this ansatz, we solve the HEH problem by minimizing the average energy self-consistently with respect to the unrestricted electronic state |ψ e ⟩ and the variational parameters in U NGS ({f q }) and |ψ ph ⟩. Within each iteration, the variational non-Gaussian parameters {f q } and the phonon wavefunction |ψ ph ⟩ (here restricted to be a Gaussian state) are optimized using the imaginary-time equations of motion derived in Ref. 45 (see Fig. 1a). This is referred to as the non-Gaussian (NGS) solver, whose computational complexity scales polynomially with the system size L. On the other hand, the fully entangled electronic part |ψ e ⟩ of the wavefunction, represented by a tailored quantum circuit, can be obtained by regarding the |ψ ph ⟩ and {f q } as fixed and further minimizing the total energy. The latter step is equivalent to solving the electronic ground state of an effective Hubbard Hamiltonian H eff = ⟨ψ ph |U † NGS H HEH U NGS |ψ ph ⟩. Physically, this H eff describes the behavior of polarons, formed by phonon-dressed electrons. The phonon dressing gives rise to a heavier effective mass by modifying the hopping strengtht and mediates a long-ranged attractionṼ ij < 0 between polarons, in the form of Eq. (9) in the "Meth-ods" section. In the case of Holstein couplings, Gaussian states are an efficient representation of the phonon wavefunction |ψ ph ⟩ [44,45,69]. Gaussian states allow to represent the effective nearest-neighbor hopping and phonon-mediated interactions in closed-form, as shown in Eqs. (10) and (11) of the "Methods" section (see also Ref. 50).
Within each self-consistent iteration, the key complexity of solving the electron-phonon-coupled problem has thus been converted into solving a purely electronic Hamiltonian H eff with extended Hubbard interactions V ij . This electronic problem can be efficiently embedded on a quantum hardware by using a suitable fermionic encoding. Here, we employ the Jordan-Wigner transformation, which maps each electron with given spin orientation to one qubit [see "Methods" section]. By applying a set of parameterized rotations to these qubits, we obtain a quantum circuit representing a variational electronic wavefunction |ψ e ({θ i })⟩. The self-consistent quantum-classical iterations in VQE then optimize the variational gate parameters {θ i } to minimize the energy ⟨ψ e ({θ i })|H eff |ψ e ({θ i })⟩. The solution of VQE approximates |ψ e ({θ i })⟩ in the variational ground state of Eq. (2). Unless explicitly specified otherwise, we conduct the VQE step of the NGS-VQE iterations with exact statevector simulations using Qiskit [71]. The standard Hubbard model has been successfully studied with various quantum circuits. Here, we employ the Hamiltonian variational ansatz [72], which can be naturally extended to nonlocal interactions: inspired by the quantum adiabatic theorem, it starts from an efficient encoding of the noninteracting lattice model with translational symmetry [73][74][75][76], and then evolves the state using alternating kinetic energy and interacting terms [15,17] [see "Methods" section]. The specific ansatz for the quantum circuit is presented in the upper panel of Fig. 1a. Since the ground state of a finite-size periodic system preserves translational symmetry, we assume the quantum gates in the same layer to share the same parameters (denoted as θ n in Fig. 1a). The expressibility of the quantum circuit ansatz is controlled by the repetition number n of the evolution block (P and H gates in Fig. 1a). Depending on the hardware specific transpilation, the number of CNOT gates for a given n follows as 6+5×n. This quantum circuit represents an efficient encoding if n does not scale exponentially with the system size L. We investigate the scaling with the system size in the "Scaling in circuit depth and system size" subsection, and show that the chosen ansatz leads to a quantitatively accurate performance of the hybrid quantum algorithm. Figure 1b shows an example of the NGS-VQE simulation for a 4-site Hubbard-Holstein model with u = 10 and λ = 10. Due to the variational nature of each step in the self-consistent NGS-VQE iteration, the total energy E always decreases. Thus, this NGS-VQE iteration is ensured to converge to a local energy minimum within the variational space set by Eq. (2), which provides a good approximation for the ground state. It is worth noting that the energy evolution is not always smooth, where a sudden drop indicates a dramatic change in the variational wavefunction from one phase to another. In the example presented in Fig. 1b, the initial state is prepared by setting the phonon wavefunction to vacuum. Consequently, the first outer-loop iteration with VQE (VQE#1 in the panel) starts with a pure Hubbard model, followed by the adjustment of phonon states and NGS parameters (NGS#1 in the panel). After the entire first outer-loop iteration (VQE#1 and NGS#1), the electronic state |ψ e ⟩ lies in an antiferromagnetic (AFM) state as a solution for the pure Hubbard model, while the phonon state |ψ ph ⟩ induces a large attractive potential in the form of Eq. (11) in the "Methods" section. This phonon-mediated interaction tends to stabilize a charge density wave (CDW), which contradicts the AFM state (see discussion in "Charge and spin phases in the Hubbard-Holstein model" for details). As a result, the electronic state rapidly evolves once the second self-consistent iteration (VQE#2) starts.
In addition to the energy evolution, we parameterize the wavefunction error using the infidelity, defined as 1 − F with the fidelity F = |⟨Ψ VQE |Ψ ED ⟩| 2 . In this context, the reference ground state |Ψ ED ⟩ is chosen as the optimal solution for each VQE and NGS step of the outer-loop iteration. Figure 1c shows the convergence of the wavefunction for the same system as Fig. 1b. By comparing these two panels, one can observe that a slow energy evolution may come with a relatively rapid change of wavefunction parameters, indicative of a bar-ren plateau [14]. Therefore, the infidelity may experience significant changes in later (outer-loop) iterations when the energy is close to convergence. We emphasize that the infidelity is posterior and cannot be used as the target function of the iteration.
In contrast to solving a Hubbard model, the NGS-VQE method involves a self-consistent outer loop between electrons and phonons. Thus, the combined NGS-VQE efficiency is essential for optimal results. To mitigate optimization issues of the variational quantum circuit, for instance the barren plateau or a multitude of local minima, we employ a three step optimization. First of all, since all gates within a single ansatz layer share the same variational parameter θ i , we can reuse parameters across circuits for different system sizes L. We therefore pre-run the VQE with smaller system sizes to initialize the circuit of the target system with these converged variational parameters. Moreover, the ground state evolves adiabatically for small changes in the model parameters (u, λ, and ω) within the same phase. This is why we further recycle converged parameters if results for similar model parameters exist. Finally, we adaptively adjust the VQE convergence criterion during the outer-loop NGS-VQE iterations by increasing the number of quantum measurements. Since the initial phonon state and variational non-Gaussian parameters are far from saddle points, we start using just a small number of measurements to give a low-accuracy estimation of the electronic ground state; with the progress of NGS-VQE iterations, we gradually raise the convergence criterion for the electronic state. All these strategies help to improve the overall performance and reduce the runtime of the hybrid NGS-VQE algorithm, especially relevant in hardware implementations (see Supplementary Note 1).
B. Simulating the correlated electron-phonon systems

Charge and spin phases in the Hubbard-Holstein model
The Hubbard-Holstein model and its extension set the stage to study the interplay of electronic correlations and EPCs in quantum materials. At half-filling and in 1D, this model results in a rich phase diagram, hosting an AFM, CDW, and a narrow intermediate phase [53][54][55][56][57][58][59][60]. To demonstrate the accuracy and efficiency of the NGS-VQE algorithm, we first restrict ourselves to the pristine Hubbard-Holstein model with g ′ = 0 and simulate the spin and charge structure factors of the ground state for different model parameters. The (static) spin structure factor is defined as while the (static) charge structure factor is defined as Phase diagram of the one-dimensional Hubbard-Holstein model. a, b Charge N (π) and spin S(π) structure factors for the ground state of a 6-site Hubbard-Holstein model, simulated by the non-Gaussian solver variational quantum eigensolver (NGS-VQE) algorithm as a function of u for fixed a λ = 1.5 and b λ = 3.5. The charge-density-wave (CDW) and antiferromagnetic (AFM) regimes are marked. The phonon frequency is set as ω = 10 and the quantum circuit depth is n = 9. c, d Distribution of the static c charge and d spin structure factors in the u − λ parameter plane, for the same conditions as a and b. The dashed line indicates the anti-adiabatic phase boundary u = 2λ between AFM and CDW. e Ground state infidelity and f (absolute-value) energy error with respect to ED results for the phase diagram in c and d.
Using the half-filled system as the benchmark platform in this paper, we focus on the nesting momentum q = π for both structure factors. In the regime where electronic interactions dominate (u ≫ λ), the spin structure factor S(π) prevails over the charge structure factor, reflecting an AFM state in a finite cluster (see Fig. 2a, b for ω = 10). With the increase of u − 2λ, N (π) gradually vanishes as charge degrees of freedom are frozen with a substantial energy penalty for double occupations. In the other limit where EPCs dominate (λ ≫ u), the charge structure factor N (π) prevails over the spin structure factor S(π). This reflects the onset of a CDW state, although a more rigorous identification requires either scaling to larger system sizes or excited-state analysis. Physically, the CDW is stabilized by the energy gain through a lattice distortion, forming an alternating pattern of holons and doublons. We summarize the dependence of both spin and charge structure factors on the two interaction parameters in Figs. 2c and d. The trends of these two observables reflect the two dominant phases qualitatively consistent with physical intuition. Due to the underlying finite-size system, the two phases are separated by a crossover instead of a sharp phase boundary. Recent studies have shown the presence of an intermediate Luther-Emery liquid phase for u ≈ 2λ, whose width is controlled by the phonon frequency [60]. The discussion of this phase requires finite-size scaling and is beyond the scope of this paper.
To determine the quantitative accuracy of our hybrid quantum algorithm, we compare the converged ground state for each set of model parameters against the NGSED solutions. The accuracy of the latter has been benchmarked by ED and DQMC [45,46,77]. As shown in Fig. 2e, the infidelity map suggests that even the largest error is in the single-digit percentage range, at most 0.03. These errors do not change significantly with increasing system size, as outlined in the "Scaling in circuit depth and system size" subsection. Interestingly, the most accurate solutions (with infidelity of order 10 −5 ) are obtained near the boundary of the CDW and AFM phases, i.e., along the diagonal u ≈ 2λ. In this regime, the finite-system solution is more metallic, due to the delicate balance between the electronic repulsion and phononmediated attraction. Therefore, the true ground state of systems near the phase boundary can be efficiently captured with a Slater determinant prepared by Givens rotations [73][74][75][76]. In contrast, the infidelity increases when the system evolves into CDW or AFM states, although the NGS-VQE algorithm yields quantitatively accurate results throughout the phase diagram. This observation is contrasting the intuition that the AFM or CDW states are more classical. Instead, these states are cat states in these small and low-dimensional systems. An accurate representation of these highly-entangled states with long-range correlations, therefore, requires deeper quantum circuits and, accordingly, more gates. The depen- Local (r = 0) and nearest-neighbor (r = 1) interaction strengthsṼ for various Hubbard u values within the antiferromagnetic (AFM) phase (g ′ = 0, λ = 1) and small phonon frequencies (close to the adiabatic limit). The inset in a highlights the considered path in the phase diagram (black arrow). c, d Local (r = 0) and nearest-neighbor (r = 1) interaction strengthsṼ for the same interaction parameters as a and b, but with large phonon frequencies. The dashed lines suggest the asymptotic values in the anti-adiabatic limit. e-h Local (r = 0) and nearest-neighbor (r = 1) interaction strengthsṼ for the same conditions as (a-d), but for the extended Hubbard model with g ′ = g/ √ 5. Similarly, the dashed lines in g and h indicate the asymptotic values in the anti-adiabatic limit. dence on circuit depth will be discussed in the "Scaling in circuit depth and system size" subsection. This sensitivity of the simulation accuracy to the model parameters is also reflected by the error of the ground-state energy, as shown in Fig. 2f.
Up to now, the benchmark has been conducted with relatively large phonon frequencies ω = 10, where the competition between CDW and AFM states is primarily controlled by the effective local interaction u eff = u − 2λ after integrating phonon fields. However, phonon frequencies in typical correlated materials are usually comparable to the electronic bandwidth, if not even reaching the adiabatic limit (ω → 0). The dependence of charge and structure factors on different phonon frequencies is discussed in Supplementary Note 2. In the thermodynamic limit, smaller phonon frequencies usually lead to a steeper crossover between the two phases [54][55][56][57], with both S(π) and N (π) dropping more rapidly when approaching the phase boundary. Note, however, that this intermediate phase cannot be resolved in a small cluster.

Phonon-mediated interactions in the
Hubbard-extended-Holstein model The wavefunction ansatz in Eq. (2) allows to extract the effective model H eff in the polaronic basis while solving for the ground state. This approach has been used to quantify the recently discovered nearest-neighbor electronic attractionṼ in cuprate chains [50]. Here, we eval-uateṼ for different systems using the NGS-VQE algorithm both to provide intuition about phonon-mediated interactions in different limits and to benchmark the validity of the method under different conditions.
We first examine the Hubbard-Holstein model without g ′ . Figures 3a,c and b,d shows the simulated local and nearest-neighbor attractive interaction in the polaronic basis. Since this interaction was discovered in cuprates with strong electron correlations, we restrict ourselves to the AFM regime (u ≫ λ) and set λ = 1. In the antiadiabatic limit, the local interactionṼ ii asymptotically approaches −λ = −1 and all nonlocal interactionsṼ i̸ =j vanish (see Fig. 3c, d). This is consistent with the integration of phonons in field theory. With the decrease of ω and proximity to u, Coulomb interactions start to influence the distribution ofṼ ij . Such an influence is more obvious for large u, where the electronic correlations are strong and renormalize the phonon self-energy. This effect can be captured by the wavefunction ansatz of Eq. (8). The electron-dressing effect for the phonon selfenergy is described by Ω eff in ref. [45]. Simultaneously, the retardation effect of finite-frequency phonons mediates the effective interaction at finite distance. Therefore, the effective nearest-neighbor attractionṼ i,i+1 increases for lower phonon frequencies (see Fig. 3b, d).
With further decreasing phonon frequencies, the im- pact of phonons becomes a mean-field-like deformation potential instead of a virtual process. Such a deformation potential contributes as an overall chemical potential instead of media of a two-particle interaction. Therefore, the simulatedṼ ij becomes extremely nonlocal and evolves into an all-to-all interaction, which is equivalent to a chemical potential shift in a canonical ensemble with fixed particle number. In the ω → 0 limit, phononmediated interactions are insensitive to the Hubbard u, because the mean-field phonon distortion is determined only by the average local electron density, which is approximately fixed in the AFM phase (see Fig. 3a, b). This insensitivity is similar to the anti-adiabatic limit but has a different origin.
We now move on to discuss the impact of nearestneighbor EPC g ′ on the complexity and accuracy of the hybrid quantum algorithm. We fix the ratio g ′ /g = 1/ √ 5 to reflect the geometric relation of apical oxygens in a 1D transition-metal oxide [50]. With electrons coupled to their nearest-neighbor phonons directly, the mediation ofṼ i,i+1 can be generated without retardation effect. Therefore, the asymptoticṼ i,i+1 in the anti-adiabatic limit is no longer zero but acquires a finite value (see Fig. 3h). This asymptotic interaction strength can be analytically calculated as −2gg ′ /ω, which is more evident than the additional interactions caused by retardation effects (the latter is three orders of magnitude smaller than the former). At the same time, the phononmediated local interactionṼ ii is further strengthened by this nonlocal g ′ coupling. The anti-adiabatic value ofṼ ii approaches −g 2 /ω − 2g ′2 /ω (see Fig. 3g). Due to the geometric distance controlling g ′ /g, the strength ofṼ ii is still comparable to that of the Hubbard-Holstein model (g ′ = 0). In the adiabatic limit, the effective interactions are similar to those obtained from the Hubbard-Holstein model, asymptotically approaching an extremely delo-calizedṼ ij . Compared to the g ′ = 0 case (Figs. 3a and b), the only difference is the asymptotic value when ω approaches zero. Using the fact that, in the longwavelength limit (q = 0),Ṽ is proportional to g 2 q /ω and g q = g+2g ′ cos q, we can estimate the ratio between these two asymptoticṼ ij values (for g ′ = g/ √ 5 and g ′ = 0) to be (1 + 2/ √ 5) 2 ≈ 3.59. This ratio agrees with the simulated results in Figs. 3a, b, e, and f.

C. Scaling in circuit depth and system size
The results presented in Fig. 3 demonstrate that our hybrid quantum algorithm is able to produce quantitatively accurate results across the full phonon spectrum (see Supplementary Note 2 for a comparison to NGSED results). To further analyze the accuracy of our algorithm, we investigate the influence of different system sizes L. The depth of the quantum circuit controls the expressibility of the variational state, potentially allowing for a more accurate approximation of the electronic ground state by increasing n. However, an efficient encoding on quantum computing platforms is only possible if the depth of the employed quantum circuit does not scale exponentially with the system size L. As shown in Fig. 4a, the ground states for small-u systems can be efficiently expressed by a Slater determinant. Thus, only a few layers are needed to reach ground state energy errors below 1 × 10 −6 when compared to ED. Larger u, however, requires deeper circuits, reaching a plateau of errors of the order 10 −4 to 10 −5 in the ground state energy. In order to investigate the scaling of necessary n with the system size L, we consider a fixed error in the electron ground state energy of |E VQE − E ED |/L = 0.1t. The circuit depth required to achieve this performance as a function of L is shown in Fig. 4b, highlighting a moderate increase in depth for small and large u. Intermediate values for u, however, require significantly deeper circuits, as quantum fluctuations are larger around the crossover between metallic and AFM phase.
So far, we have considered the influence of finite-depth quantum circuits on the electronic part of the many-body ground state. The hybrid NGS-VQE algorithm combines the electron solver with a variational NGS solver for the phonon and non-Gaussian components. That being said, errors in the VQE solutions do not necessarily accumulate, but can actually be mitigated by the phonon solver. Considering again a fixed error in the electron groundstate energy of |E NGS−VQE − E ED |/L = 0.1t and corresponding circuit depth n, we investigate the accuracy of the hybrid quantum algorithm in different regions of the phase diagram. Specifically, we consider the relative error in the ground-state energy of the converged extended-Hubbard Hamiltonian [see Eq. (9)], containing the phonon dressing of kinetic hopping and long-range interactions. Figures 4c, d indicate that the NGS-VQE simulation errors are usually at least an order of magnitude smaller than those of the electronic solvers. The largest errors are obtained for small phonon frequencies (Fig. 4d, ω = 1), where the absence of quantum fluctuations hinders the phonon solver to escape local minima during the self-consistent iteration [45]. Warm-up iterations with larger phonon frequencies can help to alleviate this issue [45]. Moreover, the relative errors do not increase for systems larger than L = 4, indicating quantitatively accurate results across different system sizes and phases. The only exception appears in the CDW phase at small phonon frequencies such as ω = 1, where the relative error oscillates with L, likely due to the degeneracy of ground states. The ability to mitigate errors of the quantum solver also provides a promising path to experimental realizations. Hardware implementations, irrespective of the specific platform, suffer from decoherence [78], making noise resilient algorithms of key importance. Our hybrid quantum algorithm is able to improve VQE results over a wide range of phonon frequencies, phase regions, and noise levels, suggesting efficient hardware realizations (see Supplementary Note 3).

III. CONCLUSIONS
Our NGS-VQE method provides a general framework for performing accurate and efficient quantum simulations of electron-phonon systems with arbitrary interaction strengths. Using this method, we have studied the Hubbard-Holstein and HEH models as examples, reproduced the CDW-AFM crossover with high precision, and extracted the phonon-mediated interactions in a wide range of phonon frequencies. While we focused on paradigmatic (and experimentally relevant) cases, this method can be generally applied to any model with electronic Coulomb correlations and Fröhlich-type electronphonon couplings. The commutation between the NGS transformation [e iS with S defined in the form of Eq. (7) of the "Methods"] guarantees a closed-form effective Hamiltonian similar to Eq. (9). This hybrid algorithm can be extended to other types of electron-boson interactions (like the Su-Schrieffer-Heeger phonon and cavity QED) through the generalization of the non-Gaussian transformation U NGS and its optimization strategy. Anharmonic potentials can be tackled at the price of replacing |ψ ph ⟩ by more complicated many-body wavefunctions similar to the electronic ones. Both generalizations are accompanied by the increase of computational complexity and should be designed based on the requirements of specific models. Moreover, as demonstrated in Ref. [46], this framework can be extended to non-equilibrium dynamics, which requires a reliable quantum solver for the long-time propagation of |ψ e ⟩. The Fourier transform of non-equilibrium dynamics with two-or multi-time correlation functions further paves the way to excited-state spectra [79][80][81].
Our work shows that the phonon solver can mitigate potential errors of the quantum hardware, facilitating a future experimental implementation. The success of an experimental realization, however, relies on an efficient implementation of the required quantum circuits, respecting the connectivity of a given device. This is especially the case for nonlocal gates, like the P gates of our ansatz representing electron-electron interactions. Trapped ion and cold atom based platforms offer high qubit-connectivities [82], reducing potential swap-overheads when implementing entangling gatelayers such as the P and H gates in our ansatz. Superconducting systems, on the other hand, have a more limited qubit-connectivity but operate at much faster rates, making them favorable for two-level NGS-VQE iterative schemes. Therefore, practical implementations of the algorithm proposed in this study call for developing higher-connectivity superconducting hardware or error mitigation schemes to compensate for potential swapoverheads.

A. NGS-VQE method and effective Hamiltonian
As mentioned in the main text, the variational electron-phonon wavefunction in the NGS-VQE method is given as with the non-Gaussian transformation U NGS = e iS . As benchmarked by ED and DQMC on small clusters [45,46], it is sufficient to truncate the S operator to the lowest-order terms where we use the momentum-space electron density ρ q = i,σ n i,σ e −iqxi , and the phonon momentum operator The goal of the NGS-VQE solver is to minimize the total energy in Eq. (3) in the variational parameter space spanned by {f q }, |ψ ph ⟩, and |ψ e ⟩. Without considering anharmonicity, the phonon state to the right of U NGS should be weakly entangled and can be efficiently captured by variational Gaussian states Here, ∆ R , ξ q are variational parameters and R q = (x q , p q ) T denotes the bosonic quadrature notation with canonical position x q and momentum p q , where we adopt the reciprocal representation for the phonon displacement For convenience, we parameterize the phonon state using the linearization of U GS named S q , which satisfies U † GS (x q , p q ) T U GS = S q (x q , p q ) T . The NGS-VQE method minimizes the total energy by updating |ψ e ⟩ and |ψ ph ⟩ iteratively.
With fixed U NGS and |ψ ph ⟩, the electronic problem that the quantum machine has to solve is the ground state of an effective Hamiltonian whereẼ ph = 1 4 ω q Tr[S q S † q ] − 2 . The phonondressed hopping becomes and the effective interaction is The VQE solution of the effective Hamiltonian in Eq. (9) gives |ψ e ⟩ in Eq. (6). The iterative optimization of U NGS and |ψ ph ⟩, for fixed |ψ e ⟩, follows the imaginary time evolution in Ref. 45. It is worth noting that the charge density correlation functions ⟨ρ q ρ −q ⟩ ∝ ij σσ ′ n i,σ n j,σ ′ necessary for the imaginary time evolution of |ψ ph ⟩ appear in H eff as well. They are therefore already measured with the energy expectation value during VQE and result in no additional computational cost.

B. Quantum circuit and ansatz
To represent the effective Hamiltonian in Eq. (9) on a quantum computer, we rely on the Jordan-Wigner transformation: each electron with given spin orientation is mapped to one qubit. Specifically, it reads where the phase factors retain the fermionic anticommutation in the spin operators S. Transforming the effective model described by Eq. (9) with L sites then yields 2L spin operators. Consequently, a quantum ansatz for the electronic ground state contains 2L qubits, which represent occupied |1⟩ = c † |0⟩ or unoccupied |0⟩ fermionic states. As the effective Hamiltonian preserves occupation number and total spin, we can restrict the electron occupation to half-filling and total spin to zero. Correspondingly, we arrange the qubits representing the two spin orientations separately (see Fig. 1a), and require that gates connecting the two spin sectors cannot change their respective occupation. As outlined in the main text, the employed quantum circuit is based on the Hamiltonian variational ansatz [72]. To encode the non-interacting (U = 0, V ij = 0) electronic model, a sequence of Givens rotations G(θ), parametrized by θ, is applied to adjacent qubits [73][74][75][76].
In the basis |00⟩, |01⟩, |10⟩, |11⟩, the gate is defined as The ground state of the full effective model is then obtained by an adiabatic evolution with the Hubbard-like Hamiltonian [15,17]. It can be decomposed into kinetic hopping terms, in the basis |00⟩, |01⟩, |10⟩, |11⟩ represented as and on-site interactions The alternating sequence of phase gates (P ) and hopping gates (H) is then repeated for a number of repetitions n, controlling the expressibility of the variational ansatz.

DATA AVAILABILITY
The presented data are deposited into the public folder Figshare. Additional numerical data that support the findings of this study are available from the corresponding authors upon reasonable request.

CODE AVAILABILITY
The relevant scripts of this study are available from the corresponding authors upon reasonable request. Hardware implementations are usually limited by the decoherence time of the employed qubit platform. Consequently, limiting the depth of a quantum circuit is crucial, however, this results in less accurate VQE results when solving for the ground state of H eff . Considering the outer-loop NGS-VQE iterations, we propose an adaptive convergence scheme to circumvent this issue: starting from a lower accuracy of the ground-state energy obtained by VQE, we gradually increase the accuracy with the progress of the outer-loop iterations between NGS and VQE. A lower accuracy in the estimated ground state energy corresponds to shallower quantum circuits, which are then gradually deepened to increase the expressibility. Figure S1 shows the comparison between a fixed-accuracy convergence criterion (FCC) and adaptive convergence criterion (ACC) for various model parameters. Here, we start with a VQE accuracy of ∆ϵ = 10 −1 , i.e. we terminate the VQE (inner-loop) iterations if the energy between iterations changes by less than ∆ϵ, i.e., |E VQE i+1 − E VQE i | < ∆ϵ. The outer-loop NGS step is triggered using this less accurate electronic state obtained from VQE. Subsequently we increase ∆ϵ with each NGS-VQE iteration, until we reach numerical precision at ∆ϵ = 10 −9 . The performance of this ACC strategy depends on the phase of the ground state. In an AFM (u = 10, Fig. S1a) or intermediate metallic phase (u = 5, Fig. S1b), the ultimate electronic ground state is not significantly different from the initial guess (Néel state). Therefore, the ACC strategy perfectly reproduces the exact solutions, while consuming much less iteration steps to reach the same total-energy accuracy compared to the regular FCC strategy. In contrast, the ACC strategy fails to converge to the correct ground state in the CDW phase (u = 0, Fig. S1c). This is because the ground-state electronic configuration in this case has alternative double occupations, which is intrinsically different from the initial electronic state. Therefore, a less accurate VQE solution, in the first few outer-loop iterations, fails to drive the electronic states to a configuration similar to the ground state; this failure further delays the relaxation of phonon states to form alternating distortions and, accordingly, fails to converge to the true ground state within a reasonable number of steps.
In general, the ACC requires more outer-loop NGS-VQE iterations, as the phonon states have to relax to account for the loss in precision. Nevertheless, except for the CDW phase, the ACC strategy requires less VQE iterations in total, with the same accuracy of the electron-phonon ground states. In other words, errors in the electronic part of FIG. S1. Convergence comparison of VQE approaches. We compare a VQE algorithm with adaptive energy convergence criterion (ACC, bottom), starting at small accuracy for the first few iterations (∆ϵ = 10 −1 , 10 −1 , 10 −2 , 10 −2 , 10 −3 , 10 −3 , then 10 −9 ), with a VQE version of fixed accuracy of ∆ϵ = 10 −9 (FCC, top). Inner-loop VQE iterations are displayed in full detail, whereas NGS iterations are compressed. (a,b) For large interactions u = 10, 5, adaptive convergence speeds up the groundstate search and even improves the overall result when compared to ED. c For vanishing interaction strengths u, convergence is slower and worse for the adaptive approach. The phonon solver cannot account for the loss of accuracy in the electronic wavefunction. Consequently, at finite u, noisy or less accurate VQE outputs can still lead to precise results, an important criterion for hardware implementations. All panels are for ω = 10, λ = 2, and L = 4 systems and quantum circuit depth n = 6. the hybrid quantum algorithm can be mitigated by the phonon solver, which is conducted on a classical computer and is computationally cheaper.

SUPPLEMENTARY NOTE 2: IMPACT OF THE PHONON FREQUENCY
The anti-adiabatic limit u = 2λ boundary controls the energetically favorable single-site electronic configuration and can be used to estimate the transition between CDW and AFM states. However, simulations with larger system sizes have demonstrated that an intermediate Luther-Emery liquid phase emerges near the boundary, where these two insulating instabilities balance [57,59,60]. This intermediate state is driven by the electronic hopping t but is also controlled by the phonon frequencies. Fig. S2a shows the parameter dependence of S(π) and N (π) on the phonon frequency. With increasing phonon frequency, the two correlation functions slightly separate from each other, reflecting the underlying intermediate phase. This change is not obvious in the small 6-site system. As a comparison, we further present the simulation results obtained from larger systems (see Fig. S2). These simulations are conducted using the NGSED instead, due to the limited gates in our NGS-VQE setting. With the increase of system size, the impact of phonon retardation effects at finite frequencies becomes more obvious. The intermediate regime expands with the phonon frequency and suppresses both the charge and spin instabilities in a relatively wide range of parameters.
We further discuss the impact of phonon frequencies on the effective interaction in the Hubbard-Holstein model. As discussed in the Fig. 3 of the main text,Ṽ ij approaches −g 2 δ ij /ω in the anti-adiabatic limit (ω → ∞), while it becomes distributed into a uniform all-to-all interaction in the adiabatic limit (ω → 0). Here, we show the entire ω dependence for both the local and nearest-neighbor interactions (see Fig. S3a,b). The evolution from adiabatic to anti-adiabatic limits are continuous, without any obvious phase transitions. (Note that the x-axis changes its scale at ω = 1, which leads to the artificial "kink" in the figure.) As mentioned in the main text, the dependence on u becomes less relevant in both limits due to different origins.
In addition, we analyze the relative error of the NGS-VQE simulatedṼ ij against NGSED results. As shown in Fig. S3c-f, largest relative errors for both the local and nearest-neighbor interactions are obtained for a smaller Hubbard interaction u close to the phase transition. This error decreases with increasing u. We attribute this error to the sensitivity of the effective interactionṼ ij to the electronic state. When the system is deeply in the AFM phase, electron density fluctuations are heavily suppressed, although the magnetic configuration stays in a cat state. Since phonons couple to the electron density, the wavefunction is determined to lowest order by the density distribution. Therefore, the phonon-mediated interactionsṼ ij exhibit very small errors for large-u systems. In contrast, the u = 3 system (with λ = 1) is close to the crossover, where the electron density fluctuation is still strong. Such a fluctuation requires deeper quantum circuits to capture, especially in the shallow AFM side where most gates are used to express the spin configurations. For any strengths of the Hubbard u, the relative error for the on-site interactionsṼ ii decreases with the rise of ω, sinceṼ ij narrows in space and the absolute value of the on-site interaction increases. In contrast, the relative error for the nearest-neighbor interactionṼ i,i+1 increases with ω, which originates from the same reason (see Fig. S3e,f ).

SUPPLEMENTARY NOTE 3: NOISE SIMULATIONS
To assess the practicality of our approach with near-term quantum devices, we study how resilient it is to characteristic hardware noise. To this end, we include both statistical noise (10 5 shots per operator expectation value) and realistic device noise in our simulations, specifically, the noise model of IBM's device ibmq kolkata [83]. Scaling the average device errors across several orders of magnitude allows for a detailed understanding of their influence on the solver's accuracy. In particular, our model applies the following noise to all qubits, which are common values for the ibmq kolkata device at η = 1: T 1 = T ave 1 /η , T 2 = T ave 2 /η , e 1q = ηe ave 1q e 2q = ηe ave 2q , e ro = ηe ave ro (19) where η ∈ {0.1, 0.5, 1, 5, 10} is a scaling factor and T ave 1 = 106.1 µs, T ave 2 = 82.93 µs, e ave 1q = 3.78 × 10 −4 , e ave 2q = 1.07 × 10 −2 , and e ave ro = 2.27 × 10 −2 are the device's average relaxation time, dephasing time, one-qubit gate, twoqubit gate, and readout error, respectively.
We observe the expected behavior of an increasing simulation error with increasing levels of device noise in the AFM phase and at the phase transition (see Fig. S4b). Crucially, we obtain ∆ rel (E) = |E NGS−VQE − E NGSED |/|E NGSED | at η = 1 about one order of magnitude larger than in the statevector case (see Fig. 4 of the main text) and dropping off steeply for lower noise levels, as are anticipated with the next generation of quantum devices. For the CDW phase, however, we see a slight decrease in simulation error with increasing hardware noise. This can occur when the optimal solution is such that the noise will naturally relax the system toward this solution, e.g., the qubit ground state |0 . . . 0⟩. Additionally, as the CDW phase is dominated by electron-phonon interactions, the NGS solver becomes more important. Fluctuations have proven to be crucial to escape local minima [45], possibly also explaining the increasing accuracy of the simulation results with noise.
FIG. S3. Converged extended Hubbard interactions as a function of ω for u = 3, 6, 9 with fixed λ = 1. a On-site (r = 0) interaction strengthṼ for three u values within the AFM phase (g ′ = 0). While being insensitive to the strength of electronic correlations u for small phonon frequencies,Ṽ asymptotically approaches −λ = −1 in the anti-adiabatic limit (ω → ∞). b Nearest neighbor (r = 1) interaction strength for three u values within the AFM phase (g ′ = 0). Starting from finite interactions at small phonon frequencies,Ṽ quickly reduces in magnitude and vanishes in the anti-adiabatic limit (ω → ∞). c,e Absolute and relative error of the on-site (r = 0) interaction strengthṼ shown in a as compared to NGSED (g ′ = 0). d,f Absolute and relative error of the nearest-neighbor (r = 1) interaction strengthṼ shown in b as compared to NGSED (g ′ = 0).