Unified description of cuprate superconductors using four-band $d$-$p$ model

In the 35 years since the discovery of cuprate superconductors, we have not yet reached a unified understanding of their properties, including their material dependence of the superconducting transition temperature $T_{\text{c}}$. The preceding theoretical and experimental studies have provided an overall picture of the phase diagram, and some important parameters for the $T_{\text{c}}$, such as the contribution of the Cu $d_{z^2}$ orbital to the Fermi surface and the site-energy difference $\Delta_{dp}$ between the Cu $d_{x^2-y^2}$ and O $p$ orbitals. However, they are somewhat empirical and limited in scope, always including exceptions, and do not provide a comprehensive view of the series of cuprates. Here we propose a four-band $d$-$p$ model as a minimal model to study material dependence in cuprates. Using the variational Monte Carlo method, we theoretically investigate the phase diagram for the La$_2$CuO$_4$ and HgBa$_2$CuO$_4$ systems and the correlation between the key parameters and the superconductivity. Our results comprehensively account for the empirical correlation between $T_{\text{c}}$ and model parameters, and thus can provide a guideline for new material design. We also show that the effect of the nearest-neighbor $d$-$d$ Coulomb interaction $V_{dd}$ is actually quite important for the stability of superconductivity and phase competition.

In the 35 years since the discovery of cuprate superconductors, we have not yet reached a unified understanding of their properties, including their material dependence of the superconducting transition temperature Tc. The preceding theoretical and experimental studies have provided an overall picture of the phase diagram, and some important parameters for the Tc, such as the contribution of the Cu d z 2 orbital to the Fermi surface and the site-energy difference ∆ dp between the Cu d x 2 −y 2 and O p orbitals. However, they are somewhat empirical and limited in scope, always including exceptions, and do not provide a comprehensive view of the series of cuprates. Here we propose a four-band d-p model as a minimal model to study material dependence in cuprates. Using the variational Monte Carlo method, we theoretically investigate the phase diagram for the La2CuO4 and HgBa2CuO4 systems and the correlation between the key parameters and the superconductivity. Our results comprehensively account for the empirical correlation between Tc and model parameters, and thus can provide a guideline for new material design. We also show that the effect of the nearest-neighbor d-d Coulomb interaction V dd is actually quite important for the stability of superconductivity and phase competition.

I. INTRODUCTION
The discovery of superconductivity in cuprates has brought about the significant progress in stronglycorrelated electron systems [1]. Along with the high superconducting transition temperature T c , cuprates show anomalous phases and phenomena such as the Mott metal-insulator transition, pseudogap phenomena, and antiferromagnetic (AF) and charge-density-wave (or stripe) phases [2,3]. Recently, nematic order [4] and ferromagnetic fluctuation [5,6] have also been observed experimentally. These newly observed experiments being constantly reported with the underlying, possibly exotic, physics have continued to attract many researchers' interests over 35 years since the first high-T c cuprate superconductor was synthesized.
It has been widely believed that the strong AF spin fluctuation characteristic to the two-dimensional square lattice structure triggers the various anomalous phenomena in cuprates [7][8][9][10]. On the basis of this picture, the effective one-band Hubbard model [11] has been studied with various numerical methods extensively [12,13]. It * h-watanb@fc.ritsumei.ac.jp succeeded to describe the important physics in cuprates not only for the ground state but also for the excited state. The angle-resolved photoemission spectroscopy (ARPES) experiments show that the Fermi surface consists of only one energy band mainly derived from the Cu d x 2 −y 2 orbital [14] and thus the one-band models are justified as long as the low-energy physics is concerned.
However, recent theoretical and experimental development sheds light on the importance of the orbital degree of freedom in cuprates. For example, material dependence of T c is difficult to discuss within the one-band models. The previous studies of the one-band models generally discussed the material dependence via the shape of different Fermi surfaces originated by different hopping integrals, e.g., the nearest-neighbor hopping t and the next nearest-neighbor hopping t ′ . Indeed, such treatment allowed us to capture the overall feature for the difference between hole-doped and electron-doped systems [15]. However, the material dependence of T c for the hole-doped systems has not been able to be explained properly. Experimentally, T c becomes higher with larger |t ′ /t|, i.e., more rounded Fermi surface [16]. On the contrary, it is predicted theoretically that in most one-band models, T c is reduced because of poor nesting properties in rounded Fermi surfaces [17,18]. A clue to resolving this contradiction can be found by seriously considering the orbital degrees of freedom, especially, the Cu d z 2 orbital [18][19][20]. The study for the two-orbital model shows that the contribution of the d z 2 orbital to the Fermi surface, which is strong in low T c materials, works against the d x 2 −y 2 -wave superconductivity [18]. This picture can explain the material dependence of T c . Indeed, a recent ARPES experiment directly observed the hybridization gap between the d z 2 and d x 2 −y 2 orbitals and pointed out the importance of the multiorbital effect on superconductivity [21].
Another extension to a multiorbital model is the introduction of the O p orbitals. A minimal model that includes the Cu d x 2 −y 2 , O p x , and O p y orbitals in the CuO 2 plane is known as the (three-band) d-p model or the Emery model [22]. The three-band d-p model has been studied with various numerical methods  as much as the one-band models. In these studies, the siteenergy difference ∆ dp between the Cu d x 2 −y 2 and O p x/y orbitals is found to be an important parameter for understanding the material dependence of T c . For example, the three-band d-p model can successfully reproduce the negative correlation between ∆ dp and T c [34], namely, a smaller ∆ dp leads to a higher T c . Furthermore, the role of the O p orbitals in the superconducting pairing [37,41,43,44,47], stripe order [26,38], loop current [29,30,33], and nematic order [33,36,39,40] has been discussed with the three-band d-p model. However, it is still difficult to correctly capture the correlation between Fermi surface topology and T c within the threeband d-p model [28,34].
In this paper, we study a four-band d-p model composed of the Cu d x 2 −y 2 and d z 2 orbitals and the O p x and p y orbitals, which can be expected to resolve the problems mentioned above. As typical examples of singlelayer hole-doped cuprates, we consider La 2 CuO 4 and HgBa 2 CuO 4 systems. We construct the tight-binding model for these systems based on the first-principles calculation and examine the effect of Coulomb interaction by the variational Monte Carlo (VMC) method. We show that the material dependence of T c is well explained with two key parameters, the site energy ε d z 2 of the Cu d z 2 orbital and the site-energy difference ∆ dp between the Cu d x 2 −y 2 and O p x/y orbitals, which is consistent with the empirical relation. We thus propose that the present four-band d-p model is a minimal model that can properly describe the material dependence of cuprate superconductors. Furthermore, we also study the effect of the nearest-neighbor d-d Coulomb interaction V dd , which has not been discussed in detail previously. We show that V dd substantially affects the stability of superconductivity and the phase competition among various phases, suggesting an important parameter for the effective model of cuprates. In addition, we find that V dd induces a crossover from a Slater insulator to a Mott insulator at the undoped limit.
The rest of this paper is organized as follows. In Sec. II, a four-band d-p model on the two-dimensional square lat-tice is introduced. The VMC method and the variational wave functions are also explained in Sec. II. The numerical results are then provided in Sec. III. The tight-binding energy bands for the La 2 CuO 4 and HgBa 2 CuO 4 systems, obtained on the basis of the first-principles calculation, are first shown in Sec. III A. The material and doping dependences of a superconducting correlation function are then examined and the role of two key parameters ε d z 2 and ∆ dp are clarified in Sec. III B. The effect of the nearest-neighbor d-d Coulomb interaction V dd is investigated in Sec. III C. The phase competition among superconductivity and other phases is briefly discussed in Sec. III D. Finally, the paper concludes with a summary in Sec. IV. The details of the variational wave functions are described in Appendix.

II. MODEL AND METHOD
A. Four-band d-p model As an effective low-energy model of cuprates, we consider a four-band d-p model on the two-dimensional square lattice (see Fig. 1) defined by the following Hamiltonian: Here, the kinetic term H kin is described by where Eq. (2) is the kinetic term in an orbital representation and Eq. (3) is in a band representation. c † iασ (c iασ ) is a creation (annihilation) operator of an electron at site i with spin σ (=↑, ↓) and orbital α (= 1, 2, 3, 4) corresponding to (d x 2 −y 2 , d z 2 , p x , p y ), respectively. t αβ ij denotes a hopping integral between orbital α at site i and orbital β at site j. t αα ii is a site energy ε α for orbital α at site i. Equation (3) is obtained by diagonalizing Eq. (2), and the energy eigenvalue E m (k) is characterized by the wave vector k and the energy band index m (= 1, 2, 3, 4). a † kmσ (a kmσ ) is a creation (annihilation) operator of the corresponding energy band with spin σ. The undoped parent compounds of cuprates correspond to one hole (i.e., seven electrons) per unit cell in this model, which is conventionally referred to as half filling, and hereafter we denote the carrier density as the number δ of holes per unit cell that are introduced into the system at half filling.
The Coulomb interaction term H int is composed of eight terms, Here, n α i = n α i↑ + n α i↓ with n α iσ = c † iασ c iασ is the number operator and S α i is the spin angular momentum operator at site i with orbital α. d 1 and d 2 are abbreviations for d x 2 −y 2 and d z 2 orbitals, respectively, and and J ′ represent on-site intraorbital, interorbital, Hund's coupling, and pair-hopping interactions between d orbitals, respectively. In this study, we set J ′ = J and U d = U ′ d + 2J [48]. The on-site intraorbital Coulomb interaction within p orbitals, U p , is also introduced. The last three terms in H int take into account the intersite Coulomb interactions between nearest-neighbor orbitals, V dp , V pp , and V dd , as shown in Fig. 2, where the sum i,j runs over all pairs of nearest-neighbor orbitals located at site i and j.
In addition, the following double counting correction term H dc is introduced, are the average electron density of the d and p orbitals in the noninteracting limit and N S is the total number of unit cells. When we apply a many-body calculation method to a multiorbital model, the site energy of each orbital is shifted due to the interaction effect. However, such energy shifts have already been included in the energy band of the tight-binding model constructed from the first-principles calculation. This is a so-called double counting problem and should be treated with care especially in the d-p model [49]. Here, we subtract the term H dc from the Hamiltonian to correct the energy shift. This is one of the reasonable treatments to avoid the double counting.

B. VMC method
The effect of Coulomb interaction is treated using a VMC method [50][51][52]. The trial wave function considered here is a Gutzwiller-Jastrow type composed of four parts, |Φ is a one-body part constructed by diagonalizing the one-body Hamiltonian including the off-diagonal elements {ρ α i }, {m α i }, and {∆ αβ ij }, which induce long-range ordering of charge, spin, and superconductivity, respectively. The renormalized hopping integrals {t αβ ij } are also included in |Φ as variational parameters. The explicit forms of them are described in Appendix.
The Gutzwiller factor TABLE I. The tight-binding parameters for the La2CuO4 and HgBa2CuO4 systems in eV units estimated on the basis of maximally localized Wannier orbitals from the first-principles LDA calculation. For comparison, the tight-binding parameters for the La2CuO4 system with reference to the QSGW method are also shown (denoted as "revised"). The definitions of ti and εα are described in Appendix 1. ∆ dp = ε d x 2 −y 2 − εp, i.e., the site-energy difference between the Cu d x 2 −y 2 and O p x/y orbitals.  Table I. Fermi energy (defined as zero energy) is set to be the case of 15% hole doping (δ=0.15). The high symmetric momenta are indicated as Γ = (0, 0), M = (π, π), and X = (π, 0).
The remaining operators and are charge and spin Jastrow factors, which control longrange charge and spin correlations, respectively. s z iα is the z component of the spin angular momentum operator at site i with orbital α. We set v c iiαβ = v s iiαβ = 0 for α, β = d 1 , d 2 because the on-site correlation between the d x 2 −y 2 and d z 2 orbitals are already taken into account in P G . In this paper, we focus mainly on the superconducting correlation functions to examine where the superconductivity appears in the phase diagram. The detailed studies on other competing orders will be discussed elsewhere. The variational parameters in |Ψ are therefore They are simultaneously optimized using stochastic reconfiguration method [55]. We show results for N S =24×24=576 unit cells (and thus 576×4=2304 orbitals in total), which is large enough to avoid finite size effects. The antiperiodic boundary conditions are set for both x and y directions of the primitive lattice vectors.

A. Band structures of La2CuO4 and HgBa2CuO4
First, we discuss the material dependence of the band structure. As a typical example of single-layer holedoped cuprates, we study the La 2 CuO 4 and HgBa 2 CuO 4 systems. We construct maximally localized Wannier orbitals [56,57] from the first-principles calculation in the local-density approximation (LDA) with ecalj package [58] and fit them with the hopping integrals t i (i = 1 − 6) and the site energy of each orbital ε α . The parameter sets determined for these systems are listed in Table I and the explicit form of the tight-binding model is described in Appendix 1. Note that the estimated site energy of the Cu d z 2 orbital, ε d z 2 , and the hybridization between the Cu d z 2 and O p x/y orbitals, t 4 , depend significantly on the method of first-principles calculation. In fact, the estimated ε d z 2 is much lower in the quasiparticle self-consistent GW (QSGW) method [61][62][63] than in the conventional LDA calculation [64]. To clarify the effect, we also consider another parameter set of La system with reference to the QSGW band structure (labeled as "revised"), where the site energy ε d z 2 is lower and t 4 is also slightly smaller than those estimated on the basis of the LDA calculation (See Table I). Figure 3 shows the noninteracting tight-binding energy bands for La and Hg systems. We can notice the clear difference among them: (i) The density of states (DOS) of the d z 2 component is extended from 0 to -2 eV in the La system [ Fig. 3(a)], while it is almost localized around -2 eV in the Hg system [ Fig. 3(c)]. This is because the d z 2 orbital is hybridized with the p x/y orbital in the La system much more strongly than in the Hg system. The d z 2 electrons obtain the itinerancy through the hybridization and therefore the d z 2 component of the La system is more extended than that of the Hg system. The revised version of the La system is located somewhere in between [ Fig. 3(b)]. (ii) The site-energy difference between the d x 2 −y 2 and p x/y orbitals, ∆ dp = ε d x 2 −y 2 − ε p , is larger in the La system than in the Hg system. This affects the occupancy of each orbital and thus the strength of the electron correlation. For example, when ∆ dp (> 0) is small, the energy band crossing the Fermi energy contains more component of the p x/y orbital, in which the intraorbital Coulomb interaction is smaller.
Starting from these energy band structures, we shall investigate the ground state property of the La and Hg systems using the VMC method. We assume that the tight-binding parameters remain unchanged with hole doping. The Coulomb interaction parameters are set as (U d , U ′ d , J, U p , V dp , V pp ) = (8.0, 6.4, 0.8, 4.0, 2.0, 1.6) t 1 for both La and Hg systems with reference to Ref. [65]. Note that the values for the La system are larger in eV units because of the larger t 1 . In the following, we set t 1 as a unit of energy. We first set V dd = 0 and then discuss the effect of finite V dd . As shown in Sec. III C, we find that even small V dd can substantially affect the property of the system.

Overview
To discuss the material dependence of superconductivity, we calculate the superconducting correlation function defined as where ∆ † τ (R i ) is a creation operator of singlet pairs between nearest-neighbor d x 2 −y 2 orbitals, and τ runs over four nearest-neighbor Cu sites (τ = ±e x , ±e y ). f (dd) τ τ ′ is a form factor of a superconducting gap function with d x 2 −y 2 symmetry, namely, f (dd) τ τ ′ = 1 for τ τ ′ and −1 for τ ⊥ τ ′ . · · · denotes Ψ| · · · |Ψ / Ψ|Ψ for the optimized variational wave function |Ψ . If P dd (r) is saturated to a finite value for r = |r| → ∞, superconducting long-range order exists. Figure 4(a) shows the behavior of P dd (r = |r|) for the Hg system at a hole doping rate δ = 0.153. It shows good convergence for r 4 and reveals that the superconducting long-range order certainly exists. The sign of the d-d pairing in a real space is shown in Fig. 4(b). It is positive in the x direction and negative in the y direction, reflecting the d x 2 −y 2 -wave symmetry expected in cuprate superconductors. We also calculate the superconducting correlation function for pairing formed between the d x 2 −y 2 and p x/y orbitals P dp (r), which is defined in the same way as P dd (r) except that c † i+τ 1σ in Eq. (11) is replaced with c † i+ τ 2 3(4)σ for τ = ±e x(y) . Although the value is one order of magnitude smaller than P dd (r) [see Fig. 4(a)], P dp (r) is also saturated to a finite value, indicative of the long-range order of d-p pairing. Note that the sign of the d-p pairing changes alternatively along both x and y directions and thus shows p-like symmetry as shown in Fig. 4(c). This is due to the phase convention of p x and p y orbitals adopted here [see Fig. 1(b)] and is consistent with the d x 2 −y 2 pairing symmetry [66]. It is also compatible with the preceding study of the threeband d-p model [45]. This kind of real space or orbital representation is useful for the analysis of cuprate superconductors [43,44,47] because the pairing length is expected to be very short [67].

Material dependence: Effect of d z 2 orbital
Let us now examine the material and doping dependence of the superconducting correlation function. We take the converged value of P dd (r → ∞) as a strength of superconductivity P dd . Figure 5 shows the doping dependence of P dd for the La, La(revised), and Hg systems. For all cases, P dd displays a dome shape as a function of the hole doping rate δ. At δ = 0, the system is insulating due to the strong correlation effect and thus P dd = 0. As δ increases, mobile carriers are introduced into the system and the mobility of the Cooper pair increases. On the other hand, the strength of the d-d pairing itself is reduced by doping because the electron correlation effect is also reduced. The balance between these two factors results in the dome-shaped behavior of P dd . This picture is expected to be universal for all hole-doped cuprate superconductors.
Next, let us study the importance of the orbital character near the Fermi energy in the material dependence of superconductivity. Generally, a large DOS at the Fermi energy is favorable for superconductivity due to the large energy gain of gap opening. However, a detailed structure of the DOS, namely, the orbital character and k dependence should be carefully investigated. For this purpose, we calculate the following momentum distribution function of holes, where c † kασ (c kασ ) is a Fourier transform of c † iασ (c iασ ) in Eq. (2). Figure 6 shows n α (k) for the La system at the superconducting phase (δ = 0.153) and the paramagnetic metallic phase (δ = 0.181). The discontinuities around the X point and in the middle of the Γ-M line in Fig. 6(b) for the paramagnetic metallic phase indicate the existence of Fermi surface. Since the node of the superconducting gap runs along the Γ-M line, a discontinuity also exists in the superconducting phase, as shown in Fig. 6(a). We can observe that the d z 2 component has large weight around the X point, exhibiting a peak structure in n α (k), where the superconducting gap becomes the largest (so-called "hot spot"). This feature is expected to be unfavorable for superconductivity because the AF spin fluctuation, which promotes the d x 2 −y 2 -wave pairing, is suppressed when the Cu d z 2 orbital contributes to the formation of the Fermi surface [68].
To investigate the effect of the d z 2 component, the ratio of the d z 2 component to the total at the X point, R = n d z 2 (X)/n tot (X), is calculated in Fig 7(a). For the La system, the ratio R increases with increasing δ and shows rapid enhancement for δ > 0. 16. It coincides with the sudden disappearance of P dd for δ > 0.16 shown in Fig. 5. On the other hand, R for the La(revised) system is much smaller than that of the La system, because ε d z 2 is lower and the hybridization between the d z 2 and d x 2 −y 2 orbitals via the p x/y orbitals is smaller. As a result, the superconducting phase is extended to a larger value of δ and a smooth dome shape is observed in P dd vs. δ as shown in Fig. 5. For the Hg system, R is much more suppressed because the d z 2 -orbital based band is almost localized and detached from the d x 2 −y 2 -orbital based band [see Fig. 3(c)]. This is an ideal condition for superconductivity [18] and therefore P dd becomes largest for the Hg system.
These results conclude that superconductivity is more enhanced when the d z 2 -orbital based band is deeply sinking and its contribution to the low-energy physics is small. Therefore, the material dependence of superconductivity is understood only by incorporating the d z 2 orbital explicitly into a model such as our model, which is a remarkable advantage over the usual one-band Hubbard and t-J models, and even the three-band d-p model. P dd calculated here corresponds to the square of the superconducting order parameter and is closely related to the superconducting transition temperature T c . The critical doping rate δ c for the La system, where P dd becomes zero, seems to be too small (∼ 0.16) compared with the experimental value (δ c = 0.25 − 0.3). Furthermore, the sudden disappearance of P dd is also unrealistic. It can be inferred that the actual value of ε d z 2 is much lower than the value estimated from the LDA calculation. Indeed, as shown in Table I, the value of ε d z 2 with reference to the QSGW calculation is much lower. Furthermore, the QSGW band structure [64] can well explain the resonant inelastic X-ray scattering experiment [69]. We expect that the revised band structure shown in Fig. 3(b) properly includes the correction for the LDA calculation and gives more realistic result for the La system.

Material dependence: Effect of apical oxygen height
From the viewpoint of the actual lattice structure, ε d z 2 is governed by the apical oxygen height, i.e., the distance between the apical oxygen and the copper: The larger the apical oxygen height is, the lower ε d z 2 is with respect to ε d x 2 −y 2 , because of the crystal field effect [68,70]. In addition, a larger apical oxygen height leads to a lower site energy ε pz of the apical oxygen due to the decrease of a crystal field effect, which in turn lowers the ε d z 2 through the hybridization between p z and d z 2 orbitals despite the increase of the distance between apical oxygen and the copper. Therefore, our result suggests that a larger apical oxygen height leads to a higher T c through a lower ε d z 2 . Indeed, the experimentally observed apical oxygen height of the Hg system is larger than that of the La system. This tendency is also consistent with the socalled Maekawa's plot [71], where a lower ε pz is related to a higher T c . Although the model itself does not explicitly include the p z orbital of the apical oxygen, the present four-band d-p model properly incorporates the effect of the apical oxygen height via adjusting the site energy ε d z 2 .

Material dependence: Effect of ∆ dp
We also show in Fig. 7(b) the quasiparticle renormalization factor z estimated from the jumps in the total momentum distribution function n tot (k) along the nodal direction of the d x 2 −y 2 -wave superconducting gap. At δ = 0, the system is insulating and thus z = 0. With increasing δ, z increases according to the decrease of the electron correlation effect. We find that z for the La system is smaller than that for the Hg system, indicating that the electron correlation effect is stronger in the La system. This can be attributed to the larger ∆ dp = ε d x 2 −y 2 − ε p for the La system, which results in the larger d-orbital occupancy of holes when it is doped. The electron correlation has a dual effect: One is to enhance the superconducting pairing and the other is to reduce the mobility of Cooper pairs. In the present case, the latter effect would be dominant and thus P dd for the La system is suppressed compared with the Hg system. This tendency is consistent with the negative correlation between ∆ dp and T c as mentioned in Sec. I.

C. Effect of V dd
Now we discuss the effect of the Coulomb interaction between nearest-neighbor d orbitals, V dd . V dd is expected to be smaller than other Coulomb interactions, U d , U ′ d , U p , V dp , and V pp [65]. However, V dd directly affects the charge and spin correlations between nearest-neighbor d electrons, which dominate the properties of cuprate superconductors. As discussed in this section, we verify that the superconductivity in the model studied here is more sensitive to the value of V dd than other Coulomb interactions. Hereafter, we treat only U d and V dd as independent parameters, and set other Coulomb interactions as (U ′ d , J, U p , V dp , V pp ) = (0.8, 0.1, 0.5, 0.25, 0.2) U d with reference to Ref. [65] Figure 8(a) shows the ground state phase diagram for the Hg system at δ = 0 within the paramagnetic phase. When V dd = 0, a metal-insulator transition occurs at U d /t 1 ∼ 6.3, assuming the paramagnetic phase. The metal-insulator transition is detected by monitoring the jump in the total momentum distribution function n tot (k) as well as the long-range behavior of the charge Jastrow factor P Jc [72]. This transition is a Mott metalinsulator transition because δ = 0 corresponds to the case with one hole per unit cell, n hole = 1, where the system cannot be a band insulator. With the introduction of small but finite V dd (one order of magnitude smaller than U d ), the metallic region is substantially enlarged. This is understood because the densities of empty sites and doubly-occupied sites increase to reduce the energy loss of V dd , thus effectively weakening U d by V dd , and then the insulating phase is destabilized.
It is noteworthy that the AF insulator is always more stable than the paramagnetic insulator in the present parameter space. With increasing U d , the system undergoes the phase transition from the metallic phase to a Slatertype AF insulator [a blue line in Fig. 8(b)], followed by the crossover to a Mott-type AF insulator [a red dotted line in Fig. 8(b)]. Here, a Slater-type AF insulator is an insulator that becomes metallic without AF order, while a Mott-type AF insulator is an insulator that remains insulating without AF order. The crossover line in Fig. 8(b) is hence identical with the paramagnetic metal-insulator transition line in Fig. 8(a). As explained next, it is crucially important for the appearance of high-T c superconductivity whether the AF insulator in the parent compounds (δ=0) is Slater-type or Mott-type [31,32,72,73].
Next we study the effect of V dd on superconductivity. Figure 9 shows the superconducting correlation function P dd at U d /t 1 = 8 and V dd /t 1 = 0, 0.3, and 0.6 for the Hg system. As V dd increases, P dd is suppressed and the peak position is moved to a smaller value of δ. In particular, a significant suppression is observed when V dd /t 1 = 0.6. In this case, P dd vs. δ does not show a dome-shaped behavior but instead an almost monotonic decrease, which rather reminds us of electron-doped cuprates [74,75]. This observation of P dd vs. δ corresponds to the fact that the system is out of the Mott insulator region at δ = 0 (see Fig. 8). The view of a "doped Mott insulator" is thus no longer valid. Similar claims have been made in the study of one-band Hubbard models [76,77]. Our result suggests that it is essential to start from the Mott insulator region at δ = 0 to reproduce the domeshaped behavior observed experimentally in hole-doped cuprates. This can be considered as a good criterion for choosing the reasonable Coulomb interaction parameters of an effective model for cuprates.

D. Phase competition
Finally, we briefly discuss the competition between superconductivity and other phases. The energy comparison among various phases in the ground state is a subtle problem and depends significantly on the numerical methods. The VMC method used here tends to overestimate the magnetic long-range ordered phases, although it is much improved as compared with a mean-field type approximation. In fact, we find that the AF phase and the stripe phase with both spin and charge modulations have lower variational energies than d x 2 −y 2 -wave superconductivity for δ < 0.3. Nevertheless, we believe that the present results capture the essence of the material dependence of cuprate superconductors and the conclusion is unchanged, when the improved wave functions incorporating the quantum fluctuations suppress these overestimated competing orders. We also note that the effect of V dd is also important for the phase competition. This is because most of the competing phases including the phase separation are governed by the correlation between nearest-neighbor d electrons. This is also the case in the one-band Hubbard model [78,79]. The detailed ground state phase diagram, including various competing phases, for the four-band d-p model is left for a future study.

IV. DISCUSSION AND SUMMARY
To obtain the unified description of cuprate superconductors, we have studied the four-band d-p model for the La 2 CuO 4 and HgBa 2 CuO 4 systems. We have shown that a lower ε d z 2 with respect to ε d x 2 −y 2 and a smaller ∆ dp (> 0) lead to a higher T c . The former results in a more localized d z 2 -orbital based band that do not interfere the superconductivity. The latter results in a larger z, namely, a weaker electron correlation effect, which promotes the itinerancy of mobile carriers and thus enhances superconductivity. The present four-band d-p model covers these two factors, beyond the usual one-band and even three-band d-p models. Therefore, this model is considered to be a minimal model that can properly describe the material dependence of cuprate superconductors, and thus it can also provide a valuable guideline to design new materials with a higher T c .
The effect of V dd has also been investigated. Al-though the value of V dd is small compared with other Coulomb interaction parameters, it substantially affects the ground state property of the system. V dd weakens the effective U d and induces the paramagnetic metalinsulator transition, or the crossover from a Slater insulator to a Mott insulator. The stability of superconductivity is also affected by V dd . Considering the doping dependence, we have to start from the Mott insulator region at δ = 0 to obtain the stable superconductivity and the dome-shaped dependence of P dd as a function of δ. Therefore, the appropriate estimation of V dd is important for the modelling of cuprates, as in other stronglycorrelated electron systems where various phases compete.

ACKNOWLEDGMENTS
The authors thank K. Kuroki for useful discussions. The computation has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and the supercomputer system HOKUSAI in RIKEN. This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas "Quantum Liquid Crystals" (KAK-ENHI Grant No. JP19H05825) from JSPS of Japan, and also supported by JSPS KAKENHI (Grant Nos. JP21H04446, JP20K03847, JP19K23433, JP19H01842, and JP18H01183).

Appendix: Construction of the trial wave function
Here, we describe the details of the trial wave function together with the noninteracting tight-binding model obtained on the basis of the first-principles calculations in Sec. III A. The construction of the trial wave function is the most important part for the VMC method. Depending on the trial state, both real and k space representations are used.

Noninteracting energy band
First, we describe how to construct the noninteracting tight-binding energy band discussed in Sec. III A. The noninteracting energy band is obtained by diagonalizing the following one-body Hamiltonian: with the hopping matrix elements given as where c † kασ (c kασ ) is a creation (annihilation) operator of an electron with momentum k, spin σ (=↑, ↓), and orbital α (= 1, 2, 3, 4) corresponding to (d x 2 −y 2 , d z 2 , p x , p y ), respectively. The bonds between nearest-neighbor Cu sites are set as unit vectors (|e x | = |e y | = 1) and the bonds between nearest-neighbor Cu and O sites are one-half of them (see Fig. 1). The hopping integrals t i (i = 1 − 6) and the site energy of each orbital ε α are determined by fitting the band structures that are obtained by the LDA or QSGW calculation. The specific values for the La, La(revised), and Hg systems are listed in Table I. Equation (A.2) is obtained by diagonalizing the Hamiltonian matrix in Eq. (A.1) and is the same with Eq. (3) in Sec. II A. E m (k) is the noninteracting energy band characterized by the wave vector k and the energy band index m (= 1, 2, 3, 4) with a † kmσ (a kmσ ) being a creation (annihilation) operator of the corresponding energy band with spin σ.

Trial wave function a. Superconductivity
To construct the trial wave function for superconductivity, we employ the Bogoliubov de-Gennes (BdG) type Hamiltonian in real space [80], i.e., (A.13) Here, T αβ ijσ is obtained from the Fourier transform of the matrix in Eq. (A.1) with renormalized hopping integrals t i and also includes the chemical potential term. The chemical potential µ is set to the Fermi energy in the non-interacting limit. ∆ αβ ij corresponds to an anomalous part that describes the superconducting pairing in real space. Therefore, the variational parameters to be optimized in |Φ aret i (i = 2 − 6) and {∆ αβ ij } witht 1 = t 1 being fixed as a unit of energy. In this study, we consider the pairing between nearest-neighbor orbitals, d-d, d-p x , d-p y , p x -p x , and p y -p y , where d denotes the d x 2 −y 2 orbital. In the paramagnetic phase, we simply set ∆ αβ ij = 0. We can also construct the trial wave function with the band (k space) representation, where ∆ mn (k) ∝ a km↑ a −kn↓ is the variational parameter. However, we find that the trial wave function with the real space representation always gives lower (i.e., better) energy than that with the band representation, especially, for large Coulomb interaction parameters. This is because the Coulomb interaction depends on the orbital, not on the band, and thus the trial wave function with the real space (orbital) representation gives the better result.

b. Uniform spin AF and stripe phases
As mentioned in Sec. II B, various long-range orderings of charge and spin can be described by introducing {ρ α i } and {m α i }. A uniform spin AF phase with A and B sublattices for the orbital α is expressed with the staggered potential m α i = −s σ m α for A sublattice +s σ m α for B sublattice (A.14) where s σ = 1(−1) for σ =↑ (↓). For a stripe phase with charge and spin periodicities λ α c = 2π/q α c and λ α s = 2π/q α s , respectively, the following potentials with spatial modulation in the x direction should be introduced: where x α c and x α s control the relative phases of charge and spin orderings, respectively. For example, the stripe phase observed around δ = 1/8 in several cuprate superconductors corresponds to λ c = 4 and λ s = 8, although the orbital dependence and relative phase are still under debate.