# A More Complete Model for High-Temperature Superconductors

Tao Sun

► To cite this version:

Tao Sun. A More Complete Model for High-Temperature Superconductors. 2017. hal-01502648

HAL Id: hal-01502648

https://hal.science/hal-01502648

Preprint submitted on 5 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A More Complete Model for High-Temperature Superconductors

Tao Sun

tao@taosphysics.net

April 5, 2017

## Abstract

To date, the Hubbard model and its strong coupling limit, the t-J model, serve as the canonical model for strongly correlated electron systems in solids. Approximating the Coulomb interaction by only the on-site term (Hubbard U-term), however, may not be sufficient to describe the essential physics of interacting electron systems. We develop a more complete model in which all the next leading order terms besides the on-site term are retained. Moreover, we discuss how the inclusion of these neglected interaction terms in the Hubbard model changes the t-J model.

## 1 Introduction

It is generally believed [1] that some essential physics of strongly correlated electron systems in solids, such as high-temperature superconductors [2], can be described by the Hubbard model [3, 4]

$$H = -t \sum_{ij,\sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow}, \quad (1)$$

where $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the creation and annihilation operators for an electron of spin $\sigma$ at site $i$, and $n_{i\sigma} = c_{i\sigma}^{\dagger}c_{i\sigma}$, the corresponding number operator of electrons. In Eq (1), $\gamma_{ij} = 1$ for nearest-neighbor (n.n.) sites $i, j$, and 0 otherwise, which restricts the summation over the nearest-neighbor pairs. The t-term describes the kinetic energy due to hopping between adjacent sites, and the U-term is the on-site approximation of the Coulomb interaction between electrons. Indeed, the Hubbard model offers the simplest way to describe both the band motion and atomic features of the interacting electron systems in solids.

In the Hubbard model, there are four possible states on the lattice sites: empty state $|0\rangle$, singly occupied states $|\uparrow\rangle$ (spin up) and $|\downarrow\rangle$ (spin down), and doubly occupied state $|\uparrow\downarrow\rangle$. Apparently, the double occupancy state $|\uparrow\downarrow\rangle$ is not favorable when U is large. It is commonly believed that a model in which all the doubly occupied states are excluded can describe the low-energy physics
of the system. Technically, such a singly occupied state system (with only empty and singly occupied sites) can be obtained by the Gutzwiller projection operator approach, which projects the Hubbard model to a subspace of the original Hilbert space. For large $U$, to the order of $J = 2t^2/U$ in the perturbation calculation, this projecting scheme leads to the well-known $t-J$ model [5, 6]

$$H = -t \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + J \sum_{ij} \gamma_{ij} (S_i \cdot S_j - \frac{1}{4} \tilde{n}_i \tilde{n}_j), \quad (2)$$

where $\tilde{c}_{i\sigma}^{\dagger} = c_{i\sigma}^{\dagger}(1 - n_{i\sigma})$, $\tilde{c}_{i\sigma} = c_{i\sigma}(1 - n_{i\sigma})$, $\tilde{n}_i = \sum_{\sigma} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{i\sigma}$, and $S_i$ is the spin operator at site $i$. Note that the projection operator also generates a so-called three-site term that has been neglected in the $t-J$ model. That is, Eq (2) is in fact a model with only two-site terms. The $t-J$ model can be viewed as a large $U$ limit of the Hubbard model, and as a prototype to describe the low-energy physics of the interacting electron systems. Using the standard Green function technique, one can indeed show that the $t-J$ model does exhibit a state with nonzero $\langle c_{i\downarrow} c_{j\uparrow} \rangle$, signaling the occurrence of superconductivity. Nonetheless, it is still somewhat difficult to obtain a more complete understanding of high-temperature superconductors from the $t-J$ model.

The Hubbard model is the simplest approximation of the general Hamiltonian of interacting electron systems, in which all the interaction terms except the on-site term ($U$-term) are neglected. It is possible that the neglected interaction terms may play important roles in the understanding of the physics of the strongly correlated electron systems. Some important physics may be omitted in the Hubbard model. As a result, the $t-J$ model given in Eq (2), which is a projection to the subspace of singly occupied states of the Hubbard model, may not be sufficient to describe the essential physics of the interacting electrons.

In this paper, we will first describe the details of establishing a generalized Hubbard model, in which all of the two-site interaction terms between nearest-neighbor sites are retained. We will then apply the Gutzwiller projection scheme to this more general model to obtain an effective Hamiltonian in the subspace of singly occupied states. Moreover, we will discuss how the inclusion of these neglected interaction terms in the Hubbard model modifies the $t-J$ model.

## 2 The Model

We now present a systematic derivation of a more complete effective Hamiltonian with all two-site interaction terms between nearest-neighbor sites in the subspace of singly occupied states.
## 2.1 Bloch Representation

The general Hamiltonian describing the dynamics of the electrons in solids can be expressed as follows [3, 7, 8] in Bloch representation

$$H = \sum_{\mathbf{k},\sigma} \epsilon_{\mathbf{k}} c_{\mathbf{k}\sigma}^{\dagger} c_{\mathbf{k}\sigma} + \frac{1}{2} \sum_{\mathbf{k}_1 \mathbf{k}_2 \mathbf{k}'_1 \mathbf{k}'_2, \sigma, \sigma'} \langle \mathbf{k}_1 \mathbf{k}_2 | \frac{1}{r} | \mathbf{k}'_1 \mathbf{k}'_2 \rangle c_{\mathbf{k}_1 \sigma}^{\dagger} c_{\mathbf{k}_2 \sigma'}^{\dagger} c_{\mathbf{k}'_2 \sigma'} c_{\mathbf{k}'_1 \sigma}, \quad (3)$$

where $c_{\mathbf{k},\sigma}^{\dagger}$ and $c_{\mathbf{k},\sigma}$ are the creation and annihilation operators for electrons in the Bloch state with wave vector $\mathbf{k}$ and spin $\sigma$, $\epsilon_{\mathbf{k}}$ is the band energy, $\mathbf{k}$ runs over the first Brillouin zone, and the matrix element

$$\langle \mathbf{k}_1 \mathbf{k}_2 | \frac{1}{r} | \mathbf{k}'_1 \mathbf{k}'_2 \rangle = e^2 \int d\mathbf{x} d\mathbf{x}' \frac{\psi_{\mathbf{k}_1}^{*}(\mathbf{x}) \psi_{\mathbf{k}'_1}(\mathbf{x}) \psi_{\mathbf{k}_2}^{*}(\mathbf{x}') \psi_{\mathbf{k}'_2}(\mathbf{x}')}{|\mathbf{x} - \mathbf{x}'|}, \quad (4)$$

where $\psi_{\mathbf{k}}(\mathbf{x})$ and $\psi_{\mathbf{k}}^{*}(\mathbf{x})$ are the Bloch functions of the energy band. In Eq (3), the first term represents the kinetic energies of the electrons in the energy bands, and the second term is the Coulomb interaction energy of electrons.

## 2.2 Wannier Representation

The Hamiltonian in Eq (3) can be transformed to its Wannier representation form

$$H = \sum_{ij,\sigma} T_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + \frac{1}{2} \sum_{ijkl,\sigma\sigma'} \langle ij | \frac{1}{r} | kl \rangle c_{i\sigma}^{\dagger} c_{j\sigma'}^{\dagger} c_{l\sigma'} c_{k\sigma}, \quad (5)$$

where $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the creation and annihilation operators for an electron with spin $\sigma$ in a Wannier orbital localized at site $i$, $T_{ij}$ is the Fourier transform of the band energy $\epsilon_{\mathbf{k}}$

$$T_{ij} = \frac{1}{N} \sum_{\mathbf{k}} \epsilon_{\mathbf{k}} e^{i\mathbf{k} \cdot (\mathbf{R}_i - \mathbf{R}_j)}, \quad (6)$$

and the Wannier representation matrix element is given by

$$\langle ij | \frac{1}{r} | kl \rangle = e^2 \int d\mathbf{x} d\mathbf{x}' \frac{\phi^*(\mathbf{x} - \mathbf{R}_i) \phi(\mathbf{x} - \mathbf{R}_k) \phi^*(\mathbf{x}' - \mathbf{R}_j) \phi(\mathbf{x}' - \mathbf{R}_l)}{|\mathbf{x} - \mathbf{x}'|}, \quad (7)$$

where $\phi(\mathbf{x} - \mathbf{R}_i)$ and $\phi^*(\mathbf{x} - \mathbf{R}_i)$ are the Wannier functions localized around lattice site $i$.

In the most general cases where $i \neq j \neq k \neq l$, the Wannier matrix elements $\langle ij | \frac{1}{r} | kl \rangle$ are the so-called four-center integrals, and the corresponding terms in the series of the Coulomb interaction in Eq (5) can be referred to as the four-site terms. Since the Wannier function $\phi(\mathbf{x} - \mathbf{R}_i)$ goes to zero rapidly when $\mathbf{x}$ is away from $\mathbf{R}_i$, the matrix element (7) is not negligible only when the sites $i, j, k, l$ are close enough so that the overlap between the Wannier
functions is sufficiently large. Therefore, the Coulomb interaction energy can be conveniently approximated by a number of its leading terms of the series in the Wannier representation expression given in Eq (5). The biggest interaction term is the so-called on-site term ($i = j = k = l$)

$$ \frac{1}{2} U \sum_{i,\sigma} n_{i\sigma} n_{i\bar{\sigma}} = U \sum_{i} n_{i\uparrow} n_{i\downarrow}, \qquad (8) $$

where $U = \langle ii|_{\bar{r}}|ii\rangle$. If only this on-site term is taken into account, the Hubbard model is obtained. The next leading order terms are the ones where the set $\{i, j, k, l\}$ contains only two distinct nearest-neighbor sites. These terms may be referred to as the two-site interaction terms [9]. In this paper, we take one step further beyond the Hubbard model: In the approximation of the Coulomb interaction, we will retain all of these two-site interaction terms, in addition to the Hubbard on-site term.

## 2.3 Two-Site Interaction Approximation

We now describe a generalized Hubbard model with all of the two-site interaction terms between two nearest-neighbor sites. The two-side interaction terms are those where the set $\{i, j, k, l\}$ consists of only one pair of nearest-neighbor sites. These terms fall under two categories:

1. Among the four sites $i, j, k, l$, three of them are in fact the same site, and the fourth one is a nearest neighbor of the other three. There are four such possibilities:

a. $i = j = k$, and $l$ is the n.n. site of $i$:

$$ \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ii|_{\bar{r}} \frac{1}{r} |ij\rangle n_{i\sigma} c_{i\bar{\sigma}}^{\dagger} c_{j\bar{\sigma}} = \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ii|_{\bar{r}} \frac{1}{r} |ij\rangle n_{i\bar{\sigma}} c_{i\sigma}^{\dagger} c_{j\sigma} \qquad (9) $$

b. $i = j = l$, and $k$ is the n.n. site of $i$:

$$ \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ii|_{\bar{r}} \frac{1}{r} |ji\rangle n_{i\bar{\sigma}} c_{i\sigma}^{\dagger} c_{j\sigma} \qquad (10) $$

c. $i = l = k$, and $j$ is the n.n. site of $i$:

$$ \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ij|_{\bar{r}} \frac{1}{r} |ii\rangle n_{i\sigma} c_{j\bar{\sigma}}^{\dagger} c_{i\bar{\sigma}} = \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ji|_{\bar{r}} \frac{1}{r} |jj\rangle n_{j\bar{\sigma}} c_{i\sigma}^{\dagger} c_{j\sigma} \qquad (11) $$

d. $j = l = k$, and $i$ is the n.n. site of $j$:

$$ \frac{1}{2} \sum_{ij,\sigma} \gamma_{ij} \langle ij|_{\bar{r}} \frac{1}{r} |jj\rangle n_{j\bar{\sigma}} c_{i\sigma}^{\dagger} c_{j\sigma} \qquad (12) $$

Eq (9) is true because the spin summation index $\sigma$ can be replaced by $\bar{\sigma}$ in the expression. Similarly, the right-hand-side of Eq (11) results from interchanging
the summation indexes i and j, in addition to the change of $\sigma \to \bar{\sigma}$. From Eq (7), it can be verified that all the matrix elements in Eqs (9-12) are equal. Thus, all of the terms above can be combined into the following expression

$$X \sum_{ij, \sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} (n_{i\sigma} + n_{j\sigma}), \qquad (13)$$

where the parameter X represents the matrix elements

$$X = \langle ii|\frac{1}{r}|ij\rangle = \langle ii|\frac{1}{r}|ji\rangle = \langle ij|\frac{1}{r}|ii\rangle = \langle ji|\frac{1}{r}|ii\rangle. \qquad (14)$$

Here, i and j represent the nearest-neighbor sites.

2. Among the sites i, j, k, l, there are two identical pairs, and moreover, the said two pairs are nearest neighbors. There are three such possibilities:

a. i = j, l = k, i and l are the n.n. sites:

$$\frac{1}{2} \sum_{ij, \sigma} \gamma_{ij} \langle ii| \frac{1}{r} |jj\rangle c_{i\sigma}^{\dagger} c_{j\sigma} c_{i\bar{\sigma}}^{\dagger} c_{j\bar{\sigma}} \qquad (15)$$

b. i = k, j = l, i and j are the n.n. sites:

$$\frac{1}{2} \sum_{ij, \sigma\sigma'} \gamma_{ij} \langle ij|\frac{1}{r}|ij\rangle n_{i\sigma} n_{j\sigma'} \qquad (16)$$

c. i = l, j = k, i and j are the n.n. sites:

$$\frac{1}{2} \sum_{ij, \sigma\sigma'} \gamma_{ij} \langle ij|\frac{1}{r}|ji\rangle c_{i\sigma}^{\dagger} c_{j\sigma'}^{\dagger} c_{i\sigma'} c_{j\sigma} \qquad (17)$$

It has been argued that $\langle ii|\frac{1}{r}|jj\rangle$ and $\langle ij|\frac{1}{r}|ji\rangle$ are in the same order, but much smaller than $\langle ij|\frac{1}{r}|ij\rangle$ [3, 10]. Therefore, summing over all of the three terms in Eqs (15-17), we have the following terms in this category

$$\frac{1}{2} V \sum_{ij, \sigma\sigma'} \gamma_{ij} n_{i\sigma} n_{j\sigma'} + \frac{1}{2} Y \sum_{ij, \sigma} \gamma_{ij} (c_{i\sigma}^{\dagger} c_{j\sigma} c_{i\bar{\sigma}}^{\dagger} c_{j\bar{\sigma}} + \sum_{\sigma'} c_{i\sigma}^{\dagger} c_{j\sigma'}^{\dagger} c_{i\sigma'} c_{j\sigma}) \qquad (18)$$

where the parameters V and Y represent the matrix elements

$$V = \langle ij|\frac{1}{r}|ij\rangle, \qquad (19)$$

$$Y = \langle ii|\frac{1}{r}|jj\rangle = \langle ij|\frac{1}{r}|ji\rangle.$$

Eqs (13) and (18) contain all the two-site Coulomb interaction terms between two nearest neighbors.

Similarly, for the kinetic energy term, only the constant and nearest-neighbor hopping terms are considered for consistency

$$T_{ij} = \frac{1}{N} \sum_{\mathbf{k}} \epsilon_{\mathbf{k}} e^{i\mathbf{k} \cdot (\mathbf{R}_i - \mathbf{R}_j)} = T_0 \delta_{ij} + T_1 \gamma_{ij} + \dots, \qquad (20)$$
where

$$ 
\begin{align} 
T_0 &= T_{ii} = \frac{1}{N} \sum_{\mathbf{k}} \epsilon_{\mathbf{k}}, \\ 
T_1 &= T_{\langle ij \rangle} = \frac{1}{N} \sum_{\mathbf{k}} \epsilon_{\mathbf{k}} e^{i\mathbf{k} \cdot (\mathbf{R}_i - \mathbf{R}_j)}, \quad \text{for n.n. sites } \langle ij \rangle, 
\end{align} 
\tag{21} $$

and (...) represents the hopping terms beyond nearest-neighbor sites. Since $T_1 < 0$, we will use $t = -T_1$ in the expressions hereafter.

Substituting Eq (20) into Eq (5), and approximating the Coulomb interaction by all the terms given in Eqs (8), (13), and (18), we have

$$ 
\begin{align} 
H \approx T_0 \sum_{i,\sigma} n_{i\sigma} - t \sum_{ij,\sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow} \nonumber \\ 
+ \frac{1}{2} V \sum_{ij,\sigma\sigma'} \gamma_{ij} n_{i\sigma} n_{j\sigma'} + X \sum_{ij,\sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} (n_{i\bar{\sigma}} + n_{j\bar{\sigma}}) \tag{22} \\ 
+ \frac{1}{2} Y \sum_{ij,\sigma} \gamma_{ij} (c_{i\sigma}^{\dagger} c_{j\sigma} c_{i\bar{\sigma}}^{\dagger} c_{j\bar{\sigma}} + \sum_{\sigma'} c_{i\sigma}^{\dagger} c_{j\sigma'}^{\dagger} c_{i\sigma'} c_{j\sigma}). \nonumber 
\end{align} 
$$

Eq (22) contains all the two-site interaction terms as well as the on-site term and can be considered the most complete generalization of the Hubbard model to date. The terms parametrized by X and V have been discussed in literature [10, 11]. The X-term is the so-called bond-change interaction, which describes the density-dependent nearest-neighbor hopping of electrons. The V-term is the so-called nearest-neighbor density interaction. The matrix element Y was discussed in reference [3]. When $X = V = Y = 0$, Eq (22) reduced to the Hubbard model.

It is worth pointing out that, in contrast to the U-term, the terms parametrized by V, X, and Y retain the long-range nature of the Coulomb interaction. Mott argued that, based on the long-range nature of the Coulomb interaction, the correlated electron system undergoes a first-order metal-insulator transition triggered by electronic correlations [12]. The Hubbard model, however, only exhibits a continuous metal-insulator transition owing to the correlation effect of electrons [3]. The reason that the Hubbard model does not predict a discontinuity in the number of current carriers may be attributed to the lack of the long-range interactions [3, 12, 13]. Therefore, a discontinuous metal-insulator transition may exist in the model given in Eq (22) [14].

## 2.4 Gutzwiller Projection Technique

Similar to the Hubbard model, in the model described by Eq (22), there are four possible states on the lattice sites: empty, singly occupied with spin up, singly occupied with spin down, and doubly occupied. The Hamiltonian describing the low-energy physics should not include the doubly occupied states. We need to project the Hamiltonian (22) to a subspace of the Hilbert space where only
empty and singly occupied sites are allowed. This can be accomplished by the Gutzwiller projection approach. The Gutzwiller projection operators can be defined in the following form [4, 15]

$$ P_s = \prod_i (1 - n_{i\uparrow} n_{i\downarrow}), \qquad (23) \\ P_d = I - P_s, $$

where $n_{i\uparrow}$ and $n_{i\downarrow}$ are the electron number operators at site $i$ with spin up and down, respectively, $i$ runs over all the lattice sites, and $I$ is the identity operator. It can easily be verified that the projection operators $P_s$ and $P_d$ have the following properties

$$ \begin{gathered} (P_s)^2 = P_s, \quad (P_s)^n = P_s, \\ (P_d)^2 = P_d, \quad (P_d)^n = P_d, \qquad (24) \\ P_s P_d = 0, \end{gathered} $$

where $n$ is any positive integer. The Gutzwiller projection operators can partition the original Hilbert space into two subspaces. The first one (denoted by $S$) contains the states with only empty and singly occupied sites, while the second one (denoted by $D$) contains the states with at least one doubly occupied site. The subscripts $s$ and $d$ denote the two subspaces $S$ and $D$, respectively.

We now derive the effective Hamiltonian in the subspace $S$ using the Gutzwiller projection method. The specific treatment follows the formalism developed in Ref. [16]. The Schrödinger eigenvalue equation of the Hamiltonian $H$ given in Eq (22) can be written as

$$ H|\Psi\rangle = \varepsilon|\Psi\rangle, \qquad (25) $$

where $|\Psi\rangle$ is the wave vector function in Hilbert space and $\varepsilon$ its corresponding energy eigenvalue. Since $P_s + P_d = I$, Eq (25) may be written as

$$ H(P_s + P_d)|\Psi\rangle = \varepsilon(P_s + P_d)|\Psi\rangle. \qquad (26) $$

Operating $P_s$ and $P_d$ respectively to Eq (26) from the left, and taking the properties of $P_s$ and $P_d$ into account, we obtain

$$ \begin{bmatrix} P_s H P_s - \varepsilon & P_s H P_d \\ P_d H P_s & P_d H P_d - \varepsilon \end{bmatrix} \begin{bmatrix} |\Psi_s\rangle \\ |\Psi_d\rangle \end{bmatrix} = 0, $$

where $|\Psi_s\rangle = P_s|\Psi\rangle$ and $|\Psi_d\rangle = P_d|\Psi\rangle$, the wave vectors in the subspaces $S$ and $D$, respectively. The projected Hamiltonian matrix above can be diagonalized as

$$ \begin{bmatrix} H_s - \varepsilon & 0 \\ 0 & H_d - \varepsilon \end{bmatrix} \begin{bmatrix} |\Psi_s\rangle \\ |\Psi_d\rangle \end{bmatrix} = 0, $$

where

$$ H_s = P_s H P_s + P_s H P_d (\varepsilon - P_d H P_d)^{-1} P_d H P_s, \qquad (27) $$

$$ H_d = P_d H P_d + P_d H P_s (\varepsilon - P_s H P_s)^{-1} P_s H P_d, \qquad (28) $$
are the projected Hamiltonians in the subspaces. In the following discussion, we are only interested in $H_s$, the effective Hamiltonian in the singly occupied subspace.

Operating the projection operators $P_s$ and $P_d$ on $H$ in Eq (22) from both the left and right, after some straightforward manipulations, we can obtain the following projected Hamiltonian matrix elements

$$P_sHP_s = T_0 \sum_{i,\sigma} \tilde{n}_{i\sigma} - t \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \frac{1}{2} V \sum_{ij,\sigma\sigma'} \gamma_{ij} \tilde{n}_{i\sigma} \tilde{n}_{j\sigma'} \\ + \frac{1}{2} Y \sum_{ij,\sigma} \gamma_{ij} (c_{i\sigma}^{\dagger} c_{j\sigma} c_{j\sigma}^{\dagger} c_{i\sigma} - \tilde{n}_{i\sigma} \tilde{n}_{j\sigma}), \quad (29)$$

$$P_sHP_d = -(t-X) \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} c_{j\sigma} n_{j\sigma}, \quad (30)$$

$$P_dHP_s = -(t-X) \sum_{ij,\sigma} \gamma_{ij} n_{i\sigma} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma}, \quad (31)$$

$$P_dHP_d = -(t-2X) \sum_{ij,\sigma} \gamma_{ij} n_{i\sigma} \tilde{c}_{i\sigma}^{\dagger} c_{j\sigma} n_{j\sigma} + (U+2T_0) \sum_i n_{i\uparrow} n_{i\downarrow} \\ + (2V-Y) \sum_{ij} \gamma_{ij} n_{i\uparrow} n_{i\downarrow} (n_{j\uparrow} + n_{j\downarrow} - n_{j\uparrow} n_{j\downarrow}) \\ + \frac{1}{2} Y \sum_{ij,\sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} c_{i\sigma}^{\dagger} c_{j\sigma}, \quad (32)$$

where $\tilde{c}_{i\sigma}^{\dagger} = c_{i\sigma}^{\dagger}(1-n_{i\sigma})$, $\tilde{c}_{i\sigma} = c_{i\sigma}(1-n_{i\sigma})$, and $\tilde{n}_i = \sum_{\sigma} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{i\sigma}$ are the projected creation, annihilation, and number operators of electrons. Note that the V- and Y- terms do contribute to the zeroth-order of the effective Hamiltonian $P_sHP_s$. These are the terms neglected in the t-J model. Note also that the X and U-terms do not survive in the expression of $P_sHP_s$, indicating that they are not operators of the subspace $S$. As we will see below, however, they do contribute to the first-order terms in the perturbation expansion when the projected Hamiltonian matrix is diagonalized.

With the results above, we are now ready to obtain the effective Hamiltonian $H_s$. For convenience, we rewrite

$$H_s = H_s^{(0)} + H_s^{(1)}, \quad (33)$$

$$H_s^{(0)} = P_sHP_s, \quad (34)$$

$$H_s^{(1)} = P_sHP_d(\varepsilon - P_dHP_d)^{-1}P_dHP_s. \quad (35)$$

The zeroth-order term $H_s^{(0)}$ has been obtained in Eq (29). In order to construct the effective Hamiltonian $H_s$ in the subspace $S$, the key operator to be determined is the inverse operator $(\varepsilon - P_dHP_d)^{-1}$. Since this is an operator in the subspace $D$, using the completeness of the wave vectors $|\Psi_d\rangle$, we have

$$(\varepsilon - P_dHP_d)^{-1} = \sum_{r,r'} |\Psi_d^r\rangle\langle\Psi_d^r|(\varepsilon - P_dHP_d)^{-1}|\Psi_d^{r'}\rangle\langle\Psi_d^{r'}|, \quad (36)$$
where the superscripts r and r' are the labels of the states in the subspace D. Please note that these states in D must be virtual because this inverse operator is in the middle of $H_s^{(1)}$, which is an operator of S. As can be seen from Eq (32), the expression of $P_dHP_d$ is complicated. Nevertheless, in the case where there is one virtual doubly occupied site [15], the matrix element

$$ \langle \Psi_d^r | (\varepsilon - P_d H P_d)^{-1} | \Psi_d^{r'} \rangle = \langle \Psi_d^r | [\varepsilon - 2T_0 - U' + \sum_{ij,\sigma} \gamma_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} \left( (t-2X)n_{i\bar{\sigma}} n_{j\bar{\sigma}} - \frac{1}{2} Y c_{i\bar{\sigma}}^{\dagger} c_{j\bar{\sigma}} \right)]^{-1} | \Psi_d^{r'} \rangle, \quad (37) $$

where $U' = U + z'(2V - Y)$, $0 \le z' \le z$, and z is the coordinate number. If we consider the parameters $\varepsilon - 2T_0$, $t - 2X$, and Y to be small quantities compared with $U'$, we can perform a perturbation expansion with respect to $1/U'$. To the leading term in $1/U'$, we have

$$ \begin{aligned} & (\varepsilon - P_d H P_d)^{-1} \\ &= \sum_{r,r'} |\Psi_d^r\rangle\langle\Psi_d^r| \left[ -\frac{1}{U'} \left( 1 + O\left(\frac{\varepsilon - 2T_0}{U'}, \frac{t-2X}{U'}, \frac{Y}{U'}\right) \right)^{-1} \right] |\Psi_d^{r'}\rangle\langle\Psi_d^{r'}| \\ &= \left[ -\frac{1}{U'} \left( 1 + O\left(\frac{\varepsilon - 2T_0}{U'}, \frac{t-2X}{U'}, \frac{Y}{U'}\right) \right)^{-1} \right] \sum_{r,r'} |\Psi_d^r\rangle\langle\Psi_d^{r'}| \delta_{r,r'} \\ &= -\frac{1}{U'} \left( 1 + O\left(\frac{\varepsilon - 2T_0}{U'}, \frac{t-2X}{U'}, \frac{Y}{U'}\right) \right)^{-1} \\ &= -\frac{1}{U'} + O\left(\frac{\varepsilon - 2T_0}{U'^2}, \frac{t-2X}{U'^2}, \frac{Y}{U'^2}\right). \end{aligned} \quad (38) $$

Substituting (38) into Eq (35) and using the expressions of $P_sHP_d$ and $P_dHP_s$ in Eqs (30) and Eqs (31), we obtain for the second term of $H_s$ in the form

$$ H_s^{(1)} = -J_0 \sum_{ij,\sigma} \sum_{i'j',\sigma'} \gamma_{ij} \gamma_{i'j'} \tilde{c}_{i\sigma}^{\dagger} c_{j\sigma} n_{j\bar{\sigma}} n_{i'\bar{\sigma}} c_{i'\sigma'}^{\dagger} \tilde{c}_{j'\sigma'}, \quad (39) $$

where $J_0 = (t - X)^2/U'$. In general, Eq (39) represents a term with four sites $i, j, i', j'$ involved. In the most general case where $i \neq j \neq i' \neq j'$, it is a four-site term and an operator of the subspace D. Therefore we examine its special cases where two of the four sites are in fact the same site. Since $<ij>$ and $<i'j'>$ are two nearest-neighbor pairs, there are only four such possibilities:

1. i and i' are the same site

2. i and j' are the same site

3. j and i' are the same site

4. j and j' are the same site

Further examination reveals that cases 1 and 4 result in a zero operator, while case 2 is still an operator in the subspace D. Only case 3 with $i' = j$ provides
a meaningful contribution to the operator in the subspace S

$$H_s^{(1)} = -J_0 \sum_{ijk,\sigma} \gamma_{ij} \gamma_{jk} \tilde{c}_{i\sigma}^{\dagger} (\tilde{n}_{j\sigma} \tilde{c}_{k\sigma} + c_{j\sigma} \tilde{c}_{j\sigma}^{\dagger} \tilde{c}_{k\sigma}). \quad (40)$$

Please note that in Eq (40) the site j is the nearest neighbor of both sites i and k. Therefore, the two-site term can only result from the case of i = k in the summation. Writing the two-site and three-site terms separately, we have

$$H_s^{(1)} = (H_s^{(1)})_{2site} + (H_s^{(1)})_{3site}, \quad (41)$$

where

$$(H_s^{(1)})_{2site} = -J_0 \sum_{ij,\sigma} \gamma_{ij} (\tilde{n}_{i\sigma} \tilde{n}_{j\sigma} + c_{i\sigma}^{\dagger} c_{j\sigma} c_{j\sigma}^{\dagger} c_{i\sigma}), \quad (42)$$

$$(H_s^{(1)})_{3site} = -J_0 \sum_{i \neq k, j, \sigma} \gamma_{ij} \gamma_{jk} \tilde{c}_{i\sigma}^{\dagger} (\tilde{n}_{j\sigma} \tilde{c}_{k\sigma} + c_{j\sigma} \tilde{c}_{j\sigma}^{\dagger} \tilde{c}_{k\sigma}). \quad (43)$$

Please note that the summation in Eq (43) excludes the case i = k. This three-site term is usually ignored in the discussion of the t-J model. We will also neglect this three-site term in the subsequent discussions.

Substituting the expression of $H_s^{(0)}$ in Eq (29) and $H_s^{(1)} \approx (H_s^{(1)})_{2site}$ in Eq (42) into Eq (33), we finally obtain, to the order of two-site interactions terms, the effective Hamiltonian $H_s$

$$H_s = T_0 \sum_{i\sigma} \tilde{n}_{i\sigma} + \sum_{ij,\sigma} \gamma_{ij} \left( -t \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} - \frac{1}{2}(2J_0 - Y) c_{i\sigma}^{\dagger} c_{j\sigma}^{\dagger} c_{i\sigma} c_{j\sigma} \right. \\ \left. + \frac{1}{2}(V - Y) \tilde{n}_{i\sigma} \tilde{n}_{j\sigma} + \frac{1}{2}(V - 2J_0) \tilde{n}_{i\sigma} \tilde{n}_{j\sigma} \right), \quad (44)$$

which can be written as

$$H_s = T_0 \sum_{i\sigma} \tilde{n}_{i\sigma} + \sum_{ij,\sigma} \gamma_{ij} \left( -t \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} - \frac{1}{2} J c_{i\sigma}^{\dagger} c_{j\sigma}^{\dagger} c_{i\sigma} c_{j\sigma} \right. \\ \left. + \frac{1}{2} p \tilde{n}_{i\sigma} \tilde{n}_{j\sigma} + \frac{1}{2} q \tilde{n}_{i\sigma} \tilde{n}_{j\sigma} \right), \quad (45)$$

with the coefficients given by the following equations

$$\begin{align} p &= V - Y, \nonumber \\ q &= V - 2J_0, \quad J_0 = (t - X)^2/U', \quad U' = U + z'(2V - Y), \tag{46} \\ J &= p - q, \nonumber \end{align}$$

where $0 \le z' \le z$ as discussed above. Eq (45) is the effective Hamiltonian in the subspace of singly occupied states. Since it is projected from a more general Hamiltonian retaining all terms up to the two-site interactions, we believe that Eq (45) can describe the low-energy physics of the interacting electron systems in solids more accurately and completely.
Please note that by rescaling $H_s$ with $t$, this model has only three dimensionless independent parameters $\bar{T}_0 = T_0/t$, $\bar{p} = p/t$, and $\bar{q} = q/t$. One may expect that the system will exhibit ferromagnetism when $\bar{q}$ is sufficiently larger than $\bar{p}$. On the other hand, the system will exhibit antiferromagnetism when $\bar{p}$ is sufficiently larger than $\bar{q}$. Our studies demonstrate that this is true. More interestingly, our study also shows that this model exhibits superconductivity for some parameter range of $\bar{p} > \bar{q}$.

## 3 The t-J Model

Apparently, Eq (45) reduces to the t-J model with $J = 2J_0$ in the case of $X = Y = V = 0$. We will show that when $Y = V$, Eq (45) will reduce to the t-J model as well (with $J = 2J_0 - V$ in this case). Making use of the following operator identities

$$C_{i\uparrow}^{\dagger} C_{i\downarrow} = S_{i}^{+} = S_{i}^{x} + i S_{i}^{y}, \qquad (47)$$

$$C_{i\downarrow}^{\dagger} C_{i\uparrow} = S_{i}^{-} = S_{i}^{x} - i S_{i}^{y}, \qquad (48)$$

$$n_{i\uparrow} - n_{i\downarrow} = 2S_{i}^{z}, \qquad (49)$$

we obtain

$$\sum_{\sigma} c_{i\sigma}^{\dagger} c_{j\sigma} c_{j\bar{\sigma}}^{\dagger} c_{i\bar{\sigma}} = -2(S_{i}^{x} S_{j}^{x} + S_{i}^{y} S_{j}^{y}), \qquad (50)$$

$$\sum_{\sigma} (\tilde{n}_{i\sigma}\tilde{n}_{j\sigma} - \tilde{n}_{i\sigma}\tilde{n}_{j\bar{\sigma}}) = 4S_{i}^{z}S_{j}^{z}, \qquad (51)$$

where $S_i^\alpha$, $\alpha = x, y, z$, are the spin operators at site $i$. Furthermore, it is easy to verify that

$$\sum_{\sigma} (\tilde{n}_{i\sigma}\tilde{n}_{j\sigma} + \tilde{n}_{i\sigma}\tilde{n}_{j\bar{\sigma}}) = \tilde{n}_{i}\tilde{n}_{j}, \qquad (52)$$

where $\tilde{n}_i = \sum_{\sigma} \tilde{n}_{i\sigma}$. Eqs (50)-(52) immediately lead to the following well-known relation

$$\sum_{\sigma} (c_{i\sigma}^{\dagger} c_{j\sigma} c_{j\bar{\sigma}}^{\dagger} c_{i\bar{\sigma}} + \tilde{n}_{i\sigma}\tilde{n}_{j\bar{\sigma}}) = -2(S_i \cdot S_j - \frac{1}{4}\tilde{n}_i\tilde{n}_j). \qquad (53)$$

Writing the two terms in the summation over $\sigma'$ explicitly, the zeroth-order effective Hamiltonian Eq (29) can be written as

$$H_s^{(0)} = P_s H P_s = T_0 \sum_{i,\sigma} \tilde{n}_{i\sigma} \\ + \sum_{ij\sigma} \gamma_{ij} \left( -t \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \frac{1}{2} Y (c_{i\sigma}^{\dagger} c_{j\sigma} c_{j\bar{\sigma}}^{\dagger} c_{i\bar{\sigma}} + \tilde{n}_{i\sigma}\tilde{n}_{j\bar{\sigma}}) \right. \\ \left. + \frac{1}{2} (V - Y) (\tilde{n}_{i\sigma}\tilde{n}_{j\sigma} + \tilde{n}_{i\sigma}\tilde{n}_{j\bar{\sigma}}) \right). \qquad (54)$$
Therefore, using Eqs (52) and (53), Eq (54) takes the form

$$H_s^{(0)} = T_0 \sum_{i,\sigma} \tilde{n}_{i\sigma} - t \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \sum_{ij} \gamma_{ij} \left( -Y(S_i \cdot S_j - \frac{1}{4} \tilde{n}_i \tilde{n}_j) + \frac{1}{2}(V-Y) \tilde{n}_i \tilde{n}_j \right). \quad (55)$$

We point out again that the terms parametrized by Y and $(V-Y)$ in Eq (55) are the neglected terms in the t-J model. Similarly, for the first-order term $(H_s^{(1)})_{2site}$ given in Eq (42), we have

$$(H_s^{(1)})_{2site} = 2J_0 \sum_{ij} \gamma_{ij} (S_i \cdot S_j - \frac{1}{4} \tilde{n}_i \tilde{n}_j). \quad (56)$$

Summing up the terms in Eqs (55) and (56), the effective Hamiltonian can be expressed as in the t-J model form

$$H_s = T_0 \sum_{i\sigma} \tilde{n}_{i\sigma} - t \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \sum_{ij} \gamma_{ij} \left( (2J_0 - Y)S_i \cdot S_j - \frac{1}{4}(2J_0 + Y - 2V)\tilde{n}_i \tilde{n}_j \right). \quad (57)$$

The above equation shows that when $V = Y$, the model reduces to the t-J model with $J = 2J_0 - V$. However, it has been argued that $V \gg Y$ in the usual cases [3]. Also, since the value of $V$ is usually large, the condition of $J = 2J_0 - V > 0$ (necessary for antiferromagnetic and superconducting states) is difficult to be realized. Therefore, the condition $V = Y$ may be rarely physically satisfied. Finally, using the parameters $p$ and $q$ defined in Eq (46), the model can be written as

$$H_s = T_0 \sum_{i\sigma} \tilde{n}_{i\sigma} - t \sum_{ij,\sigma} \gamma_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \sum_{ij} \gamma_{ij} \left[ J S_i \cdot S_j + \frac{1}{4} J' \tilde{n}_i \tilde{n}_j \right], \quad (58)$$

where $J = p - q$ and $J' = p + q$. When $J' = -J$, which means $p = 0$ (i.e. the condition $V = Y$), one has the t-J model. Our Green function analysis shows that a reasonable superconductivity dome appears at a finite value of $p$ ($p/q \sim 2.5$), which again seems to indicate that the coefficient of $\tilde{n}_i \tilde{n}_j$ of the conventional t-J model, $-J/4$, is not in the physical range.

Eqs (45) and (58) are two alternative expressions of the new model, which reveal the origin of the spin correlations from two different though equivalent perspectives. For two nonempty nearest-neighbor sites, there are two different spin configurations: the spin directions are either the same or the opposite. Eq (45) indicates that the energies for the two spin configurations are, in general, different ($p \neq q$), and that the relative direction of the spins on the nearest-neighbor sites is determined by the competition between the energies of the two spin configurations. On the other hand, as shown in Eq (58), if the interaction terms are rearranged such that the amount of energy $JS_i^z S_j^z$ from the p- and
$q$-terms is combined into the $J$-term to form a complete Heisenberg interaction term [17], the energies of the two spin configurations become the same, namely, $J'n_i n_j/4$. Then the spin correlations can be attributed to the Heisenberg interaction.

# 4 Conclusion

Starting with a generalized Hubbard model with the next leading order Coulomb interaction terms, we have systematically developed a new model for the strongly correlated electrons in solids using the Gutzwiller projection approach. Our study demonstrates that the origin of the Heisenberg interaction term can be attributed to the competition between the energies of the two spin configurations of the nearest-neighbor electrons. A comparison of the conventional $t-J$ model and the new model has been discussed. It seems that the coefficient of one of the interaction terms in the $t-J$ model, i.e. $\tilde{n}_i\tilde{n}_j$, is not in the physical range, which may be the reason why the superconductivity dome has not been obtained from the $t-J$ model.

We have studied the magnetism and superconductivity of this new model. The main results of this research have been summarized in a recent paper [18]. Our study indicates that this model provides a more complete description of the physics of strongly correlated electron systems. The system is not necessarily in a ferromagnetic state as temperature $T \to 0$ at any doping level $\delta \ge 0$. The system, however, must be in an antiferromagnetic state at the origin of the doping-temperature ($\delta-T$) plane ($T \to 0$, $\delta = 0$). Moreover, the model exhibits superconductivity in a doped region at sufficiently low temperatures. The details of the magnetic and superconducting properties of this model will be presented in forthcoming papers.

# References

[1] P. W. Anderson, *Science* **235** (1987) 1196.

[2] J. G. Bednorz and K. A. Müller, *Z. Phys. B* **64**, 189 (1986).

[3] J. Hubbard, *Proc. R. Soc. London A* **276**, 238 (1963), *Proc. R. Soc. London A* **281**, 401 (1964).

[4] M. C. Gutzwiller, *Phys. Rev. Lett.* **10**, 159 (1963).

[5] J. Spalek, *t-J model then and now: A personal perspective from the pioneering times*. (https://arxiv.org/pdf/0706.4236v1.pdf).

[6] F. C. Zhang and T. M. Rice, *Phys. Rev. B* **37**, 3759 (1988).

[7] Gerald D. Mahan, *Many-Particle Physics* (Plenum Press, New York, 1990).
[8] Z. Z. Li, *Gu Ti Li Lun [Solid State Theory], in Chinese* (Higher Education Press, Beijing, 1985)

[9] From now on, we will use the term *two-site interaction terms* to mean 1) the interaction term is a two-center integral and 2) the two sites involved are the nearest neighbors.

[10] R. Strack and D. Vollhardt, *Phys. Rev. Lett.* **70**, 2637 (1993)

[11] T. Sun, *Phil. Mag. Lett.* **73**, 201 (1996)

[12] N. F. Mott, *Proc. Phys. Soc. Lond.*, **A 62**, 416 (1949); *Metal-Insulator Transitions* (London: Taylor and Francis, 1990).

[13] D. Belitz and T. R. Kirkpatrick, *Rev. Mod. Phys.*, **66**, 261 (1994).

[14] In Ref [11], we have indicated that the X-term alone can trigger the first-order metal-insulator phase transition.

[15] D. Vollhardt, in *Proceedings of the International School of Physics "Enrico Fermi" Course CXXI*, edited by R. A. Broglia and J. R. Schrieffer (North Holland, Amsterdam, 1994), p. 31.

[16] I. Hubac and S. Wilson *Brillouin-Wigner Methods for Many-Body Systems* (Springer Netherlands, 2010)

[17] A direct way to go from Eq (45) to Eq (58) is by rearranging the interaction terms with the following: 1) subtract $JS_i^zS_j^z$ from the p- and q-terms, and 2) add it to the J-term.

[18] T. Sun, *Green Function Theory of Strongly Correlated Electron Systems*, 2017 (https://arxiv.org/abs/1703.08424).