Coexistence of $s$- and $d$-wave gaps due to pair-hopping and exchange interactions

The $t$–$J$ model is one of the model Hamiltonians that form the basis of many theoretical studies of strongly correlated electron systems [1, 2]. The $t$–$J$ model can also be derived as the low-energy effective Hamiltonian of the two-dimensional (2D) multiband Hubbard model [3, 4, 5], regarded as the fundamental model Hamiltonian for the high-$T_{\mathrm{c}}$ superconducting cuprate (HTSC). Many theoretical studies of HTSC to date use the $t$–$J$ model as the model Hamiltonian [6, 7, 8, 10, 9, 11, 12, 13]. These studies often exclude the double occupancy of Cu sites by $d$ electrons, considering that the on-site Coulomb repulsion between $d$ orbitals is much larger than the transfer energy between the $d$ and $p$ orbitals. As a result, the $t$–$J$ model contains only one electron (or hole) band and a localized spin.

However, the double occupancy of Cu sites need not necessarily be excluded when the on-site Coulomb repulsion $U$ is comparable to the transfer energy. Relaxing the single occupancy constraint and explicitly considering $U$ instead results in the $t$–$J$–$U$ model that includes both the $t$–$J$ model and the single-band Hubbard model as one of its limits [14, 15, 16, 17, 18, 19, 20]. Thus, the $t$–$J$–$U$ model serves as an interpolation between the $t$–$J$ model and the single-band Hubbard model and is able to account for more properties caused by strong correlation. However, the charge transfer gap should be comparable to $t$ in the charge transfer regime. In this case, $p$ electron scattering by $d$ electrons cannot be negligible, and both $p$ and $d$ electrons must be considered.

In this paper, I derive the three-band $t$–$J$–$U$ model from the 2D three-band Hubbard model as its effective Hamiltonian by using the Schrieffer–Wolff (SW) transformation [21] and assume that double occupancy is not excluded. In my model, the pair-hopping interaction between the $d$ and $p$ bands exists separately from the exchange interaction. Treating these interactions using iterative perturbation theory (IPT) approximation, I investigate the superconductivity of the model in a strong coupling framework. The results show that the multicomponent superconductivity emerges with the hole doping, which introduces the $d$-$p$ band hybridization through exchange and pair-hopping interactions. This emergence of the superconductivity is due to the pair-hopping and exchange interactions via the Suhl-Kondo (SK) mechanism [22, 23, 24], which stabilizes the superconducting gaps with different signs in a multiband system. In the superconducting state, the extended $s$- and $d_{x^{2}-y^{2}}$-wave superconducting gaps coexist, and the $s$- and $d$-wave gaps emerge due to the pair-hopping and exchange interactions, respectively.

Consider the three-band Hubbard model [25] that expresses the Hamiltonian as ${\mathcal{H}}={\mathcal{H}}_{0}+\sum_{\alpha}{\mathcal{H}}_{1}^{\alpha}$, where

$${\mathcal{H}}_{0}=\varepsilon_{d}\sum_{j\sigma}d_{j\sigma}^{\dagger}d_{j\sigma}+\varepsilon_{p}\sum_{\alpha}\sum_{{\mathbf{k}}\sigma}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}\sigma}^{\alpha}+U\sum_{j}d_{j\uparrow}^{\dagger}d_{j\uparrow}d_{j\downarrow}^{\dagger}d_{j\downarrow}$$

(1)

and

$${\mathcal{H}}_{1}^{\alpha}=\frac{1}{\sqrt{N}}\sum_{j}\sum_{{\mathbf{k}}\sigma}\left(V_{\alpha{\mathbf{k}}}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}d_{j\sigma}+{\mathrm{H.c.}}\right).$$

(2)

Here, $\alpha\in\{x,y\}$; $d_{j\sigma}(d_{j\sigma}^{\dagger})$ is the annihilation (creation) operator for the $d$ electron of spin $\sigma$ at Cu site $j$; $p_{{\mathbf{k}}\sigma}^{\alpha}(p_{{\mathbf{k}}\sigma}^{\alpha\dagger})$ is the annihilation (creation) operator for $p^{\alpha}$ electrons of spin $\sigma$ with momentum ${\mathbf{k}}$, based on oxygen sites in real space; $\varepsilon_{d}$ and $\varepsilon_{p}$ are the $d$ and $p$ electron site energies, respectively; $U$ is the on-site Coulomb repulsion between $d$ orbitals; and $N$ is the number of k-space points in the first Brillouin zone (FBZ). The lattice constant of the square lattice of Cu sites is the length unit. Thus, $V_{x{\mathbf{k}}}=2{\mathrm{i}}t_{pd}\sin\frac{k_{x}}{2}$ and $V_{y{\mathbf{k}}}=-2{\mathrm{i}}t_{pd}\sin\frac{k_{y}}{2}$, where $t_{pd}$ is the transfer energy between the $d$ orbital and the neighbouring $p^{\alpha}$ orbital.

In order to derive the effective Hamiltonian for $\mathcal{H}$, I adopt the SW transformation as follows:

	$\displaystyle e^{\sum_{\alpha}{\mathcal{S}}^{\alpha}}{\mathcal{H}}e^{-\sum_{\beta}{\mathcal{S}}^{\beta}}$	$\displaystyle=$	$\displaystyle{\mathcal{H}}_{0}+\sum_{\alpha}{\mathcal{H}}_{1}^{\alpha}+\sum_{\alpha}\left[{\mathcal{S}}^{\alpha},{\mathcal{H}}_{0}\right]+\sum_{\alpha\beta}\left[{\mathcal{S}}^{\alpha},{\mathcal{H}}_{1}^{\beta}\right]+\frac{1}{2}\sum_{\alpha\beta}\left[{\mathcal{S}}^{\alpha},\left[{\mathcal{S}}^{\beta},{\mathcal{H}}_{0}\right]\right]+\ldots$		(3)
		$\displaystyle=$	$\displaystyle{\mathcal{H}}_{0}+\frac{1}{2}\sum_{\alpha\beta}\left[{\mathcal{S}}^{\alpha},{\mathcal{H}}_{1}^{\beta}\right]+\ldots,$

using ${\mathcal{H}}_{1}^{\alpha}+\left[{\mathcal{S}}^{\alpha},{\mathcal{H}}_{0}\right]=0$ and

$${\mathcal{S}}^{\alpha}=\frac{1}{\sqrt{N}}\sum_{j}\sum_{{\mathbf{k}}\sigma}\left(\frac{V_{\alpha{\mathbf{k}}}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}}{\mathit{\Delta}_{pd}-U}\,n_{d\,j-\sigma}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}d_{j\sigma}+\frac{V_{\alpha{\mathbf{k}}}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}}{\mathit{\Delta}_{pd}}(1-n_{d\,j-\sigma})p_{{\mathbf{k}}\sigma}^{\alpha\dagger}d_{j\sigma}\right)-{\mathrm{H.c.}}$$

(4)

Here, ${\mathit{\Delta}}_{pd}\equiv\varepsilon_{p}-\varepsilon_{d}$, $n_{d\,j\sigma}\equiv d_{j\sigma}^{\dagger}d_{j\sigma}$, and H.c. indicates the Hermitian conjugate of the terms already written. The observable $n_{d\,j\sigma}$ has $0$ or $1$ as its eigenvalue for each $j$ and $\sigma$. Using Eqs. (2) and (4), the following results:

	$\displaystyle\left[{\mathcal{S}}^{\alpha},{\mathcal{H}}_{1}^{\beta}\right]$		$\displaystyle=-\frac{\delta_{\alpha\beta}}{N}\sum_{jj^{\prime}}\sum_{{\mathbf{k}}\sigma}\left(\frac{n_{d\,j-\sigma}}{\mathit{\Delta}_{pd}-U}+\frac{1-n_{d\,j-\sigma}}{\mathit{\Delta}_{pd}}\right)V_{\alpha{\mathbf{k}}}V_{\beta{\mathbf{k}}}^{*}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}e^{{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j^{\prime}}}d_{j^{\prime}\sigma}^{\dagger}d_{j\sigma}$
			$\displaystyle+\frac{1}{N}\sum_{j}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\sigma}\left(\frac{n_{d\,j-\sigma}}{\mathit{\Delta}_{pd}-U}+\frac{1-n_{d\,j-\sigma}}{\mathit{\Delta}_{pd}}\right)V_{\alpha{\mathbf{k}}}V_{\beta{\mathbf{k}}^{\prime}}^{*}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}e^{{\mathrm{i}}{\mathbf{k}}^{\prime}\cdot{\mathbf{R}}_{j}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}^{\prime}\sigma}^{\beta}$
			$\displaystyle-\frac{1}{N}\sum_{j}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\sigma}\left(\frac{1}{\mathit{\Delta}_{pd}-U}-\frac{1}{\mathit{\Delta}_{pd}}\right)V_{\alpha{\mathbf{k}}}V_{\beta{\mathbf{k}}^{\prime}}^{*}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}e^{{\mathrm{i}}{\mathbf{k}}^{\prime}\cdot{\mathbf{R}}_{j}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}^{\prime}-\sigma}^{\beta}d_{j\,-\sigma}^{\dagger}d_{j\sigma}$
			$\displaystyle-\frac{1}{N}\sum_{j}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\sigma}\left(\frac{1}{\mathit{\Delta}_{pd}-U}-\frac{1}{\mathit{\Delta}_{pd}}\right)V_{\alpha{\mathbf{k}}}V_{\beta{\mathbf{k}}^{\prime}}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}e^{-{\mathrm{i}}{\mathbf{k}}^{\prime}\cdot{\mathbf{R}}_{j}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}^{\prime}-\sigma}^{\beta\dagger}d_{j\,-\sigma}d_{j\sigma}+{\mathrm{H.c.}}$

Hereafter, I consider only the first two terms of the right-hand side of Eq. (3), i.e., up to the second order of $t_{pd}$. Now, I assume that the distribution of the $d$ electron is spatially uniform in the ground state and that the ground state is paramagnetic. Thus, $\langle n_{d\,j\uparrow}+n_{d\,j\downarrow}\rangle_{0}\equiv n_{d}$ and $\langle n_{d\,j\uparrow}\rangle_{0}=\langle n_{d\,j\downarrow}\rangle_{0}$ for any $j$ where $n_{d}$ is a c-number equal to the number of $d$ electrons in the ground state, where $\langle\cdots\rangle_{0}$ indicates the average in the ground state. I apply this approximation to Eqs. (3) and (LABEL:eq:06) and treat $n_{d}$ as a parameter that should be determined self-consistently. When I set $\varepsilon_{p}$ to zero, i.e., $\mathit{\Delta}_{pd}=-\varepsilon_{d}$, and omit the constant terms, I obtain the effective Hamiltonian:

$${\mathcal{H}}_{\mathrm{eff}}={\mathcal{H}}_{\mathrm{HF}}+{\mathcal{H}}_{\mathrm{ex}}+{\mathcal{H}}_{\mathrm{pair}}+{\mathcal{H}}_{U}^{\prime}.$$

(6)

${\mathcal{H}}_{\mathrm{HF}}$ is the Hartree-Fock approximation of ${\mathcal{H}}_{0}$:

$\displaystyle{\mathcal{H}}_{\mathrm{HF}}$

$\displaystyle=$

$\displaystyle\sum_{{\mathbf{k}}\sigma}\varepsilon_{d{\mathbf{k}}}d_{{\mathbf{k}}\sigma}^{\dagger}d_{{\mathbf{k}}\sigma}+\sum_{\alpha\beta}\sum_{{\mathbf{k}}\sigma}\varepsilon_{\alpha\beta{\mathbf{k}}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}\sigma}^{\beta},$

(7)

where $d_{{\mathbf{k}}\sigma}^{\dagger}=\frac{1}{\sqrt{N}}\sum_{j}d_{j\sigma}^{\dagger}e^{{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}},d_{{\mathbf{k}}\sigma}=\frac{1}{\sqrt{N}}\sum_{j}d_{j\sigma}e^{-{\mathrm{i}}{\mathbf{k}}\cdot{\mathbf{R}}_{j}}$,

$$\varepsilon_{d{\mathbf{k}}}=\varepsilon_{d}+\frac{U}{2}n_{d}+t\left(v_{x{\mathbf{k}}}v_{x{\mathbf{k}}}^{*}+v_{y{\mathbf{k}}}v_{y{\mathbf{k}}}^{*}\right),$$

(8)

and

$$\varepsilon_{\alpha\beta{\mathbf{k}}}=\left(Jn_{d}-t\right)v_{\alpha{\mathbf{k}}}v_{\beta{\mathbf{k}}}^{*},$$

(9)

with $v_{x{\mathbf{k}}}={\mathrm{i}}\sin\frac{k_{x}}{2}$, $v_{y{\mathbf{k}}}=-{\mathrm{i}}\sin\frac{k_{y}}{2}$,

$$t=4\,t_{pd}^{2}\left(\frac{n_{d}}{\varepsilon_{d}+U}+\frac{1-n_{d}}{\varepsilon_{d}}\right),$$

(10)

and

$$J=2\,t_{pd}^{2}\left(\frac{1}{\varepsilon_{d}+U}-\frac{1}{\varepsilon_{d}}\right).$$

(11)

${\mathcal{H}}_{\mathrm{ex}}$ is an exchange interaction term:

$${\mathcal{H}}_{\mathrm{ex}}=\frac{J}{N}\sum_{\alpha\beta}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\sigma}\sum_{{\mathbf{q}}}v_{\alpha{\mathbf{k}}}v_{\beta{\mathbf{k}}^{\prime}}^{*}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}^{\prime}-\sigma}^{\beta}d_{{\mathbf{k}}^{\prime}+{\mathbf{q}}\,-\sigma}^{\dagger}d_{{\mathbf{k}}+{\mathbf{q}}\,\sigma}+{\mathrm{H.c.}}$$

(12)

${\mathcal{H}}_{\mathrm{pair}}$ is a pair-hopping term:

$${\mathcal{H}}_{\mathrm{pair}}=\frac{J}{N}\sum_{\alpha\beta}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\sigma}\sum_{{\mathbf{q}}}v_{\alpha{\mathbf{k}}}v_{\beta{\mathbf{k}}^{\prime}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}p_{{\mathbf{k}}^{\prime}-\sigma}^{\beta\dagger}d_{{\mathbf{k}}^{\prime}-{\mathbf{q}}\,-\sigma}d_{{\mathbf{k}}+{\mathbf{q}}\,\sigma}+{\mathrm{H.c.}}$$

(13)

${\mathcal{H}}_{U}^{\prime}$ is the Coulomb interaction term excluding the component with ${\mathbf{q}}={\mathbf{0}}$:

$${\mathcal{H}}_{U}^{\prime}=\frac{U}{N}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}}\sum_{{\mathbf{q}}\neq{\mathbf{0}}}d_{{\mathbf{k}}+{\mathbf{q}}\,\uparrow}^{\dagger}d_{{\mathbf{k}}\uparrow}d_{{\mathbf{k}}^{\prime}-{\mathbf{q}}\,\downarrow}^{\dagger}d_{{\mathbf{k}}^{\prime}\downarrow}.$$

(14)

As a consequence, ${\mathcal{H}}_{\mathrm{eff}}$ [Eq. (6)] can be characterized by the three parameters $t$ [Eq. (10)], $J$ [Eq. (11)], and $U$, and it can be regarded as the three-band $t$–$J$–$U$ model.

Here, $t$ in Eq. (10) is positive near the half-filling in the charge-transfer regime, i.e., $U>-\varepsilon_{d}>0$. For instance, in the case $\varepsilon_{d}=-U/2$, $t>0$ for $n_{d}>0.5$, and the $d$ electron band dispersion $\varepsilon_{d{\mathbf{k}}}$ in Eq. (8) is the same as that for the single-band Hubbard model on a square lattice. $J$ in Eq. (11) is always positive in the charge-transfer regime. Thus, ${\mathcal{H}}_{\mathrm{ex}}$ in Eq. (12) describes the transverse component of the antiferromagnetic exchange interaction between the $d$ and $p$ electrons, while the longitudinal component of this interaction narrows the bandwidth of $\varepsilon_{\alpha\beta{\mathbf{k}}}$ in Eq. (9) from $t$ to $t-Jn_{d}$. Further, ${\mathcal{H}}_{\mathrm{ex}}$ indicates that the $p$ electron is affected by the spin fluctuation of the $d$ electron. As will be shown later, the $d$-wave superconducting gap composed of $d$ and $p$ electrons emerges from ${\mathcal{H}}_{\mathrm{ex}}$. ${\mathcal{H}}_{\mathrm{pair}}$ in Eq. (13) appears for the first time by considering the double occupancy of Cu sites. The pair-hopping term is not included in the single-band $t$–$J$ model if double occupancy is excluded. In the model that includes the pair-hopping interaction, electrons favour pair formation [26]. This is also true in the presence of the on-site interaction [27, 28] and in the zero-bandwidth limit [29]. Thus, the pair-hopping term in my model is expected to provide superconductivity in another way.

I introduce another assumption according to the speculation about the ground state of the three-band Hubbard model [30]. In the normal ground state, the $d$ and $p$ electrons should be combined to construct coherent quasi-particles through hybridization. The matrix elements of the hybridization between the $d$ and $p$ electrons can be found in the components with ${\mathbf{q}}={\mathbf{0}}$ in Eqs. (12) and (13) as follows. Defining

	$\displaystyle h_{pd}=-\frac{{\mathrm{i}}}{N}\sum_{\alpha}\sum_{{\mathbf{k}}}\left[v_{\alpha{\mathbf{k}}}\langle p_{{\mathbf{k}}\uparrow}^{\alpha\dagger}d_{{\mathbf{k}}\uparrow}\rangle_{0}-v_{\alpha{\mathbf{k}}}^{*}\langle d_{{\mathbf{k}}\uparrow}^{\dagger}p_{{\mathbf{k}}\uparrow}^{\alpha}\rangle_{0}\right]=-\frac{{\mathrm{i}}}{N}\sum_{\alpha}\sum_{{\mathbf{k}}}\left[v_{\alpha{\mathbf{k}}}\langle p_{{\mathbf{k}}\downarrow}^{\alpha\dagger}d_{{\mathbf{k}}\downarrow}\rangle_{0}-v_{\alpha{\mathbf{k}}}^{*}\langle d_{{\mathbf{k}}\downarrow}^{\dagger}p_{{\mathbf{k}}\downarrow}^{\alpha}\rangle_{0}\right],$
			(15)

${\mathcal{H}}_{\mathrm{eff}}$ can be rewritten as

$${\mathcal{H}}_{\mathrm{eff}}={\mathcal{H}}_{0}^{\prime}+{\mathcal{H}}_{\mathrm{ex}}^{\prime}+{\mathcal{H}}_{\mathrm{pair}}^{\prime}+{\mathcal{H}}_{U}^{\prime},$$

(16)

where

$${\mathcal{H}}_{0}^{\prime}={\mathcal{H}}_{\mathrm{HF}}+{\mathrm{i}}Jh_{pd}\sum_{\alpha}\sum_{{\mathbf{k}}\sigma}\left(v_{\alpha{\mathbf{k}}}p_{{\mathbf{k}}\sigma}^{\alpha\dagger}d_{{\mathbf{k}}\sigma}-v_{\alpha{\mathbf{k}}}^{*}d_{{\mathbf{k}}\sigma}^{\dagger}p_{{\mathbf{k}}\sigma}^{\alpha}\right).$$

(17)

Here, $\langle A\rangle_{0}$ in Eq. (15) means the expectation value of $A$ in the ground state of ${\mathcal{H}}_{0}^{\prime}$. ${\mathcal{H}}_{\mathrm{ex}}^{\prime}$ and ${\mathcal{H}}_{\mathrm{pair}}^{\prime}$ indicate the exchange interaction and pair-hopping terms excluding the component with ${\mathbf{q}}={\mathbf{0}}$ from Eqs. (12) and (13), respectively. Thus, in the ground state of ${\mathcal{H}}_{0}^{\prime}$, the $d$ and $p$ electrons are combined to construct the coherent quasi-particles when $h_{pd}>0$.

Hereafter, I treat ${\mathcal{H}}_{0}^{\prime}$ as the unperturbed part of ${\mathcal{H}}_{\mathrm{eff}}$ on the assumption that $h_{pd}>0$. I diagonalize ${\mathcal{H}}_{0}^{\prime}$ and derive the unperturbed Green functions as follows:

	$\displaystyle G_{dd}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})=\frac{{\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{xx{\mathbf{k}}}-\varepsilon_{yy{\mathbf{k}}}}{({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{+})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{-})},$		(18)
	$\displaystyle G_{d\alpha}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})=\frac{{\mathrm{i}}Jh_{pd}v_{\alpha\mathbf{k}}}{({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{+})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{-})},$		(19)
	$\displaystyle G_{\alpha d}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})=\frac{-{\mathrm{i}}Jh_{pd}v_{\alpha\mathbf{k}}^{*}}{({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{+})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{-})},$		(20)

and

	$\displaystyle\left(\begin{array}[]{cc}G_{xx}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})&G_{xy}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})\\ G_{yx}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})&G_{yy}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n})\\ \end{array}\right)=\frac{1}{({\mathrm{i}}\epsilon_{n}+\mu)({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{+})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{{\mathbf{k}}}^{-})}$		(23)
	$\displaystyle\times\left(\begin{array}[]{cc}({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{d{\mathbf{k}}})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{yy{\mathbf{k}}})-J^{2}h_{pd}^{2}v_{y{\mathbf{k}}}v_{y{\mathbf{k}}}^{*}&({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{d{\mathbf{k}}})\varepsilon_{yx{\mathbf{k}}}+J^{2}h_{pd}^{2}v_{y{\mathbf{k}}}v_{x{\mathbf{k}}}^{*}\\ ({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{d{\mathbf{k}}})\varepsilon_{xy{\mathbf{k}}}+J^{2}h_{pd}^{2}v_{x{\mathbf{k}}}v_{y{\mathbf{k}}}^{*}&({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{d{\mathbf{k}}})({\mathrm{i}}\epsilon_{n}+\mu-\varepsilon_{xx{\mathbf{k}}})-J^{2}h_{pd}^{2}v_{x{\mathbf{k}}}v_{x{\mathbf{k}}}^{*}\\ \end{array}\right).$		(26)

Here, I use the fermion Matsubara frequencies, $\epsilon_{n}=\pi T(2n+1)$, with integer $n$ and temperature $T$. $\mu$ is the chemical potential and

$$\varepsilon_{{\mathbf{k}}}^{\pm}=\frac{\varepsilon_{d{\mathbf{k}}}+\varepsilon_{xx{\mathbf{k}}}+\varepsilon_{yy{\mathbf{k}}}}{2}\pm\sqrt{\left(\frac{\varepsilon_{d{\mathbf{k}}}-\varepsilon_{xx{\mathbf{k}}}-\varepsilon_{yy{\mathbf{k}}}}{2}\right)^{2}+J^{2}h_{pd}^{2}\left(v_{x{\mathbf{k}}}v_{x{\mathbf{k}}}^{*}+v_{y{\mathbf{k}}}v_{y{\mathbf{k}}}^{*}\right)}.$$

(28)

For $h_{pd}>0$, Eq. (15) can be rewritten as

$$1=\frac{2J}{N}\sum_{{\mathbf{k}}}\frac{v_{x{\mathbf{k}}}v_{x{\mathbf{k}}}^{*}+v_{y{\mathbf{k}}}v_{y{\mathbf{k}}}^{*}}{\varepsilon_{{\mathbf{k}}}^{+}-\varepsilon_{{\mathbf{k}}}^{-}}\,\left\{\theta(\varepsilon_{{\mathbf{k}}}^{+}-\mu)-\theta(\varepsilon_{{\mathbf{k}}}^{-}-\mu)\right\},$$

(29)

where $\theta(x)$ means the Heaviside step function.

In order to investigate the superconductivity in a strong coupling framework, I start with the Dyson-Gor’kov equations:

	$\displaystyle G_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})=G_{\mu\nu}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n})+G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Sigma}}_{\kappa\lambda}({\mathbf{k}},\mathrm{i}\epsilon_{n})G_{\lambda\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})+G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Phi}}_{\kappa\lambda}({\mathbf{k}},\mathrm{i}\epsilon_{n})F_{\lambda\nu}^{\dagger}(-{\mathbf{k}},-\mathrm{i}\epsilon_{n}),$
			(30)
	$\displaystyle F_{\mu\nu}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n})=G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Sigma}}_{\kappa\lambda}({\mathbf{k}},\mathrm{i}\epsilon_{n})F_{\lambda\nu}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n})+G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Phi}}_{\kappa\lambda}^{*}({\mathbf{k}},\mathrm{i}\epsilon_{n})G_{\lambda\nu}(-{\mathbf{k}},-\mathrm{i}\epsilon_{n}),$		(31)
	$\displaystyle F_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})=G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Sigma}}_{\kappa\lambda}({\mathbf{k}},\mathrm{i}\epsilon_{n})F_{\lambda\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})+G_{\mu\kappa}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}){\mathit{\Phi}}_{\kappa\lambda}({\mathbf{k}},\mathrm{i}\epsilon_{n})G_{\lambda\nu}(-{\mathbf{k}},-\mathrm{i}\epsilon_{n}).$		(32)

The orbital indices $\mu$, $\nu$, $\kappa$, and $\lambda$ run over $d$, $x$, and $y$, and I adopt the Einstein summation convention. $G_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ and $F_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ represent the normal and anomalous Green functions, respectively, and ${\mathit{\Sigma}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ and ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ correspond to the normal and anomalous self-energies, respectively. When ${\mathcal{H}}_{\mathrm{ex}}^{\prime}+{\mathcal{H}}_{\mathrm{pair}}^{\prime}+{\mathcal{H}}_{U}^{\prime}$ in Eq. (16) is treated as a perturbation, the normal self-energies up to the second order of $J$ and $U$ are evaluated by the IPT approximation as follows:

			$\displaystyle{\mathit{\Sigma}}_{dd}(\mathbf{k},\mathrm{i}\epsilon_{n})=\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[J^{2}\chi_{J}^{G}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})G_{pp}^{0}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})+U^{2}\chi_{U}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})G_{dd}^{0}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$		(36)
			$\displaystyle{\mathit{\Sigma}}_{d\alpha}(\mathbf{k},\mathrm{i}\epsilon_{n})=-v_{\alpha\mathbf{k}}^{*}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[J^{2}\chi_{J}^{G}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})G_{pd}^{0}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$
			$\displaystyle{\mathit{\Sigma}}_{\alpha d}(\mathbf{k},\mathrm{i}\epsilon_{n})=-v_{\alpha\mathbf{k}}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[J^{2}\chi_{J}^{G}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})G_{dp}^{0}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$
			$\displaystyle{\mathit{\Sigma}}_{\alpha\beta}(\mathbf{k},\mathrm{i}\epsilon_{n})=v_{\alpha\mathbf{k}}v_{\beta\mathbf{k}}^{*}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[J^{2}\chi_{J}^{G}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})G_{dd}^{0}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right].$

The IPT approximation was first applied in the study of the half-filled single-impurity Anderson model [31, 32], and it was adopted to solve the effective impurity model in the study of the $d=\infty$ Hubbard model [33, 34]. In these works, it was shown that the second order perturbation theory in large energy scale $U$ could reproduce not only the coherent band but also the lower and upper incoherent bands. In a later section, it will be shown that my approach can reproduce similar band structure to be justified as the theory for the 2D three-band $t$–$J$–$U$ model.

The anomalous self-energies up to the second order of $J$ and $U$ are evaluated as follows:

	$\displaystyle{\mathit{\Phi}}_{dd}(\mathbf{k},\mathrm{i}\epsilon_{n})=-\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[\left\{J+J^{2}\chi_{J}^{F}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})\right\}\!F_{pp}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})+\left\{U+U^{2}\chi_{U}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})\right\}\!F_{dd}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$
			(37)
	$\displaystyle{\mathit{\Phi}}_{d\alpha}(\mathbf{k},\mathrm{i}\epsilon_{n})=v_{\alpha-\mathbf{k}}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[\left\{J+J^{2}\chi_{J}^{F}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})\right\}\!F_{pd}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$		(38)
	$\displaystyle{\mathit{\Phi}}_{\alpha d}(\mathbf{k},\mathrm{i}\epsilon_{n})=v_{\alpha\mathbf{k}}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[\left\{J+J^{2}\chi_{J}^{F}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})\right\}\!F_{dp}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right],$		(39)
	$\displaystyle{\mathit{\Phi}}_{\alpha\beta}(\mathbf{k},\mathrm{i}\epsilon_{n})=-v_{\alpha\mathbf{k}}v_{\beta-\mathbf{k}}\frac{T}{N}\sum_{{\mathbf{k}}^{\prime}n^{\prime}}\left[\left\{J+J^{2}\chi_{J}^{F}({\mathbf{k}}-{\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n}-\mathrm{i}\epsilon_{n^{\prime}})\right\}\!F_{dd}({\mathbf{k}}^{\prime},\mathrm{i}\epsilon_{n^{\prime}})\right].$		(40)

Here, the orbital indices $\alpha$ and $\beta$ run over $x$ and $y$, and

	$\displaystyle\chi_{J}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})=\chi_{dd,pp}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})-\chi_{dp,dp}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})-\chi_{pd,pd}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})+\chi_{pp,dd}^{G}({\mathbf{q}},\mathrm{i}\omega_{m}),$		(41)
	$\displaystyle\chi_{J}^{F}({\mathbf{q}},\mathrm{i}\omega_{m})=\chi_{dd,pp}^{F}({\mathbf{q}},\mathrm{i}\omega_{m})-\chi_{dp,pd}^{F}({\mathbf{q}},\mathrm{i}\omega_{m})-\chi_{pd,dp}^{F}({\mathbf{q}},\mathrm{i}\omega_{m})+\chi_{pp,dd}^{F}({\mathbf{q}},\mathrm{i}\omega_{m}),$		(42)
	$\displaystyle\chi_{U}({\mathbf{q}},\mathrm{i}\omega_{m})=\chi_{dd,dd}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})+\chi_{dd,dd}^{F}({\mathbf{q}},\mathrm{i}\omega_{m}),$		(43)
	$\displaystyle\chi_{\mu\nu,\kappa\lambda}^{G}({\mathbf{q}},\mathrm{i}\omega_{m})=-\frac{T}{N}\sum_{{\mathbf{k}}n}G_{\mu\nu}^{0}({\mathbf{q}}+{\mathbf{k}},\mathrm{i}\omega_{m}+\mathrm{i}\epsilon_{n})G_{\kappa\lambda}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n}),$		(44)
	$\displaystyle\chi_{\mu\nu,\kappa\lambda}^{F}({\mathbf{q}},\mathrm{i}\omega_{m})=-\frac{T}{N}\sum_{{\mathbf{k}}n}F_{\mu\nu}({\mathbf{q}}+{\mathbf{k}},\mathrm{i}\omega_{m}+\mathrm{i}\epsilon_{n})F_{\kappa\lambda}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n}),$		(45)
	$\displaystyle G_{pp}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha\beta}v_{\alpha\mathbf{k}}^{*}v_{\beta\mathbf{k}}G_{\alpha\beta}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(46)
	$\displaystyle G_{dp}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha\mathbf{k}}G_{d\alpha}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(47)
	$\displaystyle G_{pd}^{0}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha\mathbf{k}}^{*}G_{\alpha d}^{0}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(48)
	$\displaystyle F_{pp}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha\beta}v_{\alpha\mathbf{k}}^{*}v_{\beta-\mathbf{k}}^{*}F_{\alpha\beta}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(49)
	$\displaystyle F_{dp}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha-\mathbf{k}}^{*}F_{d\alpha}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(50)
	$\displaystyle F_{pd}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha\mathbf{k}}^{*}F_{\alpha d}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(51)
	$\displaystyle F_{pp}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha\beta}v_{\alpha\mathbf{k}}v_{\beta-\mathbf{k}}F_{\alpha\beta}^{\dagger}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(52)
	$\displaystyle F_{dp}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha-\mathbf{k}}F_{d\alpha}^{\dagger}(\mathbf{k},\mathrm{i}\epsilon_{n}),{\mathrm{and}}$		(53)
	$\displaystyle F_{pd}^{\dagger}({\mathbf{k}},\mathrm{i}\epsilon_{n})=\sum_{\alpha}v_{\alpha\mathbf{k}}F_{\alpha d}^{\dagger}(\mathbf{k},\mathrm{i}\epsilon_{n}),$		(54)

using the boson Matsubara frequencies, $\omega_{m}=2m\pi T$ with integer $m$. In Eqs. (44) and (45), $\mu$, $\nu$, $\kappa$, and $\lambda$ denote $d$ or $p$, respectively. Note that $n_{d}$, $t$, $h_{pd}$, and the chemical potential $\mu$ must be determined self-consistently in the ground state of ${\mathcal{H}}_{0}^{\prime}$ through Eqs. (8)–(10), (28), and (29). To this end, I approximate $n_{d}$ by the number of $d$ electrons in the ground state of ${\mathcal{H}}_{0}^{\prime}$:

$$n_{d}=2-\frac{2}{N}\sum_{\mathbf{k}}\frac{1}{\varepsilon_{{\mathbf{k}}}^{+}-\varepsilon_{{\mathbf{k}}}^{-}}\left\{(\varepsilon_{{\mathbf{k}}}^{+}-\varepsilon_{xx{\mathbf{k}}}-\varepsilon_{yy{\mathbf{k}}})\theta(\varepsilon_{{\mathbf{k}}}^{+}-\mu)-(\varepsilon_{{\mathbf{k}}}^{-}-\varepsilon_{xx{\mathbf{k}}}-\varepsilon_{yy{\mathbf{k}}})\theta(\varepsilon_{{\mathbf{k}}}^{-}-\mu)\right\}.$$

(55)

Specifically, I regard $n_{d}$ as a given parameter and solve Eqs. (8)–(10), (28), (29), and (55) to determine $t$, $h_{pd}$, and the number of doped holes $\delta_{\mathrm{h}}^{0}$ for the ground state of ${\mathcal{H}}_{0}^{\prime}$, where

$$\delta_{\mathrm{h}}^{0}=\frac{2}{N}\sum_{\mathbf{k}}\left[\theta(\varepsilon_{{\mathbf{k}}}^{+}-\mu)+\theta(\varepsilon_{{\mathbf{k}}}^{-}-\mu)+\theta(-\mu)\right]-1.$$

Once $t$, $h_{pd}$, and $\delta_{\mathrm{h}}^{0}$ are determined for the ground state of ${\mathcal{H}}_{0}^{\prime}$, I treat $t$, $h_{pd}$, and $\delta_{\mathrm{h}}^{0}$ as temperature independent parameters, whose values do not change from those at $T=0$. Then, Eqs. (8)–(54) are solved in a fully self-consistent manner to obtain ${\mathit{\Sigma}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ and ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$. To determine the transition temperature $T_{\mathrm{c}}$, I perform these calculations in two steps. First, ${\mathit{\Sigma}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ is calculated with ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})=0$, and $\mu$ is self-consistently determined so that $\delta_{\mathrm{h}}$ obtained from $G_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ becomes equal to $\delta_{\mathrm{h}}^{0}$. In the first step, $\mu$ is correctly adjusted to compensate the temperature-dependent shift by ${\mathit{\Sigma}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ with ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})=0$. Here, $\delta_{\mathrm{h}}=n_{d\mathrm{h}}+n_{p\mathrm{h}}-1$, where

	$\displaystyle n_{d\mathrm{h}}$	$\displaystyle=$	$\displaystyle 2-2\,\frac{T}{N}\sum_{{\mathbf{k}}n}G_{dd}({\mathbf{k}},\mathrm{i}\epsilon_{n})e^{\mathrm{i}\epsilon_{n}{0^{+}}},$		(56)
	$\displaystyle n_{p\mathrm{h}}$	$\displaystyle=$	$\displaystyle 4-2\,\frac{T}{N}\sum_{{\mathbf{k}}n}\left[G_{xx}({\mathbf{k}},\mathrm{i}\epsilon_{n})+G_{yy}({\mathbf{k}},\mathrm{i}\epsilon_{n})\right]e^{\mathrm{i}\epsilon_{n}{0^{+}}},$		(57)

and $n_{d\mathrm{h}}$ and $n_{p\mathrm{h}}$ are the number of $d$ and $p$ holes, respectively. Next, using the determined $\mu$, fully self-consistent calculations are performed to obtain ${\mathit{\Sigma}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ and ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$. At this time, only the temperature-dependent shift by ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ is reflected in $\delta_{\mathrm{h}}$ obtained from $G_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$. That is, if $\delta_{\mathrm{h}}$ deviates from $\delta_{\mathrm{h}}^{0}$, ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})\neq 0$. Therefore, the temperature at which $\delta_{\mathrm{h}}$ deviates from $\delta_{\mathrm{h}}^{0}$ is $T_{\mathrm{c}}$. Also in the second step, $\mu$ can be self-consistently determined so that $\delta_{\mathrm{h}}$ obtained from $G_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})$ becomes equal to $\delta_{\mathrm{h}}^{0}$ with ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})\neq 0$. In this case, the temperature at which $\mu$ deviates from the value with ${\mathit{\Phi}}_{\mu\nu}({\mathbf{k}},\mathrm{i}\epsilon_{n})=0$ is $T_{\mathrm{c}}$, which is consistent with the temperature at which $\delta_{\mathrm{h}}$ deviates from $\delta_{\mathrm{h}}^{0}$ with fixed $\mu$.