Zero-Temperature Phases of the 2D Hubbard-Holstein Model: A Non-Gaussian Exact Diagonalization Study

We propose a novel numerical method which embeds the variational non-Gaussian wavefunction approach with exact diagonalization, allowing for efficient treatment of correlated systems with both electron-electron and electron-phonon interactions. Using a generalized polaron transformation, we construct a variational wavefunction that minimizes the entanglement between electrons and phonons; exact diagonalization is then used to treat the electronic part of the wavefunction exactly, thus taking into account high-order correlation effects beyond the Gaussian level. Keeping the full electronic Hilbert space, the complexity is increased only by a polynomial factor relative to the exact diagonalization calculation for pure electrons. As an example, we use this method to study ground-state properties of the two-dimensional Hubbard-Holstein model, providing evidence for the existence of intervening phases between the spin and charge-ordered states. In particular, we find one of the intervening phases has strong charge susceptibility and binding energy, but is distinct from a CDW-ordered state; while the other intervening phase displays superconductivity at weak couplings. This method, as a general framework, can be extended to treat excited states and dynamics, as well as a wide range of systems with both electron-electron and electron-boson interactions.


I. INTRODUCTION
Strongly correlated systems have attracted much attention in the past few decades, not only because of the unconventional physics emergent from these systems, but also because they pose important theoretical questions about the nature of interacting systems at intermediate and strong coupling. Away from weak coupling, traditional mean-field or perturbative approaches often fail to accurately describe the physics, especially in cases with competing and/or intertwined ordering tendencies.
In the condensed matter setting, models can be classified into interacting electrons, interacting bosons, and interacting electron-boson systems. Advances in unbiased numerical many-body methods, including exact diagonalization (ED) 1,2 , quantum Monte Carlo (QMC) 3,4 and density-matrix renormalization group (DMRG) 5,6 , have greatly expanded our understanding of the fermionic and bosonic Hubbard models, together with their variants. For example, recent exact solutions of the single-and three-band Hubbard models have shed light on the stripe and d-wave superconducting phases in doped cuprates 7,8 . Although the Hubbard model is often a prototype, experimental evidence suggests coupling to phonons can also play an important role in the low-energy physics of correlated materials. For example, STM measurements have shown a significant isotope effect on the secondderivative tunnelling current 9 ; spectral experiments have shown significant lattice effects in cuprates, starting from the underdoped regime 10 , to optimal 11 and overdoped regimes 12 ; phonon softening has also been observed using Raman 13 and neutron scattering 14 . These observations suggest that electron-electron (e-e) and electron-phonon (e-ph) interactions should be taken into account simultaneously in order to properly understand the rich phenomena observed in many correlated materials.
A significant barrier to understanding the low-energy physics of models with both e-e and e-ph interactions is the challenge they pose to conventional numerical methods. On the one hand, numerical many-body approaches, such as ED and DMRG, have achieved great success in analyzing correlated electronic systems in the past decades. With the improvement of both algorithms and high-performance supercomputers, these approaches not only evaluated the ground state properties precisely, but also calculated the spectroscopies and dynamics in a well-controlled way 15,16 . However, when also taking phonons into consideration, efficient numerical treatment becomes a problem. The bosonic Hilbert space is infinitedimensional, and the total allowed phonon number has to be truncated to a small value (on the order of 1-3 phonons per site). This has largely limited the study of strongly coupled e-ph systems, e.g. the Peierls CDW systems.
On the other hand, approximate methods based on variational wavefunctions provide an alternative route to analyze correlated systems. For example, variational Lang-Firsov transformations have been applied to disentangle e-ph systems in the long-wavelength limit 17,18 . A more intricate Jastrow variational wavefunction has been employed to examine the competing spin, charge, and superconducting orders via a particular mean-field decoupling of the electrons. 19,20 More recently, these variational approaches were generalized to the non-Gaussian class of wavefunctions. With the non-Gaussian transformation chosen to be a generalized polaron transformation, this method gives a good estimate of the eph ground state 21 ; with specific parity transformation, this class of wavefunctions also perfectly decouples the Kondo and Anderson models [22][23][24] . However, extending the method to systems with e-e interactions has not been straightforward, due to the fact that four-fermion interaction terms make the parameter space much more complicated. Besides, effective e-e interactions can be generated when disentangling the e-ph coupling. The absence of correlations in the Gaussian state limits the accuracy for even pure e-ph systems. This issue becomes even more crucial for the calculation of dynamics, due to greater complexity of the polaronic dressing 25 . More precise treatment of electronic correlations is therefore imperative.
To combine the merits of these two philosophies, we propose the non-Gaussian exact diagonalization (NGSED) method. Due to the polaronic non-Gaussian ansatz for the phonon part of the wavefunction (to be described in more detail below), the complexity of the ED-based electronic calculation is increased only by a polynomial factor. At the same time, the inclusion of the full electronic Hilbert space and many-body wavefunction addresses the correlation issue of pure variational approaches, reducing the bias incurred by a mean-field treatment of the electronic state. Similar embedding ideas have been attempted in a few numerical studies. For example, the iterative optimized phonon implementation has been applied to ED 26 and CPT 27 . However, even on an optimized basis, the phonon number still spans a huge Hilbert space, limiting the calculations to a 6-site chain. The classical phonon approximation 28 and the standard Lang-Firsov transformation 29 were also employed to disentangle the local interactions in QMC and ED, but ignorance of the explicit phonon wavefunction prevents an accurate description of fluctuations of both effective tunneling and interactions. A very similar idea of embedding variational Lang-Firsov transformations with ED has been attempted in the t−J model 30,31 . With a long-wavelength, homogeneous coupling parameter, these embedding calculations still failed to capture the fluctuations caused by the polaronic dressing. Therefore, a natural extension is the embedding of a variational phonon wavefunction and polaron transformation into an exact numerical technique -this forms the intuition of our NGSED method.
Though the idea of embedding non-Gaussian transformations with numerical many-body techniques can be extended to a variety of problems, we focus on the e-ph system as a concrete topic in this article. We introduce the NGSED method for a generic e-ph model and present the iterative approach to evaluate the groundstate properties. To assess the accuracy of the variational wavefunction, we benchmark the method against an exact QMC solution of the Holstein model. We then focus on the Hubbard-Holstein model, where we examine the ground-state properties and their dependence on the ee interaction, e-ph coupling, phonon energy, and doping. We observe a shift in the antiferromagnetic (AFM) phase boundary, which we explain using the form of effective e-e interactions. We identify a region between the AFM and charge-density-wave (CDW) states in which both charge and spin orders are absent. This region can further be divided: one of the subregions has enhanced charge susceptibility and considerable binding energy, possibly corresponding to a 2D analog of the Luther-Emery liquid observed in the 1D Hubbard-Holstein model 32 ; the other subregion exhibits superconductivity at the weak-coupling side but gradually becomes metallic for stronger coupling. In contrast to the intermediate phases obtained using pure variational wavefuctions 19,20 , we do not see a dramatic broadening of the superconducting phase on the weak-coupling side, consistent with unbiased QMC conclusions 33 . Complementing previous high-temperature QMC studies, truncated-phonon ED studies, and zero-temperature variational studies, this work sheds new light on the phases in such a competingorder system.
The organization of this article is as follows. We first introduce the NGSED method and relevant derivations in Sec. II. Then we apply it to the Holstein and Hubbard-Holstein models and discuss the ground-state properties in Sec. III. We conclude our method and simulations in Sec. IV, together with the outlook of this method for other systems.

II. MODEL AND DERIVATIONS
We present the derivation of relevant formulae for a generic electron-phonon system in this section, before focusing on the Hubbard-Holstein model with specific form of e-e and e-ph interactions. A generic electronphonon model can be expressed by the Hamiltonian where c kσ (c † kσ ) annihilates (creates) an electron at momentum k with spin σ, with a dispersion relation ε k and chemical potential µ. N is the overall site number. Within second quantization, c kσ takes the reciprocal representation with respect to the annihilation operator of the Wannier orbital c kσ = i e −ik·r i c iσ / √ N . Apart from the bare dispersion, the H e−e , H e−ph and H ph terms represent the contributions from e-e interactions, e-ph coupling and phonon energy, respectively.
In general, the phonon part of Hamiltonian is and the e-ph coupling part is 34 Here, the ω q describes the phonon dispersion, g q parametrizes the e-ph interaction at a wavevector q; a q annihilates a phonon at momentum q and ρ q = iσ n iσ e −iq·r i is the electron density. For convenience in subsequent derivations, we employ the bosonic quadrature notation R q = (x q , p q ) T , with the canonical position x q = a q + a † −q and momentum p q = i(a † −q − a q ) determined by the phonon annihilation operator. These canonical operators fulfill the commutation relations Thus, the parts of Hamiltonian relevant to phonons can be rewritten as with e 1 = (1, 0) T . Without loss of generality, we allow the parameters g q and ω q to vary over momentum space, but obeying the time-reversal symmetry i.e. g q = g −q and ω q = ω −q . For the electron interaction part H e−e , the only restriction we place is that it commutes with the local density operators n i . Thus we allow for any combination of density and spin operators, such as the on-site Hubbard or long-range Coulomb interactions.
To describe the e-ph entangled system in the simplest form, we consider the wavefunction ansatz Here, the right-hand-side is a direct product of electron and phonon states, where |ψ e is treated as a full manybody state while |ψ ph is treated as a coherent Gaussian state The polaron transformation U plrn creates entanglement between these two parts of the wavefunction: In the above wavefunction prototype, the ∆ R , ξ q and λ q are variational parameters. An important feature of the wavefunction ansatz in Eq. (6) is that this wavefunction gives exact solutions to the e-ph problem in both the adiabatic (ω = 0) and anti-andiabatic (ω = ∞) limits.
In the adiabatic limit, phonons can be treated as a classical field, mean-field theory becomes exact, and the Gaussian wavefunction gives an exact description of the state. In the anti-adiabatic limit, the phonon field can be intergrated out, yielding an instantaneous, attractive on-site interaction; i.e. the attractive Hubbard model. The exact diagonalization step then solves this problem exactly. As it is exact in both the adiabatic and anti-adiabatic limits, we expect Eq. (6) does not induce significant bias in realistic models with finite ω. The accuracy of this assumption will be further assessed through the comparison with exact DQMC solutions [see Sec. III A]. Note, in principle the coherent part of the phonon wavefunction |ψ ph can involve displacements for all different momenta. However, any finite value of finite-q displacement would lead to the manual breaking of translational symmetry and over-estimate the tendency of ordering. Therefore, to avoid possible biases induced by the symmetry-breaking Gaussian states, we neglect any finite-q displacements in Eq. (7). Physically, it means phonons cannot really condense at a finite momentum, though the system might exhibit a strong instability. We impose this strong assumption because neither of magnons or Cooper pairs can condense in such a smaller cluster. Therefore, to fairly study the competition between spin-and charge-density-wave states, we only discuss their susceptibilities instead of long-range order parameters. This assumption also highly reduces the Hilbert space dimension due to the momentum conservation.
The above polaron transformation generalizes the Lang-Firsov transformation 35 . Historically, the Lang-Firsov transformation has been widely exploited in electron-boson systems to disentangle the coupling and simplify the calculation. To tackle the Hubbard-Holstein model, early attempts have pushed it to a variational transformation 17,18 . These transformations have shown advantages in solving the Holstein model 36 , Hubbard-Holstein model 37,38 , Anderson-Holstein model 39,40 , and anharmonic phonons 41 . However, due to the limitation of the numerical treatment on either the phonon or electronic side, these variational transformations were restricted only to a q−independent λ q . This treatment ignores the longer-range spatial fluctuation of the effective interaction mediated by the phonon, which we will show plays a significant role near the quantum phase transition. A direct consequence of this simplification is the overestimation of the CDW instability [we will further discuss this in Sec. III]. This limitation necessitates the generalization of this transformation to a polaronic non-Gaussian transformation.
By constructing the wavefunction through Eq. (6), we can evaluate the ground state with the manifold spanned by the variational parameters and the manybody electronic wavefunctions. Variational parameters are determined by minimizing the energy Numerically, the optimization can be iteratively achieved by decomposing into the electronic and bosonic state, with coupled coefficients. Each of them can be treated as a correction to the effective Hamiltonian while optimizing the other. Thus, for an equilibrium state, we minimize the total energy along two gradient directions sequentially. With an initial guess not far from the global minimum, we expect the many-body electronic state and variational phonon state to converge to the ground state self-consistently. In the following two subsections, we describe the procedures for evolving these two parts of the state. Afterward, we describe the above self-consistent iteration in a more strict manner using notations introduced in these two subsections.

A. Electron Ground State: Exact Diagonalization
We first optimize the electronic state (minimizing the energy), keeping fixed the variational parameters in U plrn and U GS . Then, Eq. (9) becomes an unrestricted minimization of energy in the full electronic Hilbert space, where the effective electronic Hamiltonian is given by tracing over the phonon state The H eff is an operator only on the electronic Hilbert space. Since the phonon state is Gaussian, the expression for H eff can be obtained analytically: Here the variational parameters for the phonon state are rewritten as Γ q = S q S † q , with the 2×2 matrix S q represents the linearization of the U GS , i.e., The polaronic dressing is reflected in the effective kinetic energy, i.e., the renormalized band dispersionε and the effective electronic attraction In the above derivations, we have employed the assumption that [H e−e , n i ] = 0. Note that in the last step of Eq. (12), the electron density at momentum q = 0 is nothing but the total occupation N e in a microcanonical ensemble. Therefore, the energy minimization with respect to ∆ R can be done immediately, leading As will be shown later in Eq. (24), λ 0 = g 0 /ω 0 for the saddle-point solution. Therefore, for the purpose of calculating the ground state, it is convenient to set ∆ R ≡ 0.
Different from the original Lang-Firsov transformation in the atomic limit, both the kinetic and interaction energies in the effective electronic Hamiltonian are renormalized by the phonons. The variational parameters allow us to find a balance between these two effects and minimize the entanglement between electrons and phonons 29 . Moreover, different from the widely used modified Lang-Firsov transformations 17,18 , the generalized polaron transformation and phonon Gaussian state naturally give momentum fluctuations of the effective interaction V q . In later discussions, we will show that this fluctuation is crucial near the phase boundary.
Since we keep the full electronic Hilbert space, it is straightforward to diagonalize the matrix H eff and find the ground state through a standard Lanczos approach. As we will discuss below, the ground state can be obtained alternatively through a flow equation -imaginary time evolution. However, with the full Hilbert space information, computing a matrix diagonalization is much cheaper than performing a time evolution, though the latter has been widely used in variational approaches.

B. Phonon Ground State: Imaginary Time Evolution
Keeping the electronic wavefunction fixed, the energy miminization in Eq. (9) can be achieved through the imaginary time evolution This evolution guarantees the monotonic decrease of average energy while maintaining the normalization of the wavefunction. If we restrict ∆ R = 0 as mentioned above, the derivative of the variational wavefunction becomes Taking into account the orthogonality of the electronic wavefunction basis, the tangential vectors are a † q a † −q |0 ph ⊗|ψ e and a † q |0 ph ⊗ ρ q |ψ e . 42 The rotated Hamiltonian is where α = x, y denotes the direction, while δ α is a unit vector along the α-direction. To determine the evolution of the variational wavefunction, we project the rotated Hamiltonian to the above two sets of tangential vectors. On the one hand, the projection with respect to the second-order bosonic terms is 21,25 ( The renormalized phonon energy matrix is Here E 22 = e 2 e T 2 and e 2 = (0, 1) T . Transforming the scalar equation of motion Eq. (18) into a matrix form, we should fill up the missing matrix elements in an antihermitian way, which gives Absorbing the gauge freedom, we have On the other hand, the projection to the other tangential vector gives We define the modulated electronic correlation function where both Π q and Θ q are real-valued functions. Comparing the real and imaginary parts of Eq. (22), we have the equation of motion where the density correlation is C q = ρ −q ρ q .
By solving the imaginary time equations of motion Eq. (21) and (24), one obtains the variational parameters that minimize the energy for given electronic state |ψ e .

C. Non-Gaussian Exact Diagonalization Iterations
The above two subsections outline the approach to obtain the electronic ground state with fixed variational parameters, and the ground state of variational wavefunctions with the fixed electronic state. Since the energy minimization is restricted at each step, a global ground state can be obtained only through iterations. Thus, the non-Gaussian exact diagonalization algorithm works as follows: 1. Set the initial values of the variational parameters {Γ q } and {λ q }.

Construct the effective electronic Hamiltonian in
Eq. (12) and perform exact diagonalization to obtain the (i-th iteration) electronic ground state |ψ 4. Based on the electronic many-body wavefunction |ψ (i) e , calculate the renormalized phonon energy matrix {Ω q } using Eq. (19) and the correlation functions C q , Π q and Θ q using Eq. (23). The above process is sketched in Fig. 1. Before we discuss specific parameters, we would like to briefly present an example of the NGSED iterations to give an overview of how the ground state is obtained. Fig. 2(a) shows the evolution of the energy per site (E/N ) during the iteration for λ = 2 and ω = 5. The model parameters (u, λ and ω) will be introduced and discussed later in Sec. III. For all three u values, the energy drops rapidly in the first five iterations and starts to saturate.
For this sets of model parameters, it takes ∼ 30 iterations to converge with an accuracy of 10 −6 .
To analyze the ordering tendencies of the many-body state, we evaluate the charge structure factor N (q) = ρ −q ρ q /N and spin structure factor These structure factors reflect the charge and spin ordering tendencies at certain momenta. The evolution of these observables as a function of iteration number, at the nesting momentum q = (π, π), is shown in Figs. 2(b) and (c). As the variational parameters converge, observables adjust to reflect the ordering tendencies determined by the model parameters. We will discuss the detailed parameter dependence and momentum dependence of observables in Sec. III.

III. EQUILIBRIUM PROPERTIES OF THE HUBBARD-HOLSTEIN MODEL
In this section, we apply the NGSED approach to a specific strongly correlated e-ph model and study the equilibrium properties. A typical model describing correlated electrons and phonons is the Hubbard-Holstein (HH) model 43,44 , whose Hamiltonian is The HH model is a particular example of the generic e-ph system in Eq. (1). Here, we only consider the nearest-neighbor electron hopping parametrized by the integral t, and on-site Hubbard interaction U . Both the electron-phonon coupling g and the phonon energy ω are restricted to be momentum-independent in the HH model. In this case, one can define the dimensionless ee and e-ph coupling strengths u = U/t and λ = g 2 /tω, respectively. Note, this λ is distinct from the variational parameters λ q . We adopt this notation as it is standard in the HH and non-Gaussian literature. The equilibrium phases of the Hubbard-Holstein model have been studied using different methods. Early studies have examined the equilibrium properties of the 1D HH model, using ED with optimized phonon basis 26,45 , local Lang-Firsov transformation 29 , QMC 32,46-48 , DMRG 49-51 , cluster perturbation theory 27 , and density-matrix embedding method 52 . The common results indicated CDW/AFM competition on either side of the anti-adiabatic limit u = 2λ and an intermediate region between the ordered phases. This intermediate region was originally claimed to be superconducting 46 , but more recently confirmed to be a Luther-Emery liquid with quasi-long-range charge and superconducting correlations 32 . The other extreme limit of infinite dimensions has been studied extensively using DMFT 53,54 . These studies indicate the absence of an intermediate phase.
In the context of correlated high-T c materials, the study of two-dimensional systems is more relevant. However, due to the limitations of numerical techniques, the study of the 2D Hubbard-Holstein model is relatively rare. Using determinant QMC (DQMC), Nowadnick et al. studied the phase diagram of the 2D Hubbard-Holstein model at high temperature (lower temperatures being restricted by the fermion-sign problem) and characterized the metallic phase between the competing ordered phases 55,56 . These studies were followed by ED studies at zero temperature. However, due to the infinite Hilbert-space dimensions, these ED studies of HH model was restricted to a one-phonon truncation at small clusters 57,58 or single-phonon-mode simplification [59][60][61] . These calculations, though also exact, are highly restricted by the coupling strength and fillings due to the model simplification. The variational local Lang-Firsov transformation was also applied to 2D t-J-Holstein models, with either Gutzwiller approximation 62 or exact treatment of the electrons. As mentioned above, this transformation is already close to the generic polaron transformation employed in this work, but the ignorance of the spatial fluctuations makes crucial differences in this context. More recently, the phases of the 2D Hubbard-Holstein model were examined using variational Monte Carlo (VMC), where a superconducting phase was identified in the weak-coupling limit 19,20 . However, the nature of the variational wavefunction biased the system toward superconductivity and the results have been challenged by unbiased QMC studies 33 .
With the NGSED approach, we push the ED calculation to a relatively large cluster -a 4 × 4 system, where vital high-symmetry momenta are included. Although the phonon part of the wavefunction is variational, we minimize bias by treating the electronic part as a full many-body wavefunction. We benchmark the method by comparing with DQMC in a parameter regime where the fermion-sign problem is absent. We use the parallel Arnoldi method 63 to determine the ground state wavefunction and the Runge-Kutta Dormand-Prince 5 method to solve the imaginary time evolution. In the following subsections, we first benchmark the NGSED method in the Holstein model with only e-ph interaction. Then we discuss the ground-state properties of the half-filled Hubbard-Holstein model at a fixed phonon frequency. We conclude by briefly examining the impact of phonon frequency and carrier doping in this system. we first benchmark our NGSED method with the pure Holstein model, i.e. for u = 0, because it is a case where DQMC can give exact ground-state solutions (with extrapolation to T = 0). Technically, DQMC is an unbiased numerical method for correlated fermionic models and is most efficient at high temperatures. The evaluation of low-temperature properties is usually bottlenecked by the fermion-sign problem, where the Boltzmann weight is not positive-definite. As an NP-hard problem, there are only a few models where the fermion-sign issue can be evadable, and the Holstein model is one example. Additional details regarding DQMC for the Holstein model are included in the Appendix A.
We compare the ground-state results for four different phonon frequencies -ω = 5t, t, 0.5t and 0.2t -obtained from NGSED and DQMC. As shown in Fig. 3, the charge structure factor is monotonically increasing. With large λs, the charge susceptibility approach N = 16, which is the theoretical maximal value one can reach on a 4 × 4 cluster.
In the thermodynamic limit, this charge susceptibility should always diverge with the presence of long-range charge order. For those parameter regimes accesible by DQMC, both the charge structure factor N (q) and the average energy E/N match well between these two methods. For small ω, DQMC becomes challenging at strong couplings due to prohibitively long autocorrelation times. Therefore, we compare the NGSED results with the mean-field theory (MFT) predictions for ω = 0, where MFT becomes exact. [see Appendix B for the derivations]. We find the small-ω results asymptotically approach the MFT adiabatic predictions. Interestingly, the ground-state energy, with both electrons and phonons considered, is almost independent of the phonon frequency.
Another limit of the Holstein model is the antiadiabatic limit, where the phonon frequency ω = ∞. In this limit, the phonon degrees of freedom can be integrated out, leading to an instantaneous attraction between electrons. Unlike the phonon-mediated electronic interaction V q in Eq. (14), the attraction in the antiadiabatic limit is V = 2λ, independent of momentum q 64,64 . Therefore, it leads to an on-site attraction in real spacel. Due to the infinite phonon frequency, the dressing effect becomes a virtual process, indicating that the dressing correction to the kinetic energy vanishes. Therefore, in the ω = ∞ limit, the problem exactly maps to the attractive Hubbard model with U = −2λ. It is where the Lang-Firsov transformation can exactly decouple the e-ph system. Since the phonon frequencies are all smaller or comparable to the electron bandwidth (W = 8t) and we explicitly evaluate the phonon dressing effects, the ground state properties are far away from the anti-adiabatic limit. Comparison with the anti-adiabatic limit provides intuition for the ordering tendencies and will be further discussed in the context of the Hubbard-Holstein model.
The benchmarks with exact solutions obtained from DQMC and extreme limits in the Holstein model demonstrate that the NGSED method can adequately evaluate the coupling to phonons, though both the non-Gaussian transformation and phonon states are restricted to a variational subspace of the entire Hilbert space. The full wavefunction ansatz Eq. (6) does not produce significant bias.

B. Phase Diagram of the Hubbard-Holstein Model
Having confirmed the accuracy of the method in the u = 0 limit, we move on to finite u and discuss the ground-state properties of the Hubbard-Holstein model. We first focus on the results with high-frequency phonons ω = 5t, which is relatively cheaper in computation due to less dressing effects. A brief overview of the iteration processes for a few ω = 5t systems have been shown in Fig. 2 of Sec. II C: the ground state converges to two different phases for u λ and u λ, as is seen from the structure factors. Here we present the detailed properties of this system and different phases.
Let us first look at these extreme phases, leaving the phases for u ∼ 2λ for later discussions. With the same set of parameter as Fig. 2 [λ = 4 and ω = 5], the momentum distribution of the ground-state spin and charge structure factors are shown in Fig. 4. For the udominant regime [here u = 10 λ for Figs. 4(a,b)], the system is dominant by the spin ordering, reflected by the large S(q) compared with N (q). More specifically, the spin correlation sharply peaks at q = (π, π) momentum. This reflects the tendency toward antiferromagnetism in the thermodynamic limit.At the same time, the system displays almost no charge fluctuations since the charge degrees of freedom are frozen at equilibrium. In the context of this article, we call this spin-dominant phase an "AFM phase" though we do not have spontaneous SU (2) symmetry breaking in a finite cluster. In pure variational methods with mean-field decoupling, this AFM phase indeed establishes a symmetry breaking and a spin order parameter 19,20 . On the contrary, in the λ-dominant regime [here u = 0 λ for Figs. 4(a,b)], the ground state exhibits significant charge correlations. Different from the AFM case, here only the (π, π) momentum exhibits strong correlations while other momenta are negligible. This is a difference between continuous and discrete symmetry breaking: the magnon fluctuations weaken the spin ordering in the AFM phase, while there is no Goldstone mode for the CDW phase. As expected, when charge ordering dominates, the ground state forms checkerboard doublons and holons, exhibiting no net spin correlation.
With the increase of e-e interaction u starting from the CDW phase for any fixed λ, the charge structure factor rapidly drops as shown in Figs. 5(a) and (b). There is a sharp transition near u ∼ 2λ (but slightly away from this value, see discussions below). Beyond this transition point, spin correlations rapidly build up, overwhelm the charge instability, and form the Mott AFM state. For various λs and us, we obtain the coarsegrained "phase diagram" of spin and charge structure factors in Figs. 5(c) and (d), indicating the regions of these two phases. Note that the difference between continuous and discrete symmetry breaking mentioned above leads to the distinct nature of the CDW and AFM phases. This is reflected by the ground-state degeneracy, or the excitation gap shown in Figs. 5(e) and (f). In the CDW phase (u 2λ), the ground state exhibits a two-fold degeneracy. That means, we have obtained the (π, π)-ordered Peierls phase in a 4 × 4 cluster. However, in the AFM phase (u 2λ), the ground state is nondegenerate. That being said, the SU (2) symmetry is not broken though the dominant spin susceptibility is §(π, π).
The two extreme phases described above are expected and understood. What we are more interested in is the behavior near the boundary u ∼ 2λ, where the two instabilities compete with each other. Interestingly, strong spin correlations start to build up already at u < 2λ, as reflected in Figs. 5(a) and (b). For example, in the λ = 4 system, S(π, π) becomes dominant at u = 7.4 instead of 8; whereas in the λ = 2 system, the AFM phase is reached for u ≥ 3.7 instead of 4. This is the case also for all λs in the phase diagram in Figs. 5(c) and (d). The fact that the boundary of the AFM phase sits on the u < 2λ side has been observed in 1D DMRG 50 and 2D QMC 55 studies, but was not reproduced in previous variational studies with long-wavelength Lang-Firsov transformation. Now, with the NGSED, we are able to interpret the origin of this phenomenon explicitly. For convenience, let us define the effective Coulomb interaction as In the anti-adiabatic limit, or Lang-Firsov picture, the sign of U (eff) q determines the local trend to form a doublon or spin-singlet. We find this local picture being approximately correct for systems far away from the phase boundary, as shown in Fig. 6(b): the U (eff) q are all negative for u = 6 while positive for u = 9. Though momentum-space fluctuation exists already in these cases, the ground state is qualitatively determined by the sign of the effective interaction. The overall sign accounts for the CDW and AFM at two extremes discussed above.
However, the ground-state solution for the polaronic dressing parameter λ q strongly varies over the first Brillouin zone [see Fig. 6(a)]. Due to the large charge suspectibility, the polaronic dressing, reflected in λ q , converges to a substantially larger value at the nesting momentum than other qs. In contrast to the uniform distribution assumed in the Lang-Firsov transformation and the U (eff) q ≡ U − 2λ consequence, such momentum fluctuations of λ q leads to the effective long-range interactions V q and, accordingly, the fluctuations of U  fluctuations of the effective interactions lead to exotic intermediate phases near the phase boundary. As indicated in Fig. 5, there are two narrow regimes between the welldetermined CDW and AFM phases. One intermediate regime (denoted as A) lies next to the CDW phase (e.g., 3.1 ≤ u ≤ 3.2 for λ = 2 and 6.9 ≤ u ≤ 7 for λ = 4), marked as green in Fig. 5. In this regime, the system still exhibits a large charge structure factor, but has lost the ground-state degeneracy, e.g., displays finite excitation gap. That being said, the system lies in a non-CDW phase with large charge structure factor in this narrow regime. From the perspective of U has changed sign in part of the Brillouin zone, though it is still negative at the nesting momentum. The situation in the non-CDW phase is very similar to the Luther-Emery liquid in 1D or quasi-1D system, which might display coexisting superconductivity and charge order 8,32 . In recent VMC studies, the entire intermediate regime was claimed to be superconducting 19,20 . However, due to the biased electronic wavefunction ansatz, the conclusion remains controversial 33 . With the full manybody wavefunction kept for the electrons, the NGSED calculation provides a more reliable characterization of the two intermediate phases. Unfortunately, we cannot examine the scaling of the charge or pair correlations and extract the correlation length in a finite cluster. As a compromise, we calculate the binding energy defined as Fig. 7 shows the evolution of E bd as a function of u for λ = 2 and λ = 4. For u < 3 (λ = 2) and u < 6.8 (λ = 4), the binding energy is sizable and negative and decreases (in magnitude) with the rise of u. In intermediate phase A we still see large binding energy. This might be an indication of coexisting Cooper pairs, although we cannot draw any rigorous conclusion about superconductivity in a finite cluster. One indication of an enhanced tendency toward superconductivity is the fact that the momentum dependence of λ q leads to a strongly dispersive effective phonon energy; see Eq. (19). As demonstrated in recent DQMC studies, phonon dispersion in the Holstein model may favor superconductivity over CDW 65 .
In contrast to intermediate phase A, intermediate phase B, which is adjacent to the AFM phase, exhibits both small charge and spin structure factors. Here, the effective interactions have been delicately balanced and become too weak to overcome the kinetic energy and localize electrons. Considering the possible tendency toward superconductivity present in the neighboring intermediate phase A, phase B is possibly superconducting. To investigate this possibility, we calculate the s-wave superconducting pair correlation function, defined as where the pairing operator is Note, the expectation value should be taken over the full wavefunction |Ψ instead of just the electronic wavefunction |ψ e . Different from the charge and spin structure factor, the pairing operator ∆ s does not commute with the non-Gaussian transformation U plrn . Thus, the expansion of the pairing correlation function, with electron wavefunction and variational parameters, is While permuting the electronic operators with the polaronic transformation U plrn , they physically represent the same operators of dressed quasiparticles. Therefore, to evaluate the BCS-type electronic pairs, one has to compute a superposition of FFLO-type quasiparticle pairs, because the polaronic dressing exchanges momentum between electrons and phonons. Fig. 8 presents P s as a function of u calculated for four different λs. The pairing correlation is close to 0.5 for u = 0, which is the expectation value for a CDW state; it is strongly suppressed in the AFM phase due to the low rate of double occupancy. We observe an enhancement of P s in the intermediate phase B. This enhancement is substantial for small λ, but is gradually suppressed with the increase of coupling strength. For λ > 3, the pairing correlation is no longer enhanced in the intermediate regime, but smoothly evolve between the CDW and AFM phases. If the intermediate phase B is indeed superconducting, the non-CDW phase A with large charge structure factor could be a crossover between CDW and superconductivity.
To summraize the evolution of superconductivity and the intermediate phases, we sketch a phase diagram in Fig. 8 Fig. 5. Though the λ−evolution through phase B is smooth, we denote the small-λ regime as a "strong superconductor"  and (b) spin structure factor as a function of dimensionless e-ph coupling λ for u = 8 and ω = t, 2t, 4t, 8t and 16t respectively. Diagram of (c) N (π, π) and (d) S(π, π) as a function of both u and λ. The dashed lines denote the anti-adiabatic critical line u = 2λ. The phonon frequency ω is set as t. and the large-λ regime as a "weak superconductor" or perhaps even a metal. In our calculations we find the superconducting phase is enlarged only slightly as λ decreases from 5 to 0.5. This is in contrast with the VMC predictions which assigned the entire u < 2λ, λ < 1 region as superconducting 19,20 . However, our results are consistent with the QMC conclusions in the same regime 33 . Considering that QMC is unbiased, this conclusion reflects the necessity of reliability treating the electronic wavefunction. A rigorous identification of the nature of the intermediate phases, including longrange scaling, might require future studies using DMRG together with a sophisticated treatment of phonons.

C. Impact of Phonon Frequencies and Doping
Having understood the phase diagram of the 2D Hubbard-Holstein model with a fixed phonon frequency, we briefly discuss the impact of various frequencies and carrier doping in this subsection.
Similar to the case of the Holstein model discussed in Sec. III A, we expect the e-ph system exhibiting steeper phase transitions with smaller phonon frequency. For a fixed λ, the smaller ω implies larger g/ω. As shown in Figs. 9(a) and (b) for the calculations with ω ranging from t to 16t, both charge and spin structure factors drop more rapidly for smaller frequencies when approaching the phase boundary. This trend is consistent with previous studies, e.g. finite-temperature DQMC results 55,66 and VMC results 19,20 . Intuitively, it can be understood as the adiabatic limit behaves similar to a mean-field theory, suppressing all quantum fluctuations which accumulates before reaching a phase transition. Here, using the language of the polaronic dressing in the non-Gaussian wavefunction, we provide the interpretation from a different perspective -the combined impact of polaronic dressing in both tunneling and interaction parameters. As is well-know, the Lang-Firsov transformation should give the same effective e-e interaction for a fixed λ in the atomic limit. However, the dressing parameter λ q , to generate the same V q , is larger for a smaller ω. That means, if one takes the tunneling terms into account, the polaronic renormalization for t α is larger. Therefore, the quantum fluctuations become effectively weaker with respect to the same interaction strength, leading to a sharper phase transition.
Varying both λ and u for ω = t, we obtain the phase diagram shown in Figs. 9(a) and (b) using NGSED. An immediate observation is the suppression of the intermediate phase, if it exists at all. This phase is invisible in the VMC studies on the strong-coupling side 19,20 , but is still present at finite temperature according to DQMC studies 55 . Although the effective interaction is more dispersive [see Appendix C], its impact on the electronic configuration becomes less critical, due to the suppression of quantum fluctuations as mentioned above. This accounts for the similarity of the phase diagram compared with previous ED calculations on a Peierls-Hubbard model [with only q = (π, π) phonon mode] 61 . It is worth mentioning that the convergence for smaller ω requires many more iterations since the lack of quantum fluctuations causes traps in local energy minima in the parameter space. The convergence speed can be improved by a few warm-up iterations, as discussed in Appendix C.
With the presence of finite doping, the competition between spin and charge order is not restricted to a single nesting momentum. Though both N (q) and S(q) spread out in momentum, the q = (π, π) component still dominates [the structure factors calculated at other qs are all smaller than 1.5, not shown here]. Fig. 10 shows calculations for 12.5% doping with λ = 2 and FIG. 10: Structure factors N (π, π) and S(π, π) calculated for 12.5% doping: for various interaction parameters u and (a) λ = 2, (b) λ = 4. The phonon frequency is ω = 5t. λ = 4. Both the charge and spin structure factors are significantly smaller than the half-filled case [see Fig. 5]. Interestingly, the ground-state charge structure factor for λ = 4 is not monotonically suppressed by the increase of u, in contrast to the situation at half-filling. For u < 4, the increase of electron correlations in fact slightly enhances the (π, π)-charge ordering. This trend may be regarded as a correlation-enhanced polaronic dressing effect 67 : the presence of electronic correlations reduces the mobility of carriers in a doped system, and therefore favors the polaronic dressing to some extent. A more rigorous confirmation of this non-monotonicity and a specific assessment of the underlying physics are beyond the scope of this work, and should be further investigated using a combination of multiple numerical methods.

IV. CONCLUSION AND OUTLOOK
We have presented NGSED, a new wavefunction-based method used to treat systems with both e-e and eph interactions, taking advantage of both variational non-Gaussian transformations and exact diagonalization. The variational part of the wavefunction avoids the combinatorial phonon degrees of freedom, while the full many-body electronic state minimizes bias and allows for the complexities associate with electronic correlations. We presented the formalism for this method using a generic e-ph system, where the e-ph coupling is g q , the e-e interaction is U q and the phonon energy is ω q are allowed to be momentum dependent. We applied the NGSED method to the Hubbard-Holstein model, where we compare with various other approaches. To assess the bias incurred by our variational ansatz we have benchmarked against numerically exact DQMC results on the Holstein model. The consistency with DQMC results justifies the correctness of NGSED, at least for the Holstein type of e-ph coupling.
With this new method, we have examined the groundstate properties of the 2D Hubbard-Holstein model. While the results at extreme parameter limits are consistent with known conclusions, we found interesting and delicate structures near the transition. We show that the boundary of the AFM phase is on the u < 2λ side, which is consistent with some known exact results in 1D and variational results in 2D, but has not been completely explained yet. With the information of the entangled e-ph wavefunction, we provided an intuitive picture of this boundary shift from the effective e-e interaction point of view. We demonstrate that the traditional local Lang-Firsov transformation overestimates the impact of phonons by neglecting their momentum fluctuation. The advantage of the NGSED method is reflected in part by the capability to capture these fluctuations and physically address the origin of the boundary shift.
In addition to the boundary shift, we also find two narrow intermediate phases between the CDW and AFM phases. One of them exhibits strong superconductiv-ity for relatively weak interactions, while the other is a non-CDW phase with strong charge correlations and significant binding energy. Both phases reside within the superconducting phase suggested by VMC studies 19,20 . However, the intermediate phases obtained in our NGSED calculations are much narrower and do not intersect with u = 0, a result that is supported by unbiased QMC calculations.
With the capability to adequately capture the phonon dressing, the NGSED method combines the merits of variational and exact approaches in many-body systems: it addresses the issue of both the large phonon Hilbert space and the lack of correlations in the pure variational approach. Thus, it provides a general prototype for a variety of problems involving e-e and e-ph interactions: by allowing the coupling strength g q and phonon energy ω q to vary in momentum space, it can be applied to more realistic e-ph systems like those with forward scattering, B 1g or acoustic phonons; by a rotation of the fermionic basis via the U plrn similar to Eq. (31), it can also be employed to calculate other instantaneous observables involving high-order correlations. More importantly, as a wavefunction-based method, the NGSED can be extended to excited states by sequentially projecting out the low-energy manifold; excited-state information is crucial for the evaluation of spectroscopies or thermal averages. Through the projection of equations of motion for the real-time dynamics 24 combined with advanced Krylov-subspace techniques, the NGSED method can be generalized to investigate the out of equilibrium physics in the pump-probe electron-phonon system.
The polaron transformation provides the lowest order decoupling between electrons and bosons. Extending to more intricate forms of non-Gaussian transformations, the NGSED method can be employed to decouple the interaction between electrons and other bosonic excitations, such as excitons, plasmons, and magnons. The non-Gaussian transformations have been used to study impurity models like the Kondo and Anderson models [22][23][24] , and some models in lattice gauge theory, like the 1D Schwinger model 68 , paving the way for application to Kondo-Hubbard and Anderson-Hubbard models, as well as the lattice gauge theory in higher dimensions. The study of these electron-boson or impurity problems would help to elucidate the collective and local properties of correlated materials.
More generally, numerical methods involving non-Gaussian wavefunctions offer opportunities to extend electronic structure theory. The traditional ab initio electronic structure theory is constructed on top of Gaussian states (Slater determinants), evolving into post-Hartree-Fock methods (configuration interaction, coupled cluster etc.) and multi-reference methods (CAS-SCF, CAS-CI etc.). Using the non-Gaussian wavefunctions as the fundamental basis, one can embed quantum entanglement at the outset. The NGSED method, as an analog of the full configuration interaction, can be regarded as the first building block in a non-Gaussian-based post-mean-field class of methods. Relevant post-mean-field methods constructed on this set of bases include the embedding with other many-body approaches. For example, with the same formalism handling the phonon wavefunction, the non-Gaussian transformation can be embedded with DMRG or tensor networks, self-consistently transforming a fermion-boson problem into one of quasiparticles with long-range interactions. Since the non-Gaussian transformation has rotated the many-body basis from electrons to quasiparticles, it might be helpful to reduce the fermion-sign issue in DQMC. Moreover, the multireference framework can also be extended to a non-Gaussian wavefunction basis, through the construction of superpositions of non-Gaussian wavefunctions or even NGSED.
Particle Physics (Springer, 2008), pp. 357-366. We present a brief introduction and supplementary data about the DQMC technique.
As a standard technique to many-body systems, especially the Holstein model, a detailed introduction to DQMC can be found in literatures e.g., Ref. 66,69,70. We emphasize that the notorious fermion sign-problem is absent in the Holstein model because the phonon field couples in the same way to both spin-up and spin-down electrons, so that the probability measure is proportional to the fermion determinant squared and is therefore nonnegative. The absence of a sign problem allows us to access relatively low temperatures. However, at exceedingly low temperatures DQMC calculations for the Holstein model are still limited due to prohibitively long autocorrelation times 71 . To partially mitigate this issue, we employ a combination of local and global updates, as explained in Ref. 66.
In Fig. A1 we show how the energy density E/N and structure factor N (π, π) approach their asymptotic T = 0 values for representative parameters λ = 1 and ω = 0.2, 1, 5. For all the DQMC data reported in the main text, we use the values of E/N and N (π, π) at our lowest temperature, where they have ceased to change appreciably, to approximate the value at T = 0. We note that that the requisite temperatures to probe the T = 0 limit become lower as ω increases. This trend can be understood from the fact that the Holstein model maps to the negative-U Hubbard model in the limit ω → ∞, which has a vanishing T c for coexisting SC and CDW respectively. Compared to the ω = 5 results in Fig. 5, the small-frequency system exhibit a steeper transition near the phase boundary, due to the adiabatic reasoning mentioned in the main text.
It is worth to mention that the convergence is much harder for ω = t compared t larger frequencies. It typically takes 100-200 iterations even without reaching the close proximity of the phase boundary, while the ω = 5t systems converge within 30 iterations. This is because the retardation of phonons drives the system away from an effective electronic model. The electron and phonon states have to exchange information many times to adjust to the optimal configuration. Near the phase boundary, the convergence can even be trapped by some local minima within the numerical accuracy 10 −6 , as shown in Fig. C1. Due to the suppression of quantum fluctuations, the local minima barrier becomes steeper. To overcome this issue, we add "warm-up" iterations for larger ω but with the same λ [i.e. using g = √ αg and ω = αω where α is a scaling factor much larger than 1]. These iterations are relatively faster and give raw approximations for the ground-state configurations at small frequencies, avoiding possible local minima. Fig. C1 shows the results for the ground-state energy and structure factors obtained using and without using "warm-up" iterations. For systems near a phase transition, inappropriate treatment of the convergence may lead to a completely incorrect phase near the transition, though close in energy. The results in Fig. 9 was obtained by asymptotically tuning the scaling factor from 16, 8, 4, 2 to 1. In addition to the ground-state energy and structure factors, we also present the effective interaction U (eff) q in Fig. C2. The interaction is more dispersive compared to large ωs, indicating the effective interactions mediated by the phonon become longer-range in the adiabatic limit. However, as mentioned in the main text, the suppression of the quantum fluctuations occurs exponentially; therefore, these interactions become semi-classical and lead to a sharp mean-field-like transition near the phase boundary.

Appendix D: Comparison with NGS-GS Method
In the limit of the Holstein model, previous studies have shown that the non-Gaussian transformation well describes the ground-state properties [19][20][21] . To make a specific comparison, we present the calculation based on a pure variational ansatz Ψ(t) = U plrn (t)|ψ GS ph ⊗ |ψ GS e .

(D1)
Here, the Gaussian phonon wavefunction |ψ GS ph and the non-Gaussian transformation U plrn are the same as the definition in Eq. (7) and (8). In contrast to the full many-body wavefunction |ψ e , the electronic part is also replaced by a Gaussian wavefunction The ground-state energies calculated using this non-Gaussian + Gaussian ansatz are summarized in Fig. D1. For most frequencies and coupling strengths, this ansatz is consistent with the results of NGSED, indicating that the electronic state indeed forms CDW orders in these cases. Only on the small ω and large λ limit, the NGS+GS ansatz starts to deviate (slightly) from the NGSED. This can be attributed to the fact that the dressing factor λ q , in this case, becomes huge and causes stronger fluctuations.  Fig.3(b)] and the g NGS + Gaussian wavefunction ansatz (blue squares), for phonon energies ω = 0.2t, 0.5t, t and 5t.