Unified description of cuprate superconductors using 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_{\text{kin}}$ is described by

where Eq. (2) is the kinetic term in an orbital representation and Eq. (3) is in a band representation. $c^{\dagger}_{i\alpha\sigma}$ ($c_{i\alpha\sigma}$) is a creation (annihilation) operator of an electron at site $i$ with spin $\sigma\,(=\uparrow,\downarrow)$ and orbital $\alpha\,(=1,2,3,4)$ corresponding to ($d_{x^{2}-y^{2}}$, $d_{z^{2}}$, $p_{x}$, $p_{y}$), respectively. $t^{\alpha\beta}_{ij}$ denotes a hopping integral between orbital $\alpha$ at site $i$ and orbital $\beta$ at site $j$. $t^{\alpha\alpha}_{ii}$ is a site energy $\varepsilon_{\alpha}$ for orbital $\alpha$ at site $i$. Equation (3) is obtained by diagonalizing Eq. (2), and the energy eigenvalue $E_{m}(\textbf{k})$ is characterized by the wave vector k and the energy band index $m\,(=1,2,3,4)$. $a^{\dagger}_{\textbf{k}m\sigma}$ ($a_{\textbf{k}m\sigma}$) is a creation (annihilation) operator of the corresponding energy band with spin $\sigma$. 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 $\delta$ of holes per unit cell that are introduced into the system at half filling.

The Coulomb interaction term $H_{\text{int}}$ is composed of eight terms,

Here, $n^{\alpha}_{i}=n^{\alpha}_{i\uparrow}+n^{\alpha}_{i\downarrow}$ with $n^{\alpha}_{i\sigma}=c^{\dagger}_{i\alpha\sigma}c_{i\alpha\sigma}$ is the number operator and $\textbf{S}^{\alpha}_{i}$ is the spin angular momentum operator at site $i$ with orbital $\alpha$. $d_{1}$ and $d_{2}$ are abbreviations for $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals, respectively, and $n^{d}_{i}=n^{d_{1}}_{i}+n^{d_{2}}_{i}$. $U_{d},U_{d}^{\prime},J,$ and $J^{\prime}$ represent on-site intraorbital, interorbital, Hund’s coupling, and pair-hopping interactions between $d$ orbitals, respectively. In this study, we set $J^{\prime}=J$ and $U_{d}=U^{\prime}_{d}+2J$ [48]. The on-site intraorbital Coulomb interaction within $p$ orbitals, $U_{p}$, is also introduced. The last three terms in $H_{\text{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 $\sum_{\left<i,j\right>}$ runs over all pairs of nearest-neighbor orbitals located at site $i$ and $j$.

In addition, the following double counting correction term $H_{\text{dc}}$ is introduced,

where $\langle n^{d}\rangle_{0}=\frac{1}{N_{\text{S}}}\sum_{i}\langle n^{d}_{i}\rangle$ and $\left<n^{p}\right>_{0}=\frac{1}{N_{\text{S}}}\sum_{i}\langle(n^{p_{x}}_{i}+n^{p_{y}}_{i})\rangle$ are the average electron density of the $d$ and $p$ orbitals in the noninteracting limit and $N_{\text{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_{\text{dc}}$ from the Hamiltonian to correct the energy shift. This is one of the reasonable treatments to avoid the double counting.

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,

$\left|\Phi\right>$ is a one-body part constructed by diagonalizing the one-body Hamiltonian including the off-diagonal elements $\{\rho^{\alpha}_{i}\}$, $\{m^{\alpha}_{i}\}$, and $\{\Delta^{\alpha\beta}_{ij}\}$, which induce long-range ordering of charge, spin, and superconductivity, respectively. The renormalized hopping integrals $\{\tilde{t}^{\alpha\beta}_{ij}\}$ are also included in $\left|\Phi\right>$ as variational parameters. The explicit forms of them are described in Appendix.

The Gutzwiller factor

is extended to two-orbital systems [53, 54]. In $P^{(2)}_{\text{G}}$, possible 16 patterns of charge and spin configuration of the $d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ orbitals at each site $\left|\gamma\right>$, i.e., $\left|0\right>=\left|0\;0\right>$, $\left|1\right>=\left|0\uparrow\right>$, $\cdots$, $\left|15\right>=\left|\uparrow\downarrow\;\uparrow\downarrow\right>$, are differently weighted with $\text{e}^{-g_{\gamma}}$ and $\{g_{\gamma}\}$ are optimized as variational parameters.

The remaining operators

and

are charge and spin Jastrow factors, which control long-range charge and spin correlations, respectively. $s^{z}_{i\alpha}$ is the $z$ component of the spin angular momentum operator at site $i$ with orbital $\alpha$. We set $v^{\text{c}}_{ii\alpha\beta}=v^{\text{s}}_{ii\alpha\beta}=0$ for $\alpha,\beta=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^{(2)}_{\text{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 $\left|\Psi\right>$ are therefore $\{\tilde{t}^{\alpha\beta}_{ij}\}$, $\{\Delta^{\alpha\beta}_{ij}\}$, $\{g_{\gamma}\}$, $\{v^{\text{c}}_{ij\alpha\beta}\}$, and $\{v^{\text{s}}_{ij\alpha\beta}\}$. They are simultaneously optimized using stochastic reconfiguration method [55]. We show results for $N_{\text{S}}$=24$\times$24=576 unit cells (and thus 576$\times$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.

First, we discuss the material dependence of the band structure. As a typical example of single-layer hole-doped cuprates, we study the La₂CuO₄ and HgBa₂CuO₄ 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 $\varepsilon_{\alpha}$. The parameter sets determined for these systems are listed in Table 1 and the explicit form of the tight-binding model is described in Appendix A.1. Note that the estimated site energy of the Cu $d_{z^{2}}$ orbital, $\varepsilon_{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 $\varepsilon_{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 $\varepsilon_{d_{z^{2}}}$ is lower and $t_{4}$ is also slightly smaller than those estimated on the basis of the LDA calculation (See Table 1).

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, $\Delta_{dp}=\varepsilon_{d_{x^{2}-y^{2}}}-\varepsilon_{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 $\Delta_{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^{\prime}_{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.3, we find that even small $V_{dd}$ can substantially affect the property of the system.

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^{\prime}_{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^{\prime}_{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 $\delta=0$ within the paramagnetic phase. When $V_{dd}=0$, a metal-insulator transition occurs at $U_{d}/t_{1}\sim 6.3$, assuming the paramagnetic phase. The metal-insulator transition is detected by monitoring the jump in the total momentum distribution function $n^{\text{tot}}(\textbf{k})$ as well as the long-range behavior of the charge Jastrow factor $P_{\text{J}_{\text{c}}}$ [72]. This transition is a Mott metal-insulator transition because $\delta=0$ corresponds to the case with one hole per unit cell, $n_{\text{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 Slater-type 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_{\text{c}}$ superconductivity whether the AF insulator in the parent compounds ($\delta$=0) is Slater-type or Mott-type [31, 32, 73, 72].

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 $\delta$. In particular, a significant suppression is observed when $V_{dd}/t_{1}=0.6$. In this case, $P^{dd}$ vs. $\delta$ 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. $\delta$ corresponds to the fact that the system is out of the Mott insulator region at $\delta=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 $\delta=0$ to reproduce the dome-shaped 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.

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 $\delta<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.