Phase diagram and single-particle spectrum of CuO$_2$ layers within a variational cluster approach to the 3-band Hubbard model

We carry out a detailed numerical study of the three-band Hubbard model in the underdoped region both in the hole- as well as in the electron-doped case by means of the variational cluster approach. Both the phase diagram and the low-energy single-particle spectrum are very similar to recent results for the single-band Hubbard model with next-nearest-neighbor hoppings. In particular, we obtain a mixed antiferromagnetic+superconducting phase at low doping with a first-order transition to a pure superconducting phase accompanied by phase separation. In the single-particle spectrum a clear Zhang-Rice singlet band with an incoherent and a coherent part can be seen, in which holes enter upon doping around $(\pi/2,\pi/2)$. The latter is very similar to the coherent quasi-particle band crossing the Fermi surface in the single-band model. Doped electrons go instead into the upper Hubbard band, first filling the regions of the Brillouin zone around $(\pi,0)$. This fact can be related to the enhanced robustness of the antiferromagnetic phase as a function of electron doping compared to hole doping.


Introduction
The three-band Hubbard (3BH) Hamiltonian contains a minimal set of relevant orbitals to describe the physics of the copper-oxide layers in the high-temperature superconducting cuprates, namely the Copper d x 2 −y 2 and the oxygen p x and p y orbitals [1].
By eliminating the oxygen degrees of freedom, one can further approximately reduce this effective Hamiltonian to the single-band Hubbard and to the t − J model [2], which already contain essential ingredients to describe a correlated system. Having less degrees of freedom, these latter models can be more easily treated via numerical methods such as Exact Diagonalisation (ED) or Quantum Monte Carlo (QMC). Moreover, they already describe appropriately a number of relevant physical properties of the high-Tc cuprates (see, for ex., Refs. [3,4,5,6,7]) Nevertheless, because of the approximate reduction, it is not clear whether the 3BH model contains important physics that is not described by a single-band Hubbard (1BH) Hamiltonian. The nature of the insulating state at half filling is fundamentally different in the two models in the relevant parameter range: while the 1BH model describes a Mott insulator, in the 3BH model one has a charge-transfer insulator. More differences can be seen upon doping. In a 3BH model there is a qualitative difference between hole and electron doping: while doped holes go into the oxygen orbitals forming the famous Zhang-Rice singlet [2], doped electrons go into the Cu d x 2 −y 2 orbitals. This is believed to be the reason for the better stability of the antiferromagnetic (AF) phase when doping with electrons with respect to doping with holes: while introducing electrons in the Cu-orbitals merely produces a dilution of the spins, whereas holes on oxygen sites produce a ferromagnetic coupling between neighboring Cu orbitals, which is more effective in destroying the AF order. While it is argued that one can incorporate this intrinsic electron-hole asymmetry of the 3BH model into the effective 1BH model by appropriate hopping terms [3,4], a detailed comparison between the 3BH model and such approximation is, so far, lacking.
Quite generally, the phase diagram of the high-temperature superconductors (HTSC) displays a variety of competing phases at low temperatures or energies, i. e., AF behavior, stripes, pseudogap behavior and d x 2 −y 2 pairing. These phases are nearly degenerate in energy. Therefore, in the theoretical modeling a change such as the occurrence of p − d charge fluctuations, included in the 3BH model but not in the 1BH model, can tip the balance. Furthermore, the balance can be tipped between the competing orders by using different sets of parameters for a given model and/or by using different techniques for solving this model. This is documented in numerical studies already for the 1BH and t-J models, which have shown how delicately balanced these models are between the nearly degenerate phases [5]. For example, altering the next-nearest-neighbor hopping t nnn or the strength of the on-site correlation U can favor d x 2 −y 2 pairing correlations over stripes [6]. Similarly, the delicate balance is also reflected in the different results obtained using different numerical techniques for the same model. For example, the cluster lattice sizes and the boundary conditions may frustrate stripe formation [7]. For the 3BH model early QMC studies have produced a set of model parameters that consistently describe salient features of the HTSC, such as the magnitude of the charge-transfer gap and the T -and doping dependence of the normal-state magnetic susceptibility [8,9]. On the other hand, the well-known minus sign problem embedded in the QMC calculations did not allow for simulations searching, in particular, for a superconducting (SC) state in the low-temperature limit for arbitrary fillings in the Hubbard models.
Cluster techniques provide a controlled way to systematically approach the infinitesize (and, thereby, low-energy) limit. Recently, progress has been obtained with the variational cluster approach (VCA), which was proposed and used by Potthoff and our group [10,11]. This approach provides a rather general and controlled way to go to the infinite-sized lattice fermion system at low temperatures and at T = 0, in particular. The ground-state phase diagram of the two-dimensional 1BH model was calculated within VCA by Senechal et al. [3], and by our group [12,4]. There are technical differences in the two works, but the overall results are similar. For the cluster sizes used (up to 10 sites) in the VCA, the T = 0 phase diagram of the 1BH model correctly reproduces salient features of the HTSC, such as the AF and d-wave SC ground state in doping ranges, which are qualitatively in agreement with both electron-and hole-doped cuprates. So it appears that much of the difference between electron-and hole-doped cuprates, in particular, the different stability of the AF phase, can be accounted for by a simple 1BH model, in which particle-hole symmetry is broken by a next-nearest-neighbor hopping term, i. e. t nnn .
A central question to be studied in the present paper is what is the role of the oxygen degrees of freedom, of p − d charge fluctuations etc. on the T = 0 phase diagram and the corresponding single-particle excitations. If the results of this 3BH study are similar to the above 1BH model case, one is tempted to further argue in favor of the dominant role played by spin fluctuations (which to a first approximation are equally well captured in both models) over the role of charge fluctuations. Up to what doping does this hold?
Correspondingly, in this paper, we carry out a detailed calculation of the phase diagram and of the single-particle spectrum of the 3BH model from the hole-to the electron-doped region. We use the above recently proposed cluster variational method which is appropriate to treat strongly-correlated systems, i. e., the VCA [10,11]. We allow both for a AF and for a SC phase, as well as for a mixed AF+SC one, as obtained for the 1BH model [3,12]. In the different doping range from underdoped to optimally and overdoped in both electron and hole-doped regions we evaluate the single-particle spectrum. The main outcome of our calculations is that indeed the physics of the 3BH and of the 1BH model is very similar concerning both the phase diagram as well as the single-particle spectrum. More specifically, at low doping we obtain a mixed AF+SC phase in which the transition to the pure SC phase is first order as a function of chemical potential µ h , and, therefore, accompanied by phase separation. The critical doping at which the AF+SC phase disappears is qualitatively consistent with experimental results being larger in the electron-doped case. Also the single-particle spectrum in the vicinity of the gap looks very similar to the one obtained for the 1BH model, whereby the lower Hubbard band is replaced by the dispersive Zhang-Rice band.
Related recent works address the asymmetry between electron-and hole-doping in the three-band Hubbard model without oxygen-oxygen interaction: In Ref. [13] an asymmetry between the quasiparticle mass enhancement in the electron and hole-doped region was found, which however gets compensated in the infrared optical spectral weight. In Ref. [14], the asymmetry between electron and hole doping, and the presence of waterfall structures was discussed within a Local Density Approximation (LDA)+Dynamical Mean Field Theory approach to the three-band Hubbard model. Ref. [15] contains a combined local density-functional theory and DCA (dynamical cluster approach) for 3BH-models for hole-doped cupraes. One of the important outcomes of this study is that the occurrence of SC transition depends rather sensitively on the "down-folding", i.e. whether the orbital basis set is more or less localized (with the latter reproducing the LDA bands over a larger energy window). A very relevant paper, similar in spirit albeit at finite temperature, is the work by Macridin et. al. [16], addressing the low-energy physics of cuprates within a two-band Hubbard model and again checking for the validity of the one-band Hubbard model. Here, only oxygen states that hybridize directly with Cu are considered. Nevertheless, the conclusion is similar, in that a t − t ′ − U single-band Hubbard model captures the basic low energy (below ∼ 0.5eV ) physics. For further work on the 3BH model see also, e. g., Refs. [17,18,19,20,21,22,23,13,14,15].
Our paper is organized as follows: In Sec. 2 the model we are using is presented and the VCA method summarized. Results for the phase diagram and for the single-particle spectrum are shown in Sec. 3. Finally, in Sec. 4 we draw our conclusions.

Model and method
It is commonly believed that the relevant physics of high-Tc superconductors takes place mainly on the copper-oxide layers. A minimal model for these layers contains the d x 2 −y 2 orbital as well as a p x and p y orbital per unit cell. The electron dynamics in these orbitals is described by the three-band (Emery) model, which, in hole notation, reads [1] Here, d † Rσ and p † Sσ create a hole in the copper 3d x 2 −y 2 and in the oxygen 2p δ orbital (δ = x, y depending on the position with respect to the Cu site), respectively, and n d Rσ , n p Sσ are the corresponding occupations. R denote copper and S oxygen lattice sites. RS limits the sum over next-neighbor Cu-O and SS ′ over next-neighbor O-O lattice sites. ǫ d and ǫ p are the on-site energies of the copper and the oxygen orbitals, respectively, and µ h is the hole chemical potential. t pd is the Cu-O and t pp is the direct O-O hopping amplitude. We set the zero of the energy such that we have ǫ d = 0 on the copper site, and introduce the charge-transfer gap ∆ = ǫ p − ǫ d . The Coulomb repulsion between two holes is taken into account by the terms U dd , U pp , and U pd , when the two holes sit on the same Cu orbital, on the same Oxygen orbital, or on two neighboring Cu and O orbitals, respectively. α RSδ and α ′ SS ′ describe the sign of the phases of the Cu-O and O-O hopping due to the d and p symmetries of the Cu and O orbitals, respectively (according to the convention in Fig. 1). In our calculation, we take typical values of the parameters, as obtained by constrained density-functional calculations [24], and consistent with earlier cluster calculations [25], as well as extensive Quantum-Monte-Carlo studies [8,9]. In units of the Cu-O hopping t pd ≡ t, we thus have t pp = 0.5, ∆ = 3.0, U dd = 8.0, U pd = 0.5, and U pp = 3.0.
Numerically exact methods such as Quantum Monte Carlo or exact diagonalisation are restricted to relatively high temperatures or small cluster sizes [17]. The self-energyfunctional approach (SFA) [26] provides a variational scheme to use the information from the solution of an exactly solvable "reference system" (for example a small cluster) to approximate the properties of an infinite lattice. With the help of this reference system, which must have the same interaction as the original one (the infinite lattice), one can evaluate the grand potential Ω of the original system exactly within a restricted set of self-energies. The "best" solution thus is given by finding the saddle point of Ω within this set of self-energies. Within the variational-cluster approach [10,11], the reference system is obtained by splitting the original lattice into small clusters and by adding single-particle t ′ terms which parametrize the variation of the self-energy. The SFA potential Ω of the original system can be evaluated exactly for the self energies Σ accessible to the reference systems and is given by Here, G 0 is the free Green's function of the model given by Eq. (1), Ω ′ , Σ, and G ′ are the grand canonical potential, the self-energy and the Green's function of the cluster reference system which depends on the the single-particle parameters t ′ . In the present study we consider the 12-sites cluster (i. e., a 2×2 unit cell) shown in Fig. 1 as a reference system and search for the stationary solution characterized by the condition ∂Ω/∂t ′ = 0. This stationary point provides a good approximation to the exact solution for the system in the thermodynamical limit provided the self-energy is sufficiently "short ranged", i. e., sufficiently localized within the cluster. In the limit in which the self-energy is exactly localized within the cluster, the VCA becomes an exact solution.
Since we want to describe the transition from the AF phase at low doping to the SC phase at higher doping, the reference system used in the VCA consists of the 2 × 2 CuO 2 unit-cell, displayed in Fig. 1 Here, r = ±e x , ±e y (we take the size of the unit cell a = 1), and the d-wave phase factor η d r = +1 for r = ±e x , and η d r = −1 for r = ±e y . The AF order is induced by η AF R = exp(iQ · R) with Q = (π, π). The fields controlling the variation of the selfenergy in the reference system are the AF (h AF ) and the the d-wave SC (h SC ) "Weiss" fields [3,12], as well as a shift in the on-site energy (∆ǫ), which is required in order to describe consistently the particle density [4]. We stress that these Weiss fields are only present in the reference system and are used in order to optimize the self-energy to be used in the original lattice model. In principle, one could introduce additional variational parameters such as a different on-site energy shift and a different SC Weiss field for d and p orbitals. However, too many variational parameters make the numerical search for the saddle point of Ω too difficult and time consuming.
Notice that the intercluster interaction terms cannot be treated exactly by VCA, so that-in principle-we are not able to treat U pd interactions beyond the cluster of Fig. 1. However, to neglect them completely would be a quite rough approximation, since, for example, O orbitals at the cluster boundary would be influenced by U pd much less than the ones in the center. In order to overcome this problem, we treat nonlocal interactions by periodic boundary conditions. This turned out to be quite accurate, whenever the system is not close to a charge-density wave phase transition [27].

Results
We first determine the zero-temperature (T = 0) phase diagram in the hole-and electron-doped region by determining the optimal values of the variational parameters h AF , h SC and ∆ǫ as a function of the chemical potential µ h . In general, there is always more than one solution for a given µ h : we adopt the usual criterion of taking the solution with the lowest grand-canonical energy Ω. In Fig. 2 we plot the staggered magnetization (red color) the d-wave superconducting order parameter (blue) phase at low doping both in the hole-as well as in the electron-doped case. The AF phase ends at a first order transition as a function of chemical potential, which is accompanied by phase separation. In Ref. [30] it was discussed that this phase separation is likely to be persistent for the 1BH in the hole-doped case, while it seems to disappear when  considering larger clusters in the electron-doped case. For the 3BH studied here, the critical doping at which the AF phase is destroyed is 3.5% in the hole-doped and 6.5% in the electron-doped case. This is again in qualitative agreement with the experimental situation, in that the AF phase is more stable for electron doping. It is remarkable to observe such similarities between the 1BH and the 3BH model (cf. Refs. [3,12]) despite of the fact that doping, and -in particular -the asymmetry between electron-and hole-doping is fundamentally different in the two models.
In order to explore the relation between single-particle excitations, their doping evolution, and the phase diagram of Figs. 2 and 3, we plot in Figs. 4 and 5 the spectral function A( k, ω) of the 3BH model (in a grayscale plot) for different dopings in both the hole- (Fig. 4) and electron-doped ( 5) case.
The spectrum in Fig. 4 shows the well-known features of the 3BH model. Starting  from the top of the spectrum, we first find the upper Hubbard band at ω/t ≈ 2.0. Around the Fermi energy we can see the Zhang-Rice singlet (ZRS) band [2], with a total width of about 1 (in units of t, see also Ref. [31]). This band describes a dispersive singlet state (ZRS) formed between a hole on the plaquette of the O orbitals and one on the Cu site. It consists of a coherent part crossing the Fermi energy and of an incoherent part at higher binding energies. The separation between these two structures can be associated with the so-called "waterfall" structure observed in several high-Tc compounds [32,33,34,35], see also Ref. [14]. Between the upper Hubbard band and the ZRS there is an optical gap of about 1.5t, which is consistent with the experimental value [36]. At higher binding energies ω/t < −3 we find the two oxygen bands: a weakly dispersive "almost-non-bonding" band (the weak dispersion is due to t pp ), and the dispersive antibonding band. The lower Hubbard band is located approximately at the same binding energy. At low energies, the similarity with the spectrum of the 1BH model is remarkable, as can be seen by comparing figures 4 and 5 with the spectrum of the 1BH model plotted in Fig. 6 for comparison (see also Refs. [3,12,4]). Notice that when comparing the band structure between the 1BH and the 3BH models, the ZRS band takes here the role of the coherent quasi-particle band in the 1BH model. On the other hand, the upper Hubbard band of the 1BH model is replaced with the upper Hubbard band in the 3BH model. In the hole-doped case (Fig. 4), holes first enter the ZRS band around (π/2, π/2), and form hole pockets at low doping. Here, metallic screening from nodal electrons makes the AF solution quickly unstable. We find a different situation for electron doping: The additional electrons are first doped around (π, 0)in the upper Hubbard band (Fig. 5). Here, the AF solution remains more stable, since there is a SC gap with maximal size at the anti-nodal point. Thus the screening is less effective, until -as a function of doping electrons start to reach the (π/2, π/2) point. Again these results are consistent wit the 1BH model, and also with experiments (for a review see Ref. [37]).
In order to analyze the evolution of the Fermi surface (FS) as a function of doping, we plot in Fig. 7 and Fig. 8 the low-energy spectrum in the Brillouin zone for the holeand electron-doped case, respectively. This is obtained by integrating the spectrum within an energy window of width 0.2t around the Fermi energy. As discussed above, in the hole-doped case, hole pockets start to build around (π/2, π/2)at low doping, while a large FS centered around (π, π) develops at larger doping. This is in qualitative agreement with the experiments (see, e.g. Ref. [37]) and with previous results on the 1BH model. In contrast, in the electron-doped case (Fig. 8), hole pockets first build around (π, 0), and the FS becomes connected and again closes around (π, π) at larger dopings.

Summary and Conclusions
In this paper, we have carried out an analysis of the three-band Hubbard model in the hole-and in the electron-doped region by means of the Variational Cluster Approach, whereby we allowed for an antiferromagnetic, and a superconducting Weiss field, as well as for a on-site energy as variational parameters for the cluster reference system. Results for the single-particle spectrum and for the phase diagram are qualitatively, and to some extent even quantitatively, very similar to the ones obtained for the oneband Hubbard model with a next-nearest-neighbor hopping included in order to break particle-hole symmetry [3,12,4]. Concerning the phase diagram, we obtained a mixed AF+SC phase at low doping with a transition to a pure SC phase accompanied by phase separation. On the basis of these results we can confirm, as already argued in previous work, that low-energy single-particle properties of the three-band Hubbard model can be quite well captured by a single-band Hubbard model with appropriate longer-range hopping parameters both in the hole as well as in the electron-doped case. Of course, we may expect that two particle excitations, especially charge-transfer but also spin susceptibilities may behave differently in the 1BH and the 3BH models. Investigations along this line are in progress, by means of a recently developed extension of the VCA to treat dynamical two-particle correlation functions [38].