Density-matrix embedding theory study of the one-dimensional Hubbard-Holstein model

We present a density-matrix embedding theory (DMET) study of the one-dimensional Hubbard-Holstein model, which is paradigmatic for the interplay of electron-electron and electron-phonon interactions. Analyzing the single-particle excitation gap, we find a direct Peierls insulator to Mott insulator phase transition in the adiabatic regime of slow phonons in contrast to a rather large intervening metallic phase in the anti-adiabatic regime of fast phonons. We benchmark the DMET results for both on-site energies and excitation gaps against density-matrix renormalization group (DMRG) results and find excellent agreement of the resulting phase boundaries. We also compare the fully quantum treatment of phonons against the standard Born-Oppenheimer (BO) approximation. The BO approximation gives qualitatively similar results to DMET in the adiabatic regime, but fails entirely in the anti-adiabatic regime, where BO predicts a sharp direct transition from Mott to Peierls insulator, whereas DMET correctly shows a large intervening metallic phase. This highlights the importance of quantum fluctuations in the phononic degrees of freedom for metallicity in the one-dimensional Hubbard-Holstein model.


INTRODUCTION
The interplay of competing interactions is a central theme of quantum many-body physics. In quantum materials, the electron−electron (el−el) and electron−phonon (el−ph) interactions naturally compete against each other. This is most easily understood by noting that el−el Coulomb repulsion is generically repulsive, whereas el−ph interactions can lead to effectively attractive el−el interactions, as highlighted by the Cooper pairing mechanism in conventional superconductors. 1 In strongly correlated low-dimensional materials, the competition between el−el and el−ph interactions has led to longstanding debates, such as the origin of high-temperature superconductivity and the anomalous normal states observed in entire classes of materials. 2,3 At the same time, competing interactions lead to competing ground states and phase transitions, as exemplified by the complex phase diagrams in correlated oxides. 4 This competition and sensitivity to parameter changes poses a major roadblock on the way toward reliable numerical solutions for the quantum many-body electron−phonon problem. The simplest generic el−ph model is the Hubbard−Holstein Hamiltonian. 5,6 Advanced numerical methods for its solution have been developed over the past decades that are accurate in certain regimes but cannot be easily applied in other cases. Among these are Quantum Monte Carlo (QMC) 7−13 and variational wave function 14−27 schemes and the density-matrix renormalization group (DMRG), 28 as well as dynamical mean-field theory (DMFT). 29−31 A new, promising method that has recently been added to this arsenal is density-matrix embedding theory (DMET), 32,33 which in some sense bridges DMRG and DMFT-related methods. DMET has the advantage of being numerically less demanding than DMRG and at the same time being a good descriptor for one-dimensional systems. DMET has been benchmarked against other methods for the 2D Hubbard model, 34,35 and recently, a systematic extension of DMET toward DMFT was proposed. 36 However, for the Hubbard−Holstein model, only one DMET study has been published to the best of our knowledge so far, which compared ground state energies for the 1D Hubbard−Holstein model against DMRG results. 37 In this work, we perform extensive comparisons of DMET against both DMRG and Born−Oppenheimer (BO) results for the 1D Hubbard−Holstein model. In particular, we extend the previous comparisons to excitation gaps and the difference of the electron density between neighboring sites, which indicates the existence of a charge-density wave (CDW). This allows us to construct DMET phase diagrams that are compared directly against DMRG.
The paper is organized as follows: In Section 2, we explain the extension of the DMET method for electron−phonon systems, following the derivation presented in ref 37. In Section 3, we discuss the Hubbard−Holstein model, its phases, and possible observables to determine these phases. In Section 4, we discuss our DMET results and benchmark them against a DMRG calculation. Moreover, we then benchmark a semiclassical Born−Oppenheimer calculation against DMET and show when the Born−Oppenheimer approximation fails.

DMET FOR COUPLED ELECTRON-BOSON SYSTEMS
When trying to solve the Schrodinger equation for a given general Hamiltonian, the well-known problem of the exponential wall of quantum mechanics occurs, making the costs of the calculation grow exponentially with system size. Even though there are wave function methods that scale this problem down to polynomial growth, it is still a fact that wave function methods grow fast with the size of the regarded system, making it hard to describe large systems. One possible way to circumvent this problem is the embedding idea that is used in different methods including DMET: Instead of solving the Schrodinger equation for the whole system, a small subsystem is chosen, which is small enough to be solved efficiently. The idea of DMET is then to include the interactions of the rest of the system with this subsystem, which will from now on be called "impurity", in the embedding step. This is equivalent to a complete active space calculation in quantum chemistry, assuming that the impurity part of the system is in the active space. This way, the system is in two unentangled parts: the socalled embedded system, consisting of the impurity and those parts of the rest of the system interacting with the impurity, and the environment that consists of those parts of the system that do not interact with the impurity directly. Then, by solving the embedded system, the physics of the impurity, including interactions with the rest of the system (and with that, also finite size effects and the influence of the boundaries) can be computed effectively with an accurate wave function method, since the embedded system is typically much smaller than the originally considered system. Mathematically, the projection from an original basis to the impurity plus active space basis can be formulated with the help of the Schmidt decomposition:  corresponds of the dimension of the Fock space defined on the impurity part of the system and b  corresponds to the dimension of the Fock space defined on the rest of the system. Every wave function can be decomposed into two parts |̃⟩ a and |̃⟩ b , where the former is defined on the impurity and the latter, on the rest of the system. Both parts are coupled to each other by the transition matrix Ψ ab . Performing a singular value decomposition of this matrix leads then to a new basis consisting of the many-body states | ⟩ i and | ⟩ i . The number of many-body states describing the whole system is then significantly less than before, as the sum over i in eq 2 is only going over the dimension of the impurity Fock space.
Knowing the many-body states | ⟩ i and | ⟩ i , a projection can be defined, that projects the Hamiltonian onto a new basiŝ Ĥe mb is then a many-body Hamiltonian of dimension × 2 i  , defined on a subspace of the original Fock space. For small N imp , it can be diagonalized efficiently with a chosen wave function method.
Unfortunately though, in order to find the active space states | ⟩ i , the whole transition matrix Ψ ij , that is, the whole wave function |Ψ⟩, needs to be known. This is why, instead of using the exact projection, we have to approximate it in order to find the embedding Hamiltonian. Note that we only approximate the projection of the Hamiltonian into a new basis and not the Hamiltonian itself.
In this paper, we discuss the DMET method for the coupled electron−phonon systems whose Hamiltonian is of the form: where Tê l and Ûe l describe the electronic kinetic energy and the electron−electron interactions. Tp h describes the phononic kinetic energy, and Ûe l−ph is the electron−phonon interaction. To obtain the projection, we treat both the electronic and the phononic part of the Hamiltonian on a mean-field level such that there is no coupling between electronic and phononic degrees of freedom.
As the ground state wave function of eq 8 is a product state of the electronic and the phononic wave function, also the projection of eq 5 can be factorized into an electronic and a phononic projection. We thus have to find two separate projections, one for the electrons and one for the phonons (visualized in Figure 1). As the procedure for obtaining the electronic projection has already been widely discussed in refs 32 and 33, we will only very briefly describe it here, while we will

Journal of Chemical Theory and Computation
Article present the derivation of the phononic projection in detail. Subsequently, we derive the coupled electron−phonon embedding Hamiltonian eq 6 from the original Hamiltonian in eq 7 with the two projections.
2.1. Fermionic Projection. As an initial guess for the electronic mean-field Hamiltonian, we simply choose electronic kinetic energy: =Ĥ T el mf el (9) which will later be improved in a self-consistent manner by adding a nonlocal potential V̂e l . The ground state of Tê l can be written in terms of a Slater determinant, which can be further decomposed to 33 Here, is the dimension of the electronic impurity Fock space and is the dimension electronic environment Fock space. i  is significantly smaller than j  as we choose the impurity to be much smaller than the total system. Neglecting the environment part of eq 10 as this is the part that does not couple to the impurity region of the system, we obtain 2.2. Phononic Projection. In order to obtain the phononic projection, we choose as initial guess a Hamiltonian describing a shifted harmonic oscillator, which has the same form as the phononic part of the electron−phonon Hamiltonian of eq 7: Equivalently to the electronic case, this Hamiltonian will also be improved self-consistently. The ground state wave function of this Hamiltonian is the product state of coherent states on each lattice site i: where z i = ⟨aî † + aî⟩ is the shift of the phonon mode from the initial position on the lattice site i. Note that the state |z i ⟩ is a superposition of all possible phononic Fock numbers states on site i. In the original Hamiltonian Ĥe l−ph , defined in eq 7, due to the coupling term between electrons and phonons, the total number of phonons is not conserved (which makes sense as they are only quasi-particles describing the lattice vibrations of the solid). As the coherent state defined here in eq 13 also does not obey phononic particle number conservation, it is well suited to describe our problem. Similar to the electronic case, 32,33 we determine the phononic projection by splitting the phononic wave function up into three parts: Here, the |A i ph ⟩ is again just defined on the impurity region of the lattice, |B l ph ⟩ is composited by those coherent states that have an overlap over the impurity region, and |B̃j ph ⟩ is composited by the coherent states that do not have an overlap over the impurity region. The , similar to the electronic case, determine the dimension of the phononic Fock spaces on the impurity i  and on the rest of the system j  . Due to their bosonic nature, the number of phononic basis functions is in principle indefinitely N ph = ∞. In the numerical calculation, we restrict ourselves to only a few phononic basis functions, as otherwise the computation would not be feasible. Neglecting again the part of the wave function that does not have an influence on the impurity, we define the phononic projection as

Projection of the Full Hamiltonian.
Knowing the electronic as well as the phononic projection, we are now able to find the embedding Hamiltonian of the coupled system. The coupled electron−phonon embedding Hamiltonian eq 6 is obtained by projecting the purely electronic part of the original Hamiltonian (Tê l and Ûe l ) with the electronic projection Pê l and the purely phononic part Tp h with the phononic projection Pp h . The coupling term Ûe l−ph emb needs to be projected with both the electronic and the phononic projection.  Figure 1. Visualization of the decomposition of the electron−phonon system via the projection P: Starting with a 1D lattice in real space that on each site has both electronic (blue) and phononic (red) degrees of freedom, we choose one part of the system that is from then on called impurity, whereas the rest is the environment. The electronic and the phononic sites are treated equally in this scheme by defining the full projection as a product of one projection for the electrons and one for the phonons. This projection then leaves the basis on the (electronic and phononic) impurities the same (a real space lattice), whereas it projects the (electronic and phononic) environment degrees of freedom into a new basis set whose physical meaning is abstract. Within this new basis set, the environmental degrees of freedom can be divided into those interacting with the impurity and those not interacting with the impurity, called environment. The system containing the impurity and the basis sites interacting with the impurity is called embedded system. In our calculation, we discard the environment system and only calculate the embedded system.

Journal of Chemical Theory and Computation
Article ̂=̂̂ † T P T P ph emb ph ph ph (20) =̂̂̂̂− † † − † U P P U P P el ph ph el el ph emb ph el (21) 2.4. Self-Consistently Improving the Projections. The initial guess for the projections is not necessarily very good, as in addition to assuming a noninteracting active space for both the electrons and the phonons, it also assumes a product state between electronic and phononic degrees of freedom.
We self-consistently optimize the electronic and the phononic projection by adding nonlocal potentials to the mean-field Hamiltonian eq 8: For the electronic case, as was explained in detail in previous works, 32,33 we optimize the nonlocal potential such that the electronic reduced one-particle density matrix of the interacting and the noninteracting system is minimized: For the phonons on the other hand, the mean-field Hamiltonian contains two terms, one for the kinetic energy (Tp h ) and one for the shift of each harmonic oscillator from its initial position (Ĉp h ). This is why we also have to add two nonlocal potentials (V̂p h and Ŵp h ) to the mean-field Hamiltonian. While (26) is optimized such that it minimizes the phononic reduced oneparticle density matrices of the interacting and the noninteracting system, (27) minimizes the difference between the shift of the interacting and the noninteracting system.
While V̂p h depends on the phononic reduced one-particle density matrix ⟨aî † aĵ⟩, Ŵp h depends on the shift of the phonons from their rest position ⟨aî † + aî⟩. The potentials are again found by minimizing the difference between the properties of the interacting and the noninteracting system: emb emb emb emb emb (28) where |Z⟩ is the ground state wave function of the Hamiltonian defined in eq 23 and |Ψ emb ⟩ is the ground state wave function of the full embedding Hamiltonian defined in eq 17. The whole DMET procedure is visualized in Figure 2.
The embedding problem is then solved using MPS-DMRG. 38−40 We obtain the optimal site ordering 41 from an initial approximate calculation. 42 This site ordering is then used in a second higher-precision calculation. In both cases, we construct the Hamiltonian as described in ref 43. Electronic and phononic sites were kept separate, and the DMRG3S algorithm 40 was used to achieve linear scaling with the relatively large local dimension.

HUBBARD−HOLSTEIN MODEL
The Hubbard−Holstein model is described by the Hamiltonian Here, cî σ ( †) is the electronic particle creation (annihilation) operator on site i, where σ (= ↑, ↓) is the spin degree of freedom, is the spin-up (spin-down) particle number operator, and aî ( †) is the phononic particle creation (annihilation) operator. The kinetic energy of the electrons is approximated as a simple next-neighbor hopping term T̂e l , and the electron−electron interaction is assumed to be a purely local Hubbard U term Ûe l . The phonons are approximated by harmonic oscillators T̂p h which are bilinearly coupled to the density of the electrons Ûe l−ph . One phonon mode is considered per electronic site.
3.1. Physics. In the one-dimensional Hubbard−Holstein model, three competing forces can be found: first, the electron hopping strength t, that leads to mobilization of the electrons and will put the system in a metallic phase; second, the electron−electron interaction U that, if dominant, leads to an immobilized spin wave for the electronic degrees of freedom, that is, a Mott phase; third, the electron−phonon coupling g that, if dominant, leads to a Peierls phase, which is the position of the electrons on the lattice if distorted from the initial position, forming a charge density wave.

Journal of Chemical Theory and Computation
Article Due to strong quantum fluctuations of the phonons, the Peierls phase can be destroyed and can lead to a metallic phase, if the electron−electron interactions are not too strong to prevent this. This is why we expect a distinct metallic phase when considering high phonon frequencies ω 0 in comparison to the itineracy of the electrons t. In contrast, if the phonon frequency is small compared to the electronic hopping, the metallic phase should, if existent, be smaller than in the anti-adiabatic limit.
3.2. Observables Monitoring the Phase Transition. In order to describe the phase transition of the one-dimensional Hubbard−Holstein model, we need to define observables that unambiguously show which phase the system is in.
We choose here three different observables, namely, the double occupancy ⟨n i↑ n i↓ ⟩, the electronic density difference between neighboring sites ⟨n i ⟩ − ⟨n i+1 ⟩, and the energy gap Δc defined in eq 30.
The double occupancy and the electronic density difference between neighboring sites are local properties and can simply be calculated on one arbitrary (impurity) site. Unfortunately though, the double occupancy only gives a rough estimate of the phase, and the electronic density difference between neighboring sites only indicates the transition to the Peierls phase; in the Mott phase, the electronic density stays homogeneous.
The energy gap Δc indicates unambiguously whether the system is in the metallic phase (where the gap is zero) or in an insulating state, which can be either Mott or Peierls (where the gap is finite). Unfortunately though, it cannot be simply defined locally but is a global property of the whole system for different particle numbers where E 0 N/2 is the ground state energy of the Hamiltonian for half filling, E 0 N/2 − 1 is the ground state energy of the system for half filling minus one, and E 0 N/2 + 1 is the energy of the system for half filling plus one.
As our DMET calculation has only been implemented for closed shell systems, we have to approximate the calculation of the gap. Instead of doing three calculations with different particle numbers, we consider our "sophisticated mean-field" Hamiltonian which is optimized to have similar one-body properties as the interacting Hamiltonian. We calculate the (one-body) spectrum of this Hamiltonian by diagonalizing it and then approximate the gap by defining 3.3. Parameters. The phase transition depends on the itineracy of the electrons (∝ t), the electron−electron repulsion (∝ U), the electron−phonon interaction (∝ g), and the relative velocity of the phononic degrees of freedom with respect to the electrons. This is why we introduce the adiabaticity ratio 34) accounting for the relation between the kinetic hopping energy of the electrons t and the frequency of the phonons ω 0 . We also decide to discuss our results in terms of dimensionless coupling constants: In order to obtain trustworthy results, we have to make sure that our results do not depend on the implementation details, or if that is not possible, we have to quantify the impact of our implementation on the results. We therefore conducted convergence tests for DMET with respect to the total system size N, the number of considered impurity lattice sites N imp , and the number of phononic basis functions per site N ph . The results of these tests are discussed in the Appendix. Summarizing these results, we could remove the dependence from the total system size N by linear extrapolation (Finite Size Scaling). An accurate extrapolation of the results with N imp and N ph was not possible for some parameter ranges (see Number of Phononic Basis Functions and Scaling with Impurity Size, respectively). Maintaining the balance between accuracy and numerical costs, we consider an impurity system size of N imp = 6 and a total number of bosonic basis functions of N ph = 8. These parameters, together with a bond dimension of 2000 for the DMRG impurity solver, are used throughout this paper, if not stated otherwise.
4.1. DMET Results. 4.1.1. Anti-adiabatic Limit. In Figure 3, we plot the energy gap Δc/t, the electronic density difference between neighboring sites ⟨n i ⟩ − ⟨n i+1 ⟩, and the double occupancy ⟨n i↑ n i↓ ⟩ (as defined in Section 3.2) in the antiadiabatic limit (α = 5.0) for an electron−electron repulsion of u = 1.0 and for different electron−phonon coupling strengths λ.
For all three observables, we observe a Mott phase for 0 ≤ λ ≤ 0.7. With growing λ, we indeed observe a distinct metallic phase (0.7 ≤ λ ≤ 1.1) which is followed by a Peierls phase for 1.1 ≤ λ. We summarize additional data in the phase diagram for different values of u and λ in the anti-adiabatic limit in Figure 4. We observe that the system always undergoes a distinct metallic phase when transiting from Mott to Peierls phase. This behavior agrees with the assumption that strong quantum fluctuations of the phonons lead to a destruction of the Peierls phase and also prevent the emergence of the Mott phase for weak electron− electron interaction, overall resulting in a metallic behavior. The range of the metallic phase stays approximately constant.
4.1.2. Adiabatic Limit. The occurrence of a pronounced metallic phase in the anti-adiabatic limit was to be expected; it is

Journal of Chemical Theory and Computation
Article however not clear whether this phase also occurs for all electron−electron interaction strength u in the adiabatic limit, where the phonon frequency is small in comparison to the electronic hopping and, thus, the quantum fluctuations of the phonons are suspected to be smaller.
In Figure 5, we again show the energy gap Δc/t, the electronic density difference between neighboring sites ⟨n i ⟩ − ⟨n i+1 ⟩, and the double occupancy ⟨n i↑ n i↓ ⟩ (as defined in Section 3.2) in the adiabatic limit (α = 0.5) for different electron−electron repulsions (u = 0.0; 0.2; 0.4), and different electron−phonon coupling strengths λ.
When the electron−electron interaction is absent, we do not observe a Mott phase, but there is a direct transition from the metallic to the Peierls phase at λ = 0.2. This result is expected as the Mott phase is driven by the electron−electron interaction and therefore cannot occur in this limit.
For a small electron−electron interaction, u = 0.2, the Mott phase exists for very small electron−phonon interactions 0 ≤ λ ≤ 0.1. The gap indicating the Mott phase though is very small in comparison to the gap that indicates the pronounced Peierls phase for 0.3 ≤ λ. Between the Mott and the Peierls, we observe a small metallic phase for electron phonon coupling values of 0.1 ≤ λ ≤ 0.3.
When considering bigger electron−electron interactions u = 0.4, the size of the gap indicating the Mott gap grows considerably, as does the range of the Mott phase: for 0 ≤ λ ≤ 0.25, we observe a Mott phase, followed again by a narrow metallic phase for 0.25 ≤ λ ≤ 0.45. Afterward, we observe a Peierls phase, whose gap is less pronounced than for lower u but still clearly visible.
Our results for the adiabatic limit of the Hubbard−Holstein model are summarized in the phase diagram shown in Figure 6. We observe that the Mott phase, while not existent at all for u = 0.0, grows more and more pronounced for growing electron− electron interaction values u. The range of the metallic phase shrinks with increasing electron−electron interaction until it vanishes completely for u = 0.4. The position of the metallic phase shifts from small values for electron−phonon interaction λ to intermediate values. Also, the Peierls phase is getting less pronounced for growing u.

Comparison DMET and DMRG Calculations.
We benchmark the accuracy of DMET against results obtained using the DMRG method. The results are obtained with the SyTen library, which for this purpose was expanded to be able to treat coupled Fermion-boson systems. Similar to the DMET calculation, the extrapolation of DMRG results with respect to the system size was performed for all data (see Finite Size Scaling). Also for DMRG, an accurate convergence with the number of phonon basis function per site N ph was not possible for all parameter ranges, as discussed in greater detail in Number of Phononic Basis Functions. Again, maintaining the balance between accuracy and computational costs, we pick N ph = 10 for all DMRG calculations throughout the paper as well as a bond dimension of 4000.
We compare the DMRG and the DMET results for both the anti-adiabatic limit (α = 5.0) and the adiabatic (α = 0.5) limit.
In Figure 7, we compare the DMRG and the DMET results for the anti-adiabatic limit (α = 5.0) and an electron−electron repulsion of u = 1.0. Up to an electron−phonon coupling value of λ = 1.1, we observe a quantitative agreement between DMET and DMRG for all three considered quantities including the energy gap, although different approximations were invoked to calculate this property. For larger values of λ inside the Peierls phase, DMET overestimates both the size of the energy gap (computed from the mean-field value within DMET as detailed above) and the staggered charge-density compared to the respective values within DMRG. Apparently quantum fluctuations that reduce both of the reported quantities are important close to the phase transition between the metallic and the Peierls phase, and these quantum fluctuations are more pronounced in the DMRG results compared to DMET here.
However, the point of the quantum phase transition between the metallic and the Peierls phase is predicted equivalently for the DMRG and the DMET calculation within the chosen parameter grid.
In Figure 8, we compare the DMRG and the DMET results for the adiabatic limit (α = 0.5) and an electron−electron repulsion of u = 0.2. While the position of the phase transition between the Mott and the metallic phase agrees within the

Journal of Chemical Theory and Computation
Article chosen parameter grid, the position of the phase transition between the metallic and the Peierls phase slightly disagrees between DMET and DMRG. DMRG predicts the occurrence of the phase transition for higher values of λ than DMET. Similar to the above discussion of discrepancy between energy gaps and staggered charge densities in this parameter regime, it is possible that DMRG is more sensitive to quantum fluctuations here compared to DMET and, thus, predicts a more robust metallic phase compared to DMET. Another reason for the discrepancy can be that both methods do not predict the point of the phase transition accurately in this limit due to computational limitations. As discussed in Scaling with Impurity Size, an accurate extrapolation of DMET results with the impurity size N imp is not possible in the Peierls phase. Here, especially, the size of the charge gap might be overestimated and the point of the phase transition can be predicted for too small values of λ, whereas in DMRG the point of the phase transition might be predicted for too large values of λ due to the lack of convergence with N ph in this region of the phase diagram (see Number of Phononic Basis Functions).

Journal of Chemical Theory and Computation
In contrast to eq 29, here the bosonic degrees of freedom are not considered in terms of phonons but in terms of the distortion from the initial position xî of the ions. pî is the momentum of the ions.
In the BO approximation, we assume the ions to be classical particles as, due to their higher mass, they are moving much slower than the electrons. Thus, we can neglect their kinetic energy which yields Here, the remaining ionic term, ∑ω x i i 2 2 0 2 , purely depends on the distortion of the ions and can be treated as an external parameter. We treat the BO Hamiltonian with purely electronic DMET and optimize the distortion of the ions xî to minimize the total energy. The optimization of the x i values was performed on the impurity lattice sites of the system. The calculation was performed for growing impurity sizes (N imp = 2, 4, 6) and converged at N imp = 4.
In Figure 9, we compare the double occupancy ⟨n i↑ n i↓ ⟩ and the distortion of the electronic density ⟨n i ⟩ − ⟨n i+1 ⟩ for the BO system and the full quantum mechanical system in the antiadiabatic limit (α = 5.0) and for an electron−electron repulsion of u = 1.0.
We observe that, for both observables, the Born− Oppenheimer description of the phase transition is not accurate. While in the full quantum mechanical model, the transition between metallic and Mott phase occurs for a value of λ = 1.1, in the Born−Oppenheimer model, this transition already occurs for λ = 1.0. Additionally, the actual phase transition is of second order, while the Born−Oppenheimer treatment predicts a phase transition of first order.
In Figure 10, we again compare the full quantum mechanical treatment with the BO approximation, this time for the adiabatic limit and an electron−electron repulsion of u = 0.2. While still not accurate (the phase transition is predicted too early, at λ = 0.25 (BO) instead of λ = 0.3 (full)), at least the qualitative nature of the phase transition as being second order is grasped.
This result confirms our expectation that, in order to treat the quantum phase transitions of the Hubbard−Holstein model, the quantum mechanical nature both of the electrons and of the phonons needs to be taken into account. Especially, when the phononic frequency is higher than the electronic kinetic hopping, the BO approximation, which assumes the phonons to be moving much slower than the electrons, fails.

CONCLUSIONS AND OUTLOOK
In conclusion, we have benchmarked the density-matrix embedding theory against density-matrix renormalization group results for the one-dimensional Hubbard−Holstein model.

Article
We have demonstrated good agreement not only for ground state energies but also notably of excitation gaps and phase diagrams between DMET and DMRG.
An important prospect of DMET for the electron-boson problem lies in its possible extensions to electron-photon systems. Notably, recent efforts toward cavity quantumelectrodynamical engineering of materials properties 44−51 have been made. We envision that these developments will open up a whole new field in which efficient methods able to deal with correlated electron-boson lattice systems from weak to strong coupling are urgently needed. Our benchmark study helps pave the way to these new endeavors.

■ APPENDIX
In this Appendix, we will discuss how our results depend on implementational details and how we remove this dependency, if possible, or which parameters we chose for reliable calculations.
Dependence of following parameters has to be investigated for both DMET and DMRG calculations: finite size effects with respect to the size of the full system N and number of phononic basis functions per site N ph . This analysis is done in Finite Size Scaling and Number of Phononic Basis Functions, respectively.
Additionally, we analyze the dependence of DMET results on the chosen impurity size N imp in Scaling with Impurity Size.
Although not being an observable of physical interest, the energy per site is an important property to show how well two methods agree with each other. This is why, in Energies, we show the energy per site for the DMRG method and both implementations of the DMET method (with and without the Born−Oppenheimer approximation).
Finite Size Scaling. The Hubbard−Holstein model is defined in infinite space and translationally invariant. Numerically though, we are only able to consider finite systems and therefore have to consider finite size effects and the influence of the boundaries on the observables. This is why, both for the DMRG and for the DMET calculation, we do a finite size scaling.    . Dependence of the energy gap on the number of phononic basis functions in the DMET calculation for both the adiabatic limit (α = 0.5, u = 0.2) and the anti-adiabatic limit (α = 5.0, u = 1.0). In the antiadiabatic limit, the size of the gap converges with 10 phonons in all regions of the phase diagram. In the adiabatic limit, the size of the gap converges fast in the Mott and metallic phase but does not converge in the Peierls phase. The point of the phase transition from the metallic to the Peierls phase, however, converges.

Article
As the numerical costs grow quadratically with increasing total system size, we can regard very big systems and therefore did the finite size scaling with system sizes of N = 204, N = 408, and N = 816 sites, as shown in Figure 11.
As the other observables, namely, the density difference of the electrons between neighboring sites ⟨n i ⟩ − ⟨n i+1 ⟩ and the double occupancy ⟨n i↑ n i↓ ⟩ are local properties, the finite size effects and the influence of the boundaries do not influence the results anymore for system sizes bigger than N = 204.
In the DMRG calculation, opposed to the DMET calculation, we only have two sources of possible errors due to finite size effects: the system size itself and the maximal number of considered basis functions in the phononic Fock space, N phon . The numerical costs of these calculations also grow polynomially with growing system sizes. This is why, for our extrapolation, we chose to consider system sizes of N = 24, N = 48, and N = 96, as can be seen in Figure 12.
Number of Phononic Basis Functions. In Figures 13 and 14, we show the dependence of the energy gap on the number of phononic basis functions in both DMET and DMRG calculations for both the adiabatic limit (α = 0.5, u = 0.2) and the anti-adiabatic limit (α = 5.0, u = 1.0).
In both cases, we observe a convergence of the gap size in the anti-adiabatic limit for all regions of the phase diagram (DMET with 8 phonons and DMRG with 10 phonons). In the adiabatic limit, a convergence of the gap size in the Mott and the metallic phase can be observed, but neither DMET nor DMRG seem to converge in the Peierls phase. This makes predictions regarding the actual size of the energy gap in the Peierls phase and also the comparison between DMET and DMRG results difficult in this region. The position of the phase transition from the metallic to the Peierls phase converges for DMET but not for DMRG.
Scaling with Impurity Size. The results of a DMET calculation always depend on the chosen size of the impurity N imp . The finite size effects due to the size of the impurity, however, cannot be taken into account that easily, as their scaling is not linear and therefore cannot be rescaled easily. Figure 14. Dependence of the energy gap on the number of phononic basis functions in the DMRG calculation for both the adiabatic limit (α = 0.5, u = 0.2) and the anti-adiabatic limit (α = 5.0, u = 1.0). In the antiadiabatic limit, the size of the gap converges with 8 phonons in all regions of the phase diagram. In the adiabatic limit, the size of the gap converges fast in the Mott and metallic phase but does not converge in the Peierls phase. The point of the phase transition from the metallic to the Peierls phase might be predicted for too large values of λ. Figure 15. Scaling of the energy gap with growing impurity size (DMET calculation) for the adiabatic limit (α = 0.5, u = 0.2) as well as for the anti-adiabatic limit (α = 5.0, u = 1.0). We see that the discrepancy between the results gets smaller for growing impurity sizes. Figure 16. Scaling of the energy gap with the inverse impurity size (DMET calculation) for the adiabatic limit and some values of λ in the Mott and the Peierls phase. We observe a nonlinear scaling, which makes an extrapolation with impurity size not possible. Figure 17. Comparison of the energy per site E site , calculated with the DMRG and with the DMET method. In the upper graph, we show the anti-adiabatic limit (α = 5.0, u = 1.0), and in the lower graph, the adiabatic limit (α = 0.5, u = 0.2). For both limits, the results agree quantitatively.

Article
We show the scaling for the energy gap and different impurity sizes in Figure 15 for the adiabatic and the anti-adiabatic limit. With growing impurity sites, the estimation of the energy gap in the Peierls phase gets smaller for both cases and the discrepancy between the results for growing impurity sites get smaller.
To show the scaling explicitly, we plot the energy gap as a function of the inverse impurity size 1/N imp for the adiabatic limit and multiple values of λ ( Figure 16). In contrast to the scaling with the full systems size, we cannot observe a linear behavior here, making it hard to extrapolate the results or give a quantitative error estimate. We observe that overall the size of the gap in the Peierls phase decreases with increasing N imp , suggesting that in the limit of large N imp a gap closing is possible at higher values of λ. Hence, it is possible that, for N imp = 6, which has been used for all DMET calculations throughout the paper, the point of the phase transition from metallic to Peierls phase is predicted for too small values of λ.
Energies. In order to benchmark the results of the DMET calculation, we compare the results for the calculated energy per site E site with those from the DMRG calculation. In Figure 17, we show the energy per site for the anti-adiabatic (α = 5.0, u = 1.0) as well as for the adiabatic limit (α = 0.5, u = 0.2) for DMRG and DMET calculations. For both cases, the results agree on a quantitative level.
Additionally, we also compare the energies per site between the full Hubbard−Holstein model and the Hubbard−Holstein model with BO approximation in Figure 18. For the antiadiabatic limit, the energy per site shows approximately the same behavior while not agreeing quantitatively. In the adiabatic limit, a qualitative agreement can be observed.

Notes
The authors declare no competing financial interest. Figure 18. Comparison of the energy per site E site for the full electron− phonon system (DMET calculation), with the energy per site from the same system under the BO approximation. In the upper graph, we show the anti-adiabatic (α = 5.0, u = 1.0), and in the lower graph, the adiabatic limit (α = 0.5, u = 0.2). While the behavior only approximately coincides for the anti-adiabatic case, a qualitative agreement between the two methods for the adiabatic limit can be observed.

Journal of Chemical Theory and Computation
Article