Ground-state phase diagram of the three-band Hubbard model from density matrix embedding theory

We determine the ground-state phase diagram of the three-band Hubbard model across a range of model parameters using density matrix embedding theory. We study the atomic-scale nature of the antiferromagnetic (AFM) and superconducting (SC) orders, explicitly including the oxygen degrees of freedom. All parametrizations of the model display AFM and SC phases, but the decay of AFM order with doping is too slow compared to the experimental phase diagram, and further, coexistence of AFM and SC orders occurs in all parameter sets. The local magnetic moment localizes entirely at the copper sites. The magnetic phase diagram is particularly sensitive to (cid:2) pd and t pp , and existing estimates of the charge transfer gap (cid:2) pd appear too large in so-called minimal model parametrizations. The electron-doped side of the phase diagram is qualitatively distinct from the hole-doped side and we ﬁnd an unusual two-peak structure in the SC in the full model parametrization. Examining the SC order at the atomic scale, within the larger scale d x 2 − y 2 -wave SC pairing order between Cu-Cu and O-O, we also observe a local p x ( y ) [or d xz ( yz ) ] symmetry modulation of the pair density on the Cu-O bonds. Our work highlights some of the features that arise in a three-band versus one-band picture, the role of the oxygen degrees of freedom in new kinds of atomic-scale SC orders, and the necessity of re-evaluating current parametrizations of the three-band

reasons to go beyond the one-band picture to study the original three-band model directly. For instance, (a) some important physics may be lost in the reduction to the oneband approximation (such as a role for the oxygen degrees of freedom in the pseudogap phase [8]), (b) near degeneracies of competing states seen in the one-band case [7] may in fact be resolved with the additional degrees of freedom of the three-band model, and (c) the three-band model retains the atomic structure of the CuO 2 layer and thus has a direct link to the structure of real materials as well as experimental measurements of orders at the atomic scale. Previously, the three-band Hubbard model has been investigated with several numerical methods, including direct simulations of finite lattices (by exact diagonalization (ED) [9][10][11][12][13][14][15], quantum Monte Carlo (QMC) [15][16][17][18][19][20][21][22][23], density matrix renormalization group (DMRG) [22,[24][25][26], and the random phase approximation [27][28][29]) and via Green's function based embedding theories (such as dynamical mean-field theory (DMFT) and its cluster extensions [30][31][32][33][34][35], and the variational cluster approximation (VCA) [36,37]). However, due to the complexity of the model, unlike in the one-band case, a consensus on much of the physics has yet to be reached.
Over the past few years, density matrix embedding theory (DMET) [38] has emerged as a powerful cluster embedding method. The basic idea of DMET is similar to that of (cellular) DMFT in the sense that they both map an infinite lattice to an impurity model with an environment that can be described by bath degrees of freedom, and the impurity model is self-consistently improved by matching physical quantities between a single-particle lattice solution and the correlated cluster (impurity) calculation. Technically, however, DMET has a different structure to Green's function based embedding methods and is formulated without frequency dependence and with a finite set of bath orbitals (bounded by the number of impurity orbitals). The lack of frequency-dependent quantities means that DMET calculations can utilize efficient ground-state impurity solvers that can typically treat larger clusters than can be addressed by solvers that target the impurity Green's function. DMET has been applied to a wide range of fermionic lattice models [7,[38][39][40][41][42], ab initio chemical Hamiltonians [43][44][45][46][47][48], and nonfermionic systems [49,50], as well as excited states [51,52] and time-dependent problems [53]. For a detailed review of DMET, we refer to Ref. [44].
In earlier work, DMET successfully provided an accurate description of the ground-state orders of the one-band Hubbard model [40], including in the difficult underdoped region [7]. In this work, we therefore attempt to understand the more complicated three-band Hubbard model using DMET. As we shall see, we can use DMET to provide a detailed description of the ground-state phases and orders as a function of doping, including the doping asymmetry and atomic-scale orders that are new to the three-band case. Another complication of the three-band model is the much larger parameter space than the one-band case. We use both existing parametrizations that have been published in the literature, as well as explicitly model the influence of different individual parameters on the orders. Our findings provide insights into the detailed picture of magnetic and superconducting orders that is provided by three-band models.

A. Model parametrization
As a minimal atomic model of the CuO 2 layer in cuprates, the three-band model describes the on-site and nearest-neighbor interactions among the Cu d x 2 −y 2 and O p x , p y orbitals [see Fig. 1(a)]. In the hole representation, the  Hamiltonian reads where · · · denotes nearest neighbors, d ( †) iσ and p ( †) jσ destroy (create) a hole with spin σ (∈ {α, β}) on the Cu d and O p orbitals respectively, n d iσ and n p jσ are the corresponding hole particle-number operators, and the charge transfer gap pd is defined as the orbital energy difference, p − d . Similarly to in the one-band Hubbard model, the hopping term and on-site Coulomb repulsion will be denoted t and U , and the Coulomb interaction between nearest-neighbor p, d orbitals will be denoted V pd . Note that the hopping term involves a phase factor (±1) introduced by the choice of orbital orientation in the basis as shown in Fig. 1 There has been much work to determine the parameters of the three-band model; however, a consensus set does not exist [32,[54][55][56][57]. There has been particular debate about the size of the charge transfer gap pd [32,58].
In this work, we consider four sets of published model parameters, see Table I, as well as the sensitivity of orders to changing these parameters. Note that all parameter sets are given in eV, and thus all energies in this work are reported in units of eV unless otherwise specified. The first three sets include only the most essential terms, i.e., t pd , U d , and pd , and thus we refer to them as minimal parametrizations. When normalized to units of t pd , the other parameters vary within a range of 10%. The fourth set involves all terms in Eq. (1). We refer to this as a full parametrization. In the hole representation, the minimal parametrization is equivalent to the full model with t pp , U p , and V pd set to zero.

Framework
DMET approximates the expectation values in the interacting lattice by those in a quantum impurity model. The impurity model is solved simultaneously with a fictitious noninteracting lattice problem, whose ground state defines the bath sites of the impurity model via a Schmidt 043259-2 decomposition [38,59]. Self-consistency is achieved by matching the one-particle density matrix of the impurity model and noninteracting lattice ground states via a correlation potential u applied to the noninteracting lattice. The basic steps of the DMET self-consistency loop are thus (a) compute the ground state of the noninteracting lattice Hamiltonian with correlation potential u, (b) construct the bath sites and impurity model Hamiltonian, (c) solve for the ground state of the impurity model, and (d) match the one-particle density matrices of the lattice Hamiltonian and impurity model to update u. The cycle ends when the correlation potential u is converged.
In this work, we are interested in both magnetic and superconducting phases. Consequently, the correlation potential takes the form where optimizing over v σ and αβ in the self-consistency procedure allows for formation of spin polarized and singlet superconducting pairing (between two spin channels α and β) order in the lattice and impurity problems. The noninteracting lattice Hamiltonian is then of Bogoliubov-de Gennes form [60]. The corresponding ground-state solution is a mean-field Bardeen-Cooper-Schrieffer (BCS) wave function, and a set of bath orbitals that describes the environment can be constructed from the corresponding generalized density matrix. The detailed formulas for the bath construction are summarized in Appendix A and the integral transformation expressions for the BCS mean field can be found in Refs. [40,61]. These routines have been implemented in LIBDMET [62,63].

Impurity and lattice
We used a 2 × 2 impurity cluster of CuO 2 primitive cells [37] which retains the inversion and four-fold rotation symmetry of the lattice [see Fig. 1(a)]. We embedded the cluster in a 20 × 20 unit-cell (40 × 40 site-length) lattice. We performed DMET calculations for dopings x ranging between −0.8 and 0.8 (negative denotes electron doping and positive denotes hole doping). Unless otherwise specified, we initialized u with an antiferromagnetic guess and a random pairing potential.

Impurity Hamiltonian and solver
The impurity model Hamiltonian was constructed using the noninteracting DMET bath formalism [38,40], and the ground state was determined using a density matrix renormalization group (DMRG) solver [64,65], allowing for particle-number symmetry breaking and spin polarization [40]. During the DMET self-consistent cycle, we used a maximum bond dimension M = 800. Subsequent bond dimension convergence checks were performed using (up to) M = 2000. To minimize entanglement and ensure a small bond dimension M in the ground state, we rotated the impurity Hamiltonian into a basis of split-localized molecular orbitals (MOs) obtained from the self-consistent Hartree-Fock-Bogoliubov (HFB) method, where the occupied and virtual MOs were computed using the PYSCF package [66,67], and the occupied and virtual spaces were subsequently localized separately using the Edmiston-Ruedenberg procedure that maximizes the Coulomb energy of each orbital [68,69]. The standard genetic algorithm implemented in the BLOCK program [69][70][71][72] was used to order the orbitals for the DMRG calculation. The tolerance of the DMRG sweep energy was set to 10 −6 . Convergence checks on the accuracy of DMRG energies are described in Appendix C.

DMET self-consistency
We carried out DMET self-consistency using full impuritybath fitting [40,44], where the cost function measures the least-squares difference between the correlated one-particle density matrix γ corr and the noninteracting lattice density matrix projected to the full impurity problem γ mf , We minimized w using a conjugate gradient (CG) minimizer with line search. Since the gap of the noninteracting lattice model is often small (in the case of doped systems), a finite inverse temperature β = 1000 t pd was used to define the noninteracting density matrix to ensure smooth convergence (see Appendix B for further discussion and expressions for the analytic gradient of the cost function at finite temperature). We matched the particle number on the impurity sites and on the lattice exactly by separately fitting the chemical potential using quadratic interpolation [61]. Direct inversion in the iterative subspace (DIIS) [73,74] was employed to accelerate the overall DMET convergence, using the difference of u between two adjacent iterations as the error vector. We chose the convergence threshold to be 10 −4 in the correlation potential u (per site), which we observed to translate to an energy convergence per site of better than 10 −4 . We further analyze the numerical convergence and error estimates for the DMET self-consistency in Appendix C.

Order parameters
To characterize the doping dependence of the ground state, we define average AFM and d-wave SC order parameters. As usual, the AFM order parameter is chosen as the staggered magnetization, where m d i is the local magnetic moment on a Cu-d orbital, , and the η AFM is the local structure factor, The SC order parameter here is evaluated as the average of the Cu-Cu and O-O d-wave pairing components, 043259-3 where · · · limits the summation such that only the pairing between nearest Cu-d orbitals is taken into account, and similarly · · · involves only the next-nearest coupling between O-p orbitals. The d-wave superconducting structure factor η SC is defined as

III. THE THREE-BAND PHASE DIAGRAM
A. Undoped state

Charge and magnetic moments
We present the order parameters for the undoped state from DMET and from reference calculations and experimental measurements in Table II. As expected, the d orbitals are roughly half-filled and the p orbitals are roughly doubly occupied, with some charge transfer between the two due to hybridization. Comparing the full and minimal parametrizations, in the full parametrization, the Cu site is more strongly occupied by electrons, due to the t pp term which smears out the oxygen charge and effectively transfers it to copper (while the effect of U p is very small; see the discussion in Sec. III B). Unlike on the O site, the spin density on the Cu site is polarized, with a large local magnetic moment, which compares well to the experimental value 0.3 ± 0.025 μ B (0.6 ± 0.05 μ B ) [75], as well as previously computed theoretical moments of 0.29 [36] and 0.31 [36] from VCA. In addition, the magnetic moment in the full parametrization is reduced relative to the minimal parametrizations, because the increased electron density on copper dilutes the polarized spin, while the additional holes on oxygen reduce the strength of the super-exchange-based antiferromagnetic coupling. In fact, the local magnetic moments in the minimal models appear to be too large, while that of the full model is similar to experimental results. However, it is also known from one-band calculations that the magnetic moments are overestimated in 2 × 2 DMET clusters relative to the thermodynamic limit (e.g., by about 25% at U = 6; see Fig. S1 in the Supplemental Material [80]. Assuming similar finite-size errors, then the minimal parametrization may provide reasonable magnetic moments at half-filling in the thermodynamic limit (although this is not necessarily the case under doping; see below).

Band gap
As a simple estimate of the single-particle gap, we also computed the energy gap of the converged DMET noninteracting lattice Hamiltonian (DMET NI gap), i.e., E g = ε CBM − ε VBM , where C(V)BM denotes conduction (valence) band minimum (maximum). Note that although the charge and spin densities in the different parametrizations are generally similar, the DMET NI gap varies more significantly, from 2.2 to 4.4 eV. The Hybertsen and Hanke parameter sets were derived from calculations on La 2 CuO 4 (LCO), where the optical energy gap is variously reported as being in the range 1.5 to 2.0 eV [76][77][78] (note that the optical gap is generally smaller than the fundamental gap). The estimated DMET NI gap of 2.5 and 2.2 eV for the Hybertsen and Hanke full parameter sets, respectively, are thus in reasonable agreement with the experimental gap. However, the minimal Hanke parametrization seriously overestimates the gap. The Martin parameter set, obtained from calculations on finite-sized Cu-O clusters, are all systematically larger than in the other sets, and thus give the largest DMET NI gap. However, since the ratio of parameters in the Martin model remains similar to other parametrizations (and thus give rise to similar charge and spin distributions) this suggests that all energy parameters in the Martin model should simply be simultaneously rescaled downward.

Orbital resolved band structure
Unlike in the one-band Hubbard model, where the insulating gap arises between Hubbard bands, the gap in correlated insulators in the three-band model can arise from both Hubbard and charge-transfer mechanisms. In Fig. 2, we plot the projected electronic band structure and density of states from the DMET noninteracting lattice Hamiltonian, as converged for the fully parametrized Hanke model. The CBM is mainly of Cu d character (upper Hubbard band), while the VBM shows mixed character, dominated somewhat by O-p. The mixed orbital character of the valence bands around the Fermi level is consistent with the Zhang-Rice singlet (ZRS) hypothesis [4], in which hybridization between oxygen and copper orbitals induces superexchange that leads to singlets of O and Cu holes. Further support for the ZRS picture comes from the k-dependent orbital weights; that of Cu-d is greater at the point, while that of O-p is larger at the M point, consistent with earlier model analysis of the ZRS state [81] and results from VCA [36]. In total, these observations indicate that the undoped three-band model ground state is a charge-transfer insulator, with mainly a p-d type energy gap (see Ref. [82] for experimental evidence of the charge-transfer nature of the band gap). The strong k-dependent hybridization clearly poses challenges for numerical downfolding techniques to a one-band picture.
Comparing the DMET NI band structure to the Hartree-Fock mean-field description (also shown in Fig. 2), we find that the HF gap (≈4 eV) is significantly overestimated, and the d-p hybridization is significantly weaker, resulting in a VBM with dominant oxygen p character and very narrow dispersion. Thus, the reduced gap and d-p hybridization, both seen in experiment, are fluctuation-driven phenomena, whose average effect is being captured by the DMET correlation potential u.

Hole-doped phases with standard parametrizations
More interesting ground states, including those with superconducting order, appear under doping. An important difference with the one-band case is the asymmetry of the three-band model under doping. We first focus on the orders that appear under hole doping. Although our calculations are all at zero temperature, we can loosely identify the magnitude of the order parameters with transition temperatures in the phase diagram, thus allowing us to compare them to the experimental phase diagram. In Fig. 3, we plot the AFM and d-wave SC order parameters of the Hanke model as a function of hole doping (Hybertsen and Martin minimal model results are very similar to those of the Hanke minimal model, as shown in Fig. S2 of the Supplemental Material [80]). In the fully parametrized model, we find two different solutions of the DMET self-consistency, labeled solution 1 (obtained from a weakly spin polarized AFM guess) and solution 2 (obtained from a strongly polarized AFM guess).
For all parameter sets, we observe that the AFM order parameter decreases as doping increases, consistent with the general behavior of the cuprate phase diagram [82]. However, for the minimal models, the AFM order persists even up to large dopings (e.g., ≈0.15 at x = 0.3). In interpreting this discrepancy, one complication is that computation is measuring atomic scale local order, while experimental measurements are likely averages over various inhomogeneities (e.g., different orientations of stripes in different layers) which would typically lead to reduced moments. Leaving this aside, however, the overestimation of the computed moment could originate either from the remaining finite-size error in the DMET calculation or from the unphysical nature of the parametrization (e.g., the lack of doping dependence of the parameters). From our earlier work on the one-band Hubbard model [40], we know that DMET calculations using a 2 × 2 impurity (e.g., in the range U/t = 6 − 8) indeed overmagnetize not only at half-filling but also in the doped regime (see Fig. S1 of the Supplemental Material [80]). However, the one-band AFM order nonetheless vanishes at dopings larger than 0.25, more rapidly than what we observe in the minimal parametrized three-band model. In addition, the full parametrization of the three-band model also predicts a more realistic trend for the AFM order at large doping. Taken together, this suggests that the observed persistent AFM order is likely due to the oversimplified minimal model parameters. Although the AFM order in the full model does decrease to zero in the observed doping range, it vanishes between x = 0.2 and 0.3 (more similar to the one-band model). This is beyond the experimental boundary for the pure AFM phase (x < 0.1), but close to the boundary of the pseudogap region [83,84]. Like in the one-band model, we would expect longer wavelength orders (such as striped phases [7,40]) to appear in this region with larger computational clusters.
From Figs. 3 and S2, we see that d-wave superconducting order (coexisting with antiferromagnetism) appears in the phase diagram of all parameter sets. (Discussion of additional pairing orders, as well as comparisons to the one-band model, can be found further below). In the minimal models, the dwave pairing reaches a maximum at around x = 0.15-0.20.
As a result of the overestimation of AFM order discussed above, the minimal models show coexistence of AFM + SC order for all the studied dopings. However, in the full parametrization, the two coexist in the ranges 0.1 to 0.4 (for solution 2) and 0.1 to 0.3 (for solution 1), with d-wave order reaching a maximum near x ≈ 0.30-0.35, somewhat larger than seen in experiments (≈0.15-0.2) [75]. Solutions 1 and 2 coincide for x < 0.2 and x > 0.4 but are distinct in between, reflecting the known competition between orders at intermediate doping [84]; solution 1 is slightly lower in energy and displays significantly stronger superconducting order. Note that it is also possible to converge a paramagnetic SC solution [by constraining the correlation potential in Eq. (2) so that v α = v β and = † ]. In this case, the SC order is already evident at x = 0.1, since the AFM order is artificially suppressed. However, the energy of this paramagnetic state is much higher than the AFM + SC states we have discussed and is unstable if one releases the constraints on the potential. We thus believe the coexistence of AFM and SC order to be a true feature of the three-band model ground state, as has also been observed in VCA studies [36,37].

Range of reasonable parameters
In view of the significant differences between the minimal and full parametrizations, we now examine more deeply how individual parameters influence the phase diagram. To do so, we change the individual parameters appearing in the Hanke minimal model and restrict ourselves to the magnetic order for simplicity. We compare the magnetic phase diagram computed using both Hartree-Fock and DMET. While the mean-field Hartree-Fock method overestimates the magnetic moments, and the resulting AFM domes always lie above the DMET ones in the plots, it should be noted that parametrizations are often derived from mean-field calculations. Thus, the difference in sensitivity between DMET and HF to the model parameters gives some insight into the sizes of errors arising from mean-field parametrization schemes.
We first study the influence of pd , whose value is uncertain in the literature [32,58]; the magnetic phase diagram computed using both Hartree-Fock and DMET is shown in Fig. 4. On the hole-doped side, when pd 4 eV, the HF magnetic moment does not vanish even at a large doping of x ≈ 0.8, while in contrast, DMET always predicts a finite AFM region with a sharp peak at x = 0. The DMET magnetic moment m AFM increases monotonically from 0.14 to 0.44 as we increase pd , which can be understood from second-order perturbation theory: The effective d-d hopping t dd ∝ t 2 pd pd , thus a larger pd gives a smaller t dd , and thus enhances the magnetic moment. Along with the larger moments, the critical doping point where m AFM → 0 shifts to larger doping as pd increases. Given that, even accounting for finite cluster errors (see above), the minimal model appears to overestimate the magnetic moment under doping and the critical doping concentration; these results suggest one should renormalize pd to smaller values, around 2-3 eV. Finally, we see that the asymmetry with respect to electron and hole doping becomes more pronounced when pd increases, and the magnetic moment is less sensitive to hole doping rather than electron doping. Thus, the appropriate value of pd should neither be too small (as the AFM order and doping asymmetry will both be too weak; see also Ref. [32] for a discussion of the unphysical behavior with small pd ) nor too large (m AFM order will be too strong to be suppressed by doping, especially on the hole-doped side). Reference [34] suggests a range (1.2-2.6 eV) of pd for cuprates, which overlaps the range of our estimates.
We next check the effect of on-site Coulomb repulsion terms. The moment versus U d is shown in Fig. 5. Unlike pd , the influence of U d on the shape of the curves is very small, e.g., the undoped DMET m AFM only increases from 0.35 to 0.38 when U d varies from 6 to 14 eV. The influence on the curve shape is more significant for HF than it is for DMET. In the one-band Hubbard model, however, the situation is very different, where m AFM increases substantially as U is increased [40]. This observation supports viewing the three-band model as primarily a charge transfer insulator (and thus less sensitive to the change in the on-site Coulomb U d ), rather than a Mott insulator, whose magnetic moment is directly mediated by U . The situation for the on-site Coulomb repulsion U p (see Fig. 6) is very similar to that for U d : The undoped DMET m AFM only increases from 0.37 to 0.38 as U d varies from 0 to 8 eV, and the HF curves show a similarly weak sensitivity.
We finally study the effect of nearest neighbor oxygen hopping t pp (see Fig. 7). From the figure, we see that the AFM order is effectively frustrated by large t pp , similar to the effect of t in the one-band Hubbard model. It has been shown in Ref. [32] that t pp can vary substantially for different cuprates (unlike t pd , which is almost unchanged between materials). Our results here suggest that a reasonable range for this parameter is around 0.5-1.0 eV; too large a t pp will suppress the AFM order.
Overall, we find that the magnetic phase diagram is sensitive to pd and t pp but not to U d and U p . The improved results of the full model are thus likely due to the introduction of t pp , rather than U p . In particular, if we wish to have a reasonable description of the three-band Hubbard model within a minimal set of parameters, pd should be renormalized to a smaller value to take the effect of t pp into account. The Hanke parametrization of the full model yields more physical results, and thus we will only use this full model in the remainder of the discussion. However, we note that it is still not optimal with respect to choosing values of pd and t pp that match experiment. This may in part be due to the mean-field derivation of some of the parameters.

Electron doped phases in the full model
We now turn to the electron doped orders, which as mentioned above, are different from the hole-doped orders, unlike in the one-band model [26,84]. We show the AFM and SC order versus both hole and electron doping in Fig. 8 (the holedoped side corresponds to solution 1 in Fig. 3). As we dope with more electrons, the AFM order diminishes. The critical doping x c that makes m AFM vanish (0.15-0.20) is smaller than that on the hole-doped side. This is quite different from what is seen in experiments: The commonly accepted cuprate phase diagram typically shows a sudden drop of AFM order on the hole-doped side [82], with a larger region of coexistence on the electron-doped side. This likely reflects the fact that a single parameter set does not describe the electron-doped and hole-doped materials equally well.
For the SC phase, the overall d-wave pairing magnitude is smaller in the electron-doped region, similar to the lower T c s seen in experiment. Also, the SC phase on the electron doped side has an interesting "M"-shaped two-peak structure: The d-wave SC order increases first with respect to the doping but decays to a small value around the AFM critical x c , before growing to another peak after the AFM order vanishes. The first peak around x = 0.05 is very similar in shape to the peak in DMET calculations of the one-band Hubbard model, where the SC order emerges immediately after doping (see the lower panel of Fig. 8). The second peak, occurring after the disappearance of the AFM order, is similar to the hole-doped SC peak. The presence of two qualitatively different SC phases may be a hint of the types of competing orders that can arise on the electron-doped side, which to date have not been much investigated in numerical studies.

Atomic scale orders in the full model
Beyond the bulk order parameters, the three-band model and the explicit inclusion of both copper and oxygen atoms into the DMET impurity cluster allows for the possibility of studying the magnetic and superconducting order at the atomic scale. The explicit charge, spin, and pairing orders are shown in Fig. 9. We only present representative results from the Hanke full model at x = 0.0, x = 0.3 (solution 1) and x = −0.3 doping, since the results from other parametrizations and dopings are qualitatively similar. (Further plots are presented in Figs. S3-S8 of the Supplemental Material [80]). Comparing Figs. 9(a) and 9(b), we see that on doping the holes mainly occupy the oxygen sites and the hole density on copper only increases slightly. Combined with the fact that doped electrons mainly reside on Cu [see Fig. 9(e); the hole density on Cu is reduced], this reflects the particle-hole asymmetry of the three-band model [26,84]. With respect to pairing order, we see d x 2 −y 2 -wave symmetry clearly between neighboring Cu sites (i.e., it transforms according to the B representation of the C 4 group and the sign of the pairing changes on rotating by 90 • ); see Figs. 9(b) and 9(e). The Cu-Cu pairing order is the largest pairing order between the atoms. From Figs. 9(c) and 9(f), we also see d-wave order between the next-nearest O p orbitals. Although the magnitude is slightly smaller than that of the Cu-Cu pairing, it still contributes almost ≈20-40% of the bulk d-wave order in Eq. (6). We note that the O-O pairing contribution is also asymmetric with respect to doping. In particular, its contribution can be as large as ≈40% in the hole-doped side but only 20-30% in the electron-doped region. Finally, we consider the pairing order between Cu-O and the nearest O-O atoms; see Figs. 9(d) and 9(g). We see that the coupling between the nearest O-O atoms has s-wave symmetry but is quite weak, related to the incompatible orbital orientations. On the other hand, we find the pairing between Cu-O to be relatively strong (in all parameter sets). The local symmetry of Cu-O coupling has p x(y) -wave [or d xy(yz) -wave] symmetry (the pattern transforms according to the E representation of the C 4 group), which to our knowledge has not previously been reported. We note that the superconducting phase pattern between Cu and O is similar to the orbital current-current correlation patterns in Ref. [13], although the current-current correlations were reported to be extremely weak. The pattern is also similar to the asymmetry reported as a hidden order in polarized elastic neutron diffraction experiments [8]. Further investigation of these and other intriguing connections to intracell orders is left to future work.

IV. CONCLUSIONS
In summary, we have used density matrix embedding theory to characterize the ground-state phases of the three-band Hubbard model. We have calculated the charge, local magnetic moments, projected energy bands and density of states of the undoped three-band model, which support a chargetransfer insulating character at zero doping.
We also studied the doping dependence of the ground-state (phase diagram) of the model paying particular attention to the local antiferromagnetic (AFM) and superconducting (SC) orders. In a broad range of model parameters we find a decrease in AFM order upon doping and a SC dome. Unlike in the one-band picture, the models all predict a large region of coexistence of AFM + SC orders, with the AFM order decreasing quite slowly. Comparison to experimental data and earlier theoretical studies suggests that the minimal parametrized models overestimate the AFM order and lead to poorer energy gaps, relative to the full parametrizations, which also include oxygen and oxygen-copper Coulomb repulsion, and oxygen-oxygen hopping. The magnetic moment is particularly sensitive to the pd and t pp parameters, and in the minimal model, the charge transfer gap pd should be renormalized downward to better capture the experimental 043259-8 The three-band model further allowed us to study order at the atomic scale. In the SC region, we observed strong dwave pairing between Cu-Cu and the next-nearest O-O, weak extended s-wave coupling between the nearest O-O atoms, and p-(or d xz , d yz )-like symmetry pairing between Cu-O. The intriguing symmetry of the latter order, similar to that seen in some experiments, illustrates the new physics that emerges at atomic length scales in the three-band model. Exploring such physics in more detail will be the subject of future work. Given a Slater determinant lattice mean-field wave function HF , the DMET bath orbitals can be constructed in several equivalent ways, e.g., via a singular value decomposition (SVD) of the MO coefficients [38] or the environment-impurity part of the density matrix [41], or via eigenvalue decomposition of the projected overlap matrix [43] or environment-environment block of the density matrix [44]. All these methods define a set of bath orbitals, which has nonzero overlap with the impurity sites. In this work, we use SVD of the environment-impurity block of the density matrix to efficiently construct the bath in a periodic lattice system [41,47]. By taking the first unit (or super-) cell as the impurity and the remaining cells as the environment, the whole density matrix of the lattice is divided into four blocks, where γ imp−imp is the density matrix of the first cell, i.e., γ (R = 0), while γ env−imp is the coupling density matrix between the first cell and other cells, i.e., γ (R = 0). These two blocks can be easily computed in periodic systems by a Fourier transform of the density matrix in k space, where N k is the number of k points in the first Brillouin zone. The bath can then be computed by SVD of γ R =0 [47],  where B yields the bath orbital coefficients and the singular values measure the entanglement between bath and impurity orbitals. We note that the γ imp−env and γ env−env blocks are not needed for bath construction and in fact, their computation and storage would be prohibitively expensive in a periodic calculation with many k points. Therefore, using the SVD of γ R =0 is more economical [O(N k N 3 orb ) cost] than diagonalizing the γ env−env block.
In addition to the cost, there are two other advantages of using the SVD. First, it is easy to discard the noncoupled bath orbitals, i.e., bath orbitals with (almost) zero singular values. These bath orbitals are essentially core or virtual orbitals, which have little entanglement with the impurity, and they should be removed from the impurity problem for better numerical stability during the DMET self-consistency. Second, when using finite-temperature smearing, the idempotency of the lattice mean-field density matrix is slightly broken. Rigorous treatment of finite temperature requires a large number of bath orbitals than the number of impurity sites [52]. However, the SVD still provides a good first approximation of the finitetemperature bath, especially at low temperatures.
For the superconducting states HFB , the construction of the bath can still be carried out using SVD, but acting on the (env-imp block of) of the generalized density matrix [40], where the anomalous part of the density matrix κ allows for a nonzero SC order parameter. Within the singlet pairing picture, particles of one spin sector are allowed to couple with holes of the other spin. Therefore, the number of bath orbitals is effectively two times larger than in the normal state DMET.

APPENDIX B: ANALYTIC GRADIENTS OF COST FUNCTION EQ. (3) AT FINITE TEMPERATURE
Once the gradients of Eq. (3) are obtained, we can utilize efficient gradient-based numerical methods, such as CG or the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, to optimize the correlation potential. By differentiating Eq. (3) with respect to u i j we have and thus the key task in Eq. (B1) is to evaluate the response of the mean-field density matrix with respect to a perturbation, ∂γ mf kl /∂u i j . The response at zero temperature can be written in terms of orbital coefficients and energies (see, e.g., Refs. [44,61]) using first-order perturbation theory, where we have assumed the system is gapped. However, when the system becomes (nearly) gapless, this expression diverges. In such cases, the divergent gradient causes the optimization to fail, and this is a source of many convergence difficulties in DMET.
One way to ameliorate this issue is to introduce a finitetemperature smearing, similar to what is used in mean-field calculations of metals. With an inverse temperature β and a perturbation δu, the Fermi-Dirac density matrix is defined as where μ is the Fermi level for the (quasi)particles. The response of γ with respect to the correlation potential u then involves two terms, where the first term is the direct response of the density at a fixed Fermi level, while the second term reflects the contribution of the implicit change in the Fermi level due to the change in potential. The final expression for the first term in Eq. (B4) is where It is easy to check that K pq is always finite when ε p = ε q . One can also let β go to infinity and choose p (q) to label occupied (virtual) orbitals; the gradient then gives the correct zero-temperature limit in Eq. (B2) (up to a symmetrization). Usually this contribution is very small at low temperatures, compared to the direct response in Eq. (B5). However, this contribution will be important in a real finite-temperature simulation, e.g., in Ref. [52]. We summarize the derivation of Eqs. (B5)-(B7) in the Supplemental Material [80].

APPENDIX C: NUMERICAL CONVERGENCE
Here we assess the accuracy and convergence of the DMET procedure in the three-band model calculations. The error in the DMET calculations arises from three possible sources: (a) DMET self-consistency error (from incomplete convergence), (b) DMRG solver error due to the finite bond dimension, and (c) error from the finite size of the impurity. The finite-size error (c) can, in principle, be eliminated by increasing the cluster size and extrapolating to the thermodynamic limit (TDL), as performed in the one-band Hubbard model case [40]. In this work, we use a fixed 2 × 2 cluster size due to the increased computational cost of the three-band model, and thus we cannot assess the finite-size error, except via some comparisons to the 2 × 2 cluster error in the one-band model. However, the error due to (a) and (b) can be directly estimated in our framework, which we now discuss.  Figure 10 shows the overall convergence of DMET with respect to the number of DMET self-consistent iterations. We observe qualitatively different convergence in the normal and superconducting parts of the DMET phase diagram. To illustrate this, we plot the DMET energy, AFM, and (d-wave) SC order parameter for the Hybertsen model at different dopings x. (These order parameters are defined precisely in Sec. II B). We first discuss the undoped system. Here we see that the DMET cycle converges smoothly within seven iterations. For the DMET energy, a single DMET step is enough to converge to ≈10 −4 , demonstrating the utility of single-shot DMET calculations in normal (and especially nonmagnetic) states. The order parameters (density matrices) are more strongly affected by self-consistency. We find that the AFM order increases during the iterations, while the SC order is suppressed, giving a pure antiferromagnetic state at convergence. We next consider x = 0.2 doping. Here, the self-consistency cycle converges more slowly, requiring about 20 DMET iterations to reach convergence. The total energy as well as AFM order converges at around the 10th iteration, while the SC order oscillates until the 20th iteration. This in part reflects the influence of the initial guess: the AFM guess [v σ in Eq. (2)] is quite close to the converged potential, while the SC guess [ αβ in Eq. (2)] is initialized randomly and thus needs more iterations to converge. If we were to restrict the DMET optimization to only pairing potentials with d-wave symmetry (as is commonly done in most cluster DMFT [85] or VCA calculations [37]), the convergence would be much faster. However, the more general form of the correlation potential in DMET allows for the possibility of other pairing channels and orders to emerge. The remaining DMET self-consistency error can be estimated from the difference between the expectation values (e.g., DMET energy) of the last two iterations [40], e.g., δE = 1 2 |E (n − 1) − E (n)|. Consistent with our chosen convergence criterion, the typical size of the DMET selfconsistency error is less than 10 −5 (for both the energy and order parameters) in the undoped region and less than 10 −4 (for the energy) and ≈10 −3 (for the order parameters) in the doped region.
The error from the DMRG solver can be estimated using standard techniques based on the discarded weight in the DMRG calculation [69,86,87] and can be further reduced by extrapolation. The error in the impurity observables (used to evaluate the DMET energy and order parameters) is linear in the (sufficiently small) discarded weight δ and hence can be extrapolated to the exact result (δ = 0) [87]. The convergence with bond dimension M for fixed correlation potential u is shown in Fig. 11. We find that the discarded weight in the normal state (undoped model) is extremely small and usually less than 10 −8 , and thus extrapolation is unnecessary. In fact, calculations can be carried out using a bond dimension as small as M = 100 without any significant error. On the other hand, when the system becomes superconducting, the discarded weight also increases, e.g., to 3 × 10 −5 at M = 800, indicating that the system is more entangled. In such situations, extrapolation has a significant effect on the DMET expectation values. Compared to the extrapolated values, at M = 800 the error in the energy (per site) and order parameters is about 10 −3 .
In summary, from the above analysis, we find that the DMET calculations can be smoothly converged, with minimal error from either the self-consistency or from the solver.