# Green Function Theory of Strongly Correlated Electron Systems

Tao Sun

► To cite this version:

Tao Sun. Green Function Theory of Strongly Correlated Electron Systems. 2017. [hal-01491501]

HAL Id: hal-01491501

https://hal.science/hal-01491501

Preprint submitted on 17 Mar 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
Green Function Theory of Strongly Correlated
Electron Systems

Tao Sun *

March 17, 2017

Abstract

A novel effective Hamiltonian in the subspace of singly occupied
states is obtained by applying the Gutzwiller projection approach to
a generalized Hubbard model with the interactions between two nearest-
neighbor sites. This model provides a more complete description of the
physics of strongly correlated electron systems. The system is not nec-
essarily 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. We summarize the studies and provide a
phase diagram of the antiferromagnetism and the superconductivity of
the model in the $\delta-T$ plane here. Details will be presented in subsequent
papers.

# 1 Introduction

Since the discovery of high-temperature superconductivity in 1986 [1], a tremen-
dous number of studies have been devoted to the understanding of its under-
lying mechanism. It is generally believed that some essential physics of this
strongly correlated electron system can be described by the Hubbard model
[2, 3] and its strong-coupling limit, the *t*-*J* model [4, 5]. The Hubbard model
is the simplest approximation of the general Hamiltonian of the interacting
electron systems, in which all the Coulomb interaction terms except the on-
site term (*U*-term) are neglected. It is possible that the neglected interaction

<footnote>*Contact: tao@taosphysics.net</footnote>
terms may play important roles in the understanding of the physics of the strongly correlated electron systems. That is, some important physics might have been missed in the Hubbard model. As a result, the t-J model, which is a projection of the Hubbard model to the subspace of singly occupied states, may not be sufficient to describe the essential physics of the strong coupling systems. In this study, we will apply the Gutzwiller projection scheme to a more general model, in which, in addition to the Hubbard on-site term, all of the two-site interaction terms between nearest-neighbor sites are retained in the approximation of the Coulomb potential energy. Our study shows that this novel model provides a more complete understanding of the properties of high-temperature superconductors.

## 2 The Model

The general Hamiltonian describing the dynamics of electrons in Wannier representation can be expressed as [2, 6, 7]

$$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 (1)$$

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_k$

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

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 (3)$$

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 four-center integrals and the corresponding terms in the series of the Coulomb interaction in Eq (1) are the so-called 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 (3) is not negligible only when the sites $i, j, k, l$ are close enough so that the overlaps between the Wannier functions are sufficiently large. Therefore, the Coulomb interaction energy can be conveniently approximated by a number of its leading terms of the series. 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\sigma'}, \quad (4)$$
where $n_{i\sigma} = c_{i\sigma}^{\dagger}c_{i\sigma}$ and $U = \langle ii|\frac{1}{r}|ii\rangle$. If only this on-site term is taken into account, the Hubbard model is obtained. Apparently, the next leading terms are the ones where the set $\{i, j, k, l\}$ actually consists of only one pair of nearest-neighbor sites, which may be referred to as the two-site interaction terms [8]. In this work, we take one step further beyond the Hubbard model. We retain all of the terms up to the two-site interaction terms in the approximation of the Coulomb interaction. Such a Hamiltonian can be written as

$$ 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} \\ + \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\sigma} + n_{j\sigma}) \\ + \frac{1}{2} Y \sum_{ij,\sigma} \gamma_{ij} (c_{i\sigma}^{\dagger} c_{j\sigma} c_{i\sigma}^{\dagger} c_{j\sigma} + \sum_{\sigma'} c_{i\sigma}^{\dagger} c_{j\sigma'} c_{i\sigma'} c_{j\sigma}), \quad (5) $$

where $\gamma_{ij} = 1$ for nearest-neighbor sites $i, j$, and $0$ otherwise, which restricts the summation over the nearest-neighbor pairs, and

$$ \begin{aligned} V &= \langle ij | \frac{1}{r} | ji \rangle, \\ 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, \\ Y &= \langle ii | \frac{1}{r} | jj \rangle = \langle ij | \frac{1}{r} | ji \rangle. \end{aligned} \quad (6) $$

In fact, the $X$ and $V$ terms have been discussed in literature previously [6, 9], and the matrix element $Y$ has been mentioned in [2]. Please note that in Eq (5) only the constant and nearest-neighbor hopping terms of the kinetic energy are considered for consistency.

Using the Gutzwiller projection operator technique, the Hamiltonian (5) can be projected to a subspace of the Hilbert space where only empty and singly occupied sites are allowed [3, 10, 11]. To the leading order of the perturbation expansion, the effective Hamiltonian in the subspace of singly occupied states 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 (7) $$

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\sigma} = \sum_{\sigma} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{i\sigma}$, and the coefficients are given by the following equations

$$ \begin{aligned} p &= V - Y, \\ q &= V - 2J_0, \quad J_0 = (t - X)^2 / U', \quad U' = U + z'(2V - Y), \\ J &= p - q, \end{aligned} \quad (8) $$
where $0 \le z' \le z$. When $V=X=Y=0$, Eq (7) reduces to the $t-J$ model [4]. Rescaling $H_s$ with $t$, this model has only three dimensionless independent material 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 show that this is true. More interestingly, we will show that this model will exhibit superconductivity for some parameter range of $\bar{p} > \bar{q}$.

### 3 Green Function Method

We use the Zubarev Green function technique [2, 12] to study the physics of the model established above. Define the grand canonical Hamiltonian operator [13]

$$K = H_s - \mu N, \qquad (9)$$

where $N$ is the total number operator, $H_s$ is the Hamiltonian given in Eq (7), and $\mu$ is the chemical potential of electrons. The corresponding grand partition function may be written as

$$Z_G = \mathrm{Tr} e^{-\beta K} = e^{-\beta \Omega}, \quad \beta = \frac{1}{kT}, \qquad (10)$$

where $k$ is the Boltzmann's constant and $\Omega$ is the thermodynamic potential of the system. Therefore, the ensemble average of any operator $O$ can be calculated by

$$\langle O \rangle = \mathrm{Tr}[Oe^{\beta(\Omega-K)}]. \qquad (11)$$

Note here that $O$ and $K$ are operators, while $\Omega$ is a $c$-number. The ensemble average $\langle O \rangle$ is assumed to be the measurement value of the observable $O$ in experiments. The retarded Green function of two Fermion operators $A(t)$ and $B(t')$ in the Heisenberg representation can be defined as

$$\ll A(t)|B(t')\gg = -i\theta(t-t')\langle\{A(t), B(t')\}\rangle, \qquad (12)$$

where $\{A, B\} = AB + BA$, the Fermion anticommutation relation, $\theta(t)$ is the usual step function, and $\langle...\rangle$ means the ensemble average with the grand partition function $Z_G$, defined in Eq (11). Since $\ll A(t)|B(t')\gg$ is a function of $t-t'$, it is convenient to define the Fourier transform

$$\ll A|B\gg_{\omega} = \int_{-\infty}^{\infty} d(t-t') \ll A(t)|B(t')\gg e^{i\omega(t-t')}. \qquad (13)$$

It can be shown that the retarded Green function in frequency space satisfies the following equation of motion

$$\omega \ll A|B\gg_{\omega} = \langle\{A, B\}\rangle + \ll[A, K]\ll B\gg_{\omega}, \qquad (14)$$
where again {$A, B$} is the anticommutator of operators $A$ and $B$, and $[A, K]$ is the commutator of the operators $A$ and $K$. The ensemble average $\langle B(t)A(t') \rangle$ can be calculated by

$$ \langle B(t)A(t') \rangle = i \int_{-\infty}^{\infty} \frac{d\omega}{2\pi} \frac{e^{-i\omega(t-t')}}{e^{\beta\omega} + 1} \left[ \langle A|B \rangle_{\omega+i0^+} - \langle A|B \rangle_{\omega-i0^+} \right]. \quad (15) $$

Eqs (14)-(15) are applicable to all subsequent studies of the ferromagnetism, antiferromagnetism, and superconductivity of the model given in Eq (7).

Due to the complexity of the model, we only consider the simplest situation in this study: a three-dimensional square lattice system. The two-dimensional case and other complicated factors are left as future work.

## 4 Ferromagnetism

We first study the ferromagnetism of the model using the standard Zubarev Green function technique. For a uniform system, we assume that the ensemble average of the electron number operator $\langle n_{l\sigma} \rangle$ is independent of lattice site $l$, but dependent on spin $\sigma$. Denoting $n_{\sigma} = \langle n_{l\sigma} \rangle$, the average number of electrons per site and magnetization can be defined as

$$ n = n_{\uparrow} + n_{\downarrow}, \quad (16) $$

$$ m = n_{\uparrow} - n_{\downarrow}, \quad (17) $$

from which, $n_{\sigma}$ can be expressed as

$$ n_{\sigma} = \frac{1}{2}(n + \sigma m), \quad (18) $$

where $\sigma = 1(-1)$ in the expression for spin subscript $\uparrow$ ($\downarrow$). Here $n$ is a known quantity when the doping level is known, from which the chemical potential $\mu$ is determined. For electron-doped materials, $n = 1 + \delta$, while for hole-doped cases, $n = 1 - \delta$, where $\delta$ is the doping concentration. We will focus on the hole-doping only in this study.

The analysis above shows that the central task in the study of the magnetic properties is to calculate the ensemble average of the number operator of electrons, $\langle n_{l\sigma} \rangle$. According to Eq (15), $\langle n_{l\sigma} \rangle$ can be evaluated from the Green function $\langle c_{l\sigma} | c_{l'\sigma}^{\dagger} \rangle$. The retarded Green function in frequency space satisfies the following equation of motion

$$ \omega \ll c_{l\sigma} |c_{l'\sigma}^{\dagger} \rangle \gg \omega = \delta_{ll'} + \langle [c_{l\sigma}, K] | c_{l'\sigma}^{\dagger} \rangle \gg \omega, \quad (19) $$

where $[c_{l\sigma}, K]$ is the commutator of the operators $c_{l\sigma}$ and $K$, and $K$ is given in Eq (9). The Green function $\langle c_{l\sigma} | c_{l'\sigma}^{\dagger} \rangle$ is a function of $\omega$ and $R_{l'} - R_l$ for systems with translational invariance.
Carrying out the commutation relation $[c_{l\sigma}, K]$, to the lowest order of the de-coupling approximation of the sequence of Green function equations, we obtain

$$ \ll c_{l\sigma} |c_{l'\sigma}^{\dagger}| \gg_{\omega} = \frac{1}{N} \sum_{\mathbf{k}} e^{i\mathbf{k} \cdot (\mathbf{R}_l - \mathbf{R}_{l'})} \frac{1}{\omega + \mu - za_{\sigma} - b_{\sigma} \gamma(\mathbf{k})}, \quad (20) $$

where $z$ is the coordination number and the $\sigma$ dependent coefficients $a_\sigma$ and $b_\sigma$ are given by the following equations

$$ \begin{aligned} a_{\sigma} &= a_{1\sigma} + \sigma m a_{2\sigma}, \\ b_{\sigma} &= b_{1\sigma} + \sigma m b_{2\sigma}, \end{aligned} \quad (21) $$

$$ a_{1\sigma} = T_0\delta + t(1+\delta) + (p+q)\delta\left[\frac{1}{4}(1-\delta^2) + \Lambda_{\bar{\sigma}}^2 + \frac{1}{4}m^2\right], $$

$$ a_{2\sigma} = -T_0 - t + \frac{1}{4}(3p-q) - (p+q)\left(\frac{\delta^2}{4} + \Lambda_{\bar{\sigma}}^2 - \frac{1}{4}m^2\right), $$

$$ \begin{aligned} b_{1\sigma} = t &\left[ -\frac{1}{4}(1+\delta)^2 + \Lambda_{\bar{\sigma}}^2 + 2\Lambda_{\bar{\sigma}}\Lambda_{\sigma} - \frac{1}{4}m^2 \right] - (p-q)\Lambda_{\bar{\sigma}} \\ &+ \left[ -\frac{1}{2}(p-q) - 2(p+q)\left(\frac{1}{4}\delta^2 - \Lambda_{\bar{\sigma}}^2 + \frac{1}{4}m^2\right) \right] \Lambda_{\sigma}, \end{aligned} \quad (22) $$

$$ b_{2\sigma} = -\frac{1}{2}t(1+\delta) - (p+q)\delta\Lambda_{\sigma}, $$

and the dispersion relation $\gamma(\mathbf{k})$ is defined as

$$ \gamma(\mathbf{k}) = \sum_{\mathbf{a}} e^{\mathbf{k} \cdot \mathbf{a}}, \quad (23) $$

where $\mathbf{a}$ is the lattice space vector and the summation is over all the nearest-neighbor sites. In Eq (22), $\Lambda_\sigma$ represents the value of the nonzero formal ensemble average $\langle c_{l'\sigma}^{\dagger} c_{l\sigma} \rangle$ for the nearest-neighbor sites $l$ and $l'$. $\Lambda_\sigma$ will be calculated self-consistently [14].

Using Eqs (15) and (20), the ensemble average $\langle c_{l'\sigma}^{\dagger} c_{l\sigma} \rangle$ can be calculated as follows

$$ \langle c_{l'\sigma}^{\dagger} c_{l\sigma} \rangle = \frac{1}{N} \sum_{\mathbf{k}} e^{i\mathbf{k} \cdot (\mathbf{R}_l - \mathbf{R}_{l'})} \frac{1}{e^{\beta(E_{\sigma}(\mathbf{k}) - \mu)} + 1}, \quad (24) $$

where the elementary excitations energy spectrum is given by

$$ E_{\sigma}(\mathbf{k}) = za_{\sigma} + b_{\sigma}\gamma(\mathbf{k}). \quad (25) $$

Eq (24) is the basic formula for the ensemble average $\langle c_{l'\sigma}^{\dagger} c_{l\sigma} \rangle$, from which we are able to obtain the self-consistent equations of $n_\sigma$ and $\Lambda_\sigma$

$$ n_{\sigma} = z \int_{-1}^{1} dx \rho(x) \frac{1}{e^{\beta z b_{\sigma}(x-x_{0\sigma})} + 1}, \quad (26) $$

$$ \Lambda_{\sigma} = \frac{z^2}{z_1} \int_{-1}^{1} dx \rho(x) \frac{x}{e^{\beta z b_{\sigma}(x-x_{0\sigma})} + 1}, \quad (27) $$
where $x_{0\sigma} = (\bar{\mu} - a_{\sigma})/b_{\sigma}$ with $\bar{\mu} = \mu/z$, and the density of state can be expressed as

$$ \rho(x) = \frac{1}{3\pi} \sqrt{1 - x^2}, \qquad (28) $$

for three dimensional systems [2]. Here the summation over **k** has been converted to an integral over $x \to \gamma(\mathbf{k})/z$ and $z_1 = z(1 - \delta)$, the average number of the nonempty nearest-neighbor sites of the site $l$ at doping level $\delta$.

Eqs (26) and (27) are the self-consistent equations describing the ferromagnetism of the model described by Eq (7). Five unknowns, $n_{\uparrow}$, $n_{\downarrow}$, $\Lambda_{\uparrow}$, $\Lambda_{\downarrow}$, and $\mu$ need to be determined from the system of equations. When Eqs (26) and (27) are written down explicitly for each spin, one has four independent equations. The fifth equation is given in Eq (16), which indicates that there is only one independent unknown in the set of $n_{\uparrow}$ and $n_{\downarrow}$. Given $\delta$ (or $n = 1 - \delta$) and $T$ (or $\beta = 1/kT$), solving Eqs (16), (26), and (27) for $n_{\sigma}$, $\Lambda_{\sigma}$, and $\mu$, we can obtain the information about the magnetization $m$ as a function of $\delta$ and $T$. Of course, $m$ depends on the material parameters $T_0 = T_0/t$, $\bar{p} = p/t$, and $\bar{q} = q/t$ as well.

The system can be solved numerically. We briefly outline the main features of the ferromagnetism as follows:

1. An increase in $\bar{q}$ or decrease in $\bar{p}$ is favorable for the magnetization $m$.

2. The magnetization $m$ decreases with temperature.

3. The doping dependence of $m$ is complicated. The favorability of doping for $m$ depends on the values of $\bar{p}$ and $\bar{q}$.

The detailed results will be presented in a subsequent paper. Our main purpose of discussing the ferromagnetism here is to examine the singularity of the self-consistent equations and make a comparison with that of the antiferromagnetism and superconductivity.

The most interesting feature of the system of equations (26) and (27) is that, when $m = 0$, there is no singularity as $T \to 0$ (or $\beta \to \infty$). This feature persists even for the situation of $\delta = 0$ and $T \to 0$. When $\beta \to \infty$, the integrals in Eq (26) and (27) can be carried out analytically with the following results: For $b_{\sigma} > 0$, integrating gives

$$ \begin{aligned} n_{\sigma} &= \frac{1}{2} + \frac{1}{\pi} \left[ x_{0\sigma} \sqrt{1 - x_{0\sigma}^2} + \arcsin(x_{0\sigma}) \right], \\ \Lambda_{\sigma} &= -\frac{4}{\pi z_1} (1 - x_{0\sigma}^2)^{3/2}, \end{aligned} \qquad (29) $$
while for $b_σ < 0$, we have

$$ 
\begin{aligned}
n_{\sigma} &= \frac{1}{2} - \frac{1}{\pi} \left[ x_{0\sigma} \sqrt{1 - x_{0\sigma}^2} + \arcsin(x_{0\sigma}) \right], \\
\Lambda_{\sigma} &= \frac{4}{\pi z_1} (1 - x_{0\sigma}^2)^{3/2},
\end{aligned} 
\qquad (30) $$

Indeed, Eqs (29) and (30) are not singular for any value of $m$ even for the case of $\delta = 0$. This means that the model of Eq (7) does not necessarily have a ferromagnetic phase with nonzero $m$ at the condition of $T \to 0$. That is, $m = 0$ is a possible solution of the model at zero temperature. As we will see below, this is in remarkable contrast to the cases of the antiferromagnetism and superconductivity of this model.

# 5 Antiferromagnetism

We now proceed to study the antiferromagnetism of the model. As we will see below, this model exhibits much more interesting antiferromagnetic features than those of the ferromagnetism. The most striking one is that the system is necessarily in an antiferromagnetic state as $T \to 0$ at $\delta = 0$.

The formalism is similar to the discussion of the ferromagnetism. It is convenient, however, to formally divide the original lattice into two sublattices so that the nearest-neighbor sites belong to different sublattices. We use $i_1, j_1, ..., $ to label the sites in the first sublattice, and $i_2, j_2, ..., $ to label the sites in the second sublattice. The central task is to calculate the average number of electrons per site for each spin direction in both sublattices. Denoting the mean number of electrons with spin $\sigma$ in the sublattices as $n_{1\sigma} = \langle n_{i_1\sigma} \rangle$ and $n_{2\sigma} = \langle n_{i_2\sigma} \rangle$ respectively, we can define the mean number of electrons per site and magnetization for the first sublattice

$$ n_1 = n_{1\uparrow} + n_{1\downarrow}, \qquad (31) $$

$$ m_1 = n_{1\uparrow} - n_{1\downarrow}. \qquad (32) $$

Similarly, for the second sublattice, we have

$$ n_2 = n_{2\uparrow} + n_{2\downarrow}, \qquad (33) $$

$$ m_2 = n_{2\uparrow} - n_{2\downarrow}. \qquad (34) $$

Eqs (31)-(34) immediately lead to the following expressions for each sublattice

$$ n_{1\sigma} = \frac{1}{2}(n_1 + \sigma m_1), \qquad (35) $$

$$ n_{2\sigma} = \frac{1}{2}(n_2 + \sigma m_2). \qquad (36) $$
We make two important assumptions: 1) $n_1 = n_2 = n$, and 2) $m_1 = -m_2 = m$. Therefore, Eqs (35) and (36) become

$$n_{1\sigma} = \frac{1}{2}(n + \sigma m), \qquad (37)$$

$$n_{2\sigma} = \frac{1}{2}(n - \sigma m), \qquad (38)$$

where $n = 1 - \delta$ and $m$ will be calculated self-consistently. The above analysis shows that we only need to find the average electron number in one sublattice. For example, if $n_{1\sigma}$ is known, then $m = n_{1\uparrow} - n_{1\downarrow}$ is known, and therefor $n_{2\sigma}$ is known.

Similar to the ferromagnism case, the ensemble averages $n_{1\sigma} = \langle n_{i_1\sigma} \rangle$ can be calculated from the Green function $\ll c_{l_1\sigma}|c_{l'_1\sigma}^\dagger \gg$. Nevertheless, further analysis shows that the two Green functions $\ll c_{l_1\sigma}|c_{l'_1\sigma}^\dagger \gg$ and $\ll c_{l_2\sigma}|c_{l'_1\sigma}^\dagger \gg$ are coupled to each other even at the lowest order of the decoupling approximation of the Green function equations. Therefore we need to discuss the two Green functions simultaneously. In frequency space, the two Green function satisfy the following equations of motion

$$\omega \ll c_{l_1\sigma}|c_{l'_1\sigma}^\dagger \gg_{\omega} = \delta_{ll'} + \ll [c_{l_1\sigma}, K] |c_{l'_1\sigma}^\dagger \gg_{\omega}, \qquad (39)$$

$$\omega \ll c_{l_2\sigma}|c_{l'_1\sigma}^\dagger \gg_{\omega} = \ll [c_{l_2\sigma}, K] |c_{l'_1\sigma}^\dagger \gg_{\omega}, \qquad (40)$$

where the anticommutation relations $\{c_{l_1\sigma}, c_{l'_1}^\dagger\} = \delta_{ll'}$ and $\{c_{l_2\sigma}, c_{l'_1}^\dagger\} = 0$ have been used. Using Eq (9), the commutation relations $[c_{l_1\sigma}, K]$ and $[c_{l_2\sigma}, K]$ can be obtained exactly. Nonetheless, an approximation must be used to break off the chain of Green function equations. To the lowest order of the decoupling approximation, the Green functions can expressed as

$$\ll c_{l_1\sigma}|c_{l'_1\sigma}^\dagger \gg_{\omega} = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} e^{i\mathbf{k} \cdot (\mathbf{R}_{l_1} - \mathbf{R}_{l'_1})} g_{\mathbf{k}\sigma}(\omega), \qquad (41)$$

$$\ll c_{l_2\sigma}|c_{l'_1\sigma}^\dagger \gg_{\omega} = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} e^{i\mathbf{k} \cdot (\mathbf{R}_{l_2} - \mathbf{R}_{l'_1})} f_{\mathbf{k}\sigma}(\omega), \qquad (42)$$

where

$$g_{\mathbf{k}\sigma}(\omega) = \frac{A_{\mathbf{k}\sigma}}{\omega + \mu - E_1(\mathbf{k})} + \frac{B_{\mathbf{k}\sigma}}{\omega + \mu - E_2(\mathbf{k})}, \qquad (43)$$

$$f_{\mathbf{k}\sigma}(\omega) = \frac{C_{\mathbf{k}}}{\omega + \mu - E_1(\mathbf{k})} + \frac{D_{\mathbf{k}}}{\omega + \mu - E_2(\mathbf{k})}, \qquad (44)$$

with the elementary excitation spectrum

$$E_1(\mathbf{k}) = za - \sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}, \qquad (45)$$

$$E_2(\mathbf{k}) = za + \sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}, \qquad (46)$$
and the partial-fraction decomposition coefficients

$$A_{k\sigma} = \frac{1}{2} \left[ 1 - \frac{\sigma zma'}{\sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}} \right], \quad (47)$$

$$B_{k\sigma} = \frac{1}{2} \left[ 1 + \frac{\sigma zma'}{\sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}} \right], \quad (48)$$

$$C_{\mathbf{k}} = - \frac{b'\gamma(\mathbf{k})}{2\sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}} = D_{\mathbf{k}}, \quad (49)$$

where $\gamma(\mathbf{k})$ is given in Eq (23), $b' = b(1 - \delta)$, and the coefficients $a, a',$ and $b$ are given the following equations

$$a = \frac{T_0}{z} \delta + t(1 + \delta)\Lambda + (p + q)\delta \left[ \frac{1}{4}(1 - \delta^2) + \Lambda^2 + \frac{1}{4}m^2 \right], \quad (50)$$

$$a' = -\frac{T_0}{z} + t\Lambda + \frac{1}{4}(3q - p) + (p + q)\left[-\frac{1}{4}\delta^2 + \Lambda^2 + \frac{1}{4}m^2\right], \quad (51)$$

$$b = t \left[ -\frac{1}{4}(1 + \delta)^2 + 3\Lambda^2 + \frac{1}{4}m^2 \right] + \left[ -J - \frac{1}{2}(p - q) + 2(p + q)\left(-\frac{1}{4}\delta^2 + \Lambda^2 + \frac{1}{4}m^2\right) \right] \Lambda. \quad (52)$$

Here $\Lambda$ is the value of the nonzero formal ensemble average $\langle c_{l_1\sigma}^{\dagger} c_{l_2\sigma} \rangle$, where $l_1$ and $l_2$ are the nearest-neighbor sites. We have assumed that $\Lambda$ is independent of spin $\sigma$. Similar to the ferromagnetism case, $\Lambda$ will be calculated self-consistently [14]. Please note that $a, a'$, and $b$ are all independent of $\sigma$ explicitly. As can be seen from Eqs (45) and (46), two branches of the elementary excitation spectrum exist in the system, both of which are spin $\sigma$ independent. The coefficients $A_{k\sigma}$ and $B_{k\sigma}$ depend on spin $\sigma$ explicitly, while $C_k$ and $D_k$ are spin independent. Therefore $f_{k\sigma}(\omega)$ does not depend on $\sigma$ explicitly, which results from the assumption that $\Lambda$ is independent of $\sigma$.

Using Eqs (15), (41) and (42), we can obtain the following ensemble averages

$$\langle c_{l_1'\sigma}^{\dagger} c_{l_1\sigma} \rangle = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} e^{i\mathbf{k}\cdot(\mathbf{l}_1 - \mathbf{l}'_1)} \left[ \frac{A_{k\sigma}}{e^{\beta(E_1(\mathbf{k}) - \mu)} + 1} + \frac{B_{k\sigma}}{e^{\beta(E_2(\mathbf{k}) - \mu)} + 1} \right], \quad (53)$$

$$\langle c_{l_1'\sigma}^{\dagger} c_{l_2\sigma} \rangle = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} e^{i\mathbf{k}\cdot(\mathbf{l}_2 - \mathbf{l}'_1)} \left[ \frac{C_{\mathbf{k}}}{e^{\beta(E_1(\mathbf{k}) - \mu)} + 1} + \frac{D_{\mathbf{k}}}{e^{\beta(E_2(\mathbf{k}) - \mu)} + 1} \right], \quad (54)$$

where we have used $\mathbf{l}_1 - \mathbf{l}'_1$ to represent $\mathbf{R}_{\mathbf{l}_1} - \mathbf{R}_{\mathbf{l}'_1}$ for conciseness. Setting $\mathbf{l}_1 = \mathbf{l}'_1$ in Eq (53), we have

$$n_{1\sigma} = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} \left[ \frac{A_{k\sigma}}{e^{\beta(E_1(\mathbf{k}) - \mu)} + 1} + \frac{B_{k\sigma}}{e^{\beta(E_2(\mathbf{k}) - \mu)} + 1} \right]. \quad (55)$$
Noticing $n = n_{1\uparrow} + n_{1\downarrow}$ and $m = n_{1\uparrow} - n_{1\downarrow}$, and using Eqs (47)-(48), Eq (55) immediately leads to the following self-consistent equations for $\mu$ and $m$

$$n = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} \left[ \frac{1}{e^{\beta(E_1(\mathbf{k})-\mu)} + 1} + \frac{1}{e^{\beta(E_2(\mathbf{k})-\mu)} + 1} \right], \quad (56)$$

$$1 = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} \frac{-za'}{\sqrt{(zma')^2 + (b\gamma(\mathbf{k}))^2}} \left[ \frac{1}{e^{\beta(E_1(\mathbf{k})-\mu)} + 1} - \frac{1}{e^{\beta(E_2(\mathbf{k})-\mu)} + 1} \right]. \quad (57)$$

The self-consistent equation for $\Lambda$ can be obtained from Eq (54) by setting $l_2$ and $l'_1$ as the nearest neighbors

$$z_1 \Lambda = \frac{2}{N} \sum_{\mathbf{k}}^{N/2} \frac{-b'\gamma(\mathbf{k})^2}{2\sqrt{(zma')^2 + (b'\gamma(\mathbf{k}))^2}} \left[ \frac{1}{e^{\beta(E_1(\mathbf{k})-\mu)} + 1} - \frac{1}{e^{\beta(E_2(\mathbf{k})-\mu)} + 1} \right], \quad (58)$$

where Eq (49) has been used. Here $z_1 = z(1 - \delta)$, the average number of the nonempty nearest neighbors of the site $l_2$. Eqs (45) and (46) show that $E_1(\mathbf{k}) \le E_2(\mathbf{k})$. Thus from Eqs (57) and (58), we see that for meaningful solutions, 1) $a'$ must be negative, and 2) $\Lambda$ and $b$ must have the opposite sign.

Finally, introducing a new variable $x \to \gamma(\mathbf{k})/z$, and noticing that the integrand depends on $x$ through $x^2$, for the case of $b > 0$, the self-consistent equations can be written in the following integration form

$$n = 2z \int_0^1 dx \rho(x) \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} + \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \quad (59)$$

$$\Lambda = -\frac{z^2}{z_1} \int_0^1 dx \rho(x) \frac{x^2}{\sqrt{\bar{m}^2 + x^2}} \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} - \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \quad (60)$$

$$b' = -2za' \int_0^1 dx \rho(x) \frac{1}{\sqrt{\bar{m}^2 + x^2}} \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} - \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \quad (61)$$

where the density of states $\rho(x)$ is given in Eq (28) and the two branches of energy excitation spectrum take the form

$$\epsilon_1(x) = zb' \left[ -x_0 - \sqrt{\bar{m}^2 + x^2} \right], \quad (62)$$

$$\epsilon_2(x) = zb' \left[ -x_0 + \sqrt{\bar{m}^2 + x^2} \right], \quad (63)$$

where $x_0 = (\mu/z - a)/b'$ and $\bar{m}^2 = (ma'/b')^2$. Eqs (59)-(61) are the basic self-consistent equations for antiferromagnetism of the model in Eq (7). At a given condition of the doping level $\delta$ and temperature $T$, three unknowns,
$x_0$, $\Lambda$, and $m$ can be determined from the system of equations. Of course, $x_0$, $\Lambda$, and $m$ also depend on the dimensionless material parameters $\bar{T}_0 = T_0/t$, $\bar{p} = p/t$, $\bar{q} = q/t$. In the following, we discuss the case of $T_0 = 0$.

The system is in ground state at absolute zero temperature. This is the $\delta$-axis in the $\delta-T$ plane. When $T \to 0$, the integral in Eq (59) can be carried out analytically with the result

$$ \delta = \frac{2}{\pi} [\xi \sqrt{1 - \xi^2} + \arcsin(\xi)], \qquad (64) $$

where $\xi = \sqrt{x_0^2 - \bar{m}^2}$. In the case of $b > 0$, $x_0 < 0$ is applicable to the case of $n = 1 - \delta$ (hole doping), while $x_0 > 0$ is applicable to the case of $n = 1 + \delta$ (electron doping). When $T \to 0$, for the hole doping case ($x_0 < 0$), Eqs (60) and (61) can be written as

$$ \Lambda = -\frac{z^2}{z_1} \int_{\xi}^{1} dx \rho(x) \frac{x^2}{\sqrt{\bar{m}^2 + x^2}}, \qquad (65) $$

$$ b' = -2za' \int_{\xi}^{1} dx \rho(x) \frac{1}{\sqrt{\bar{m}^2 + x^2}}, \qquad (66) $$

where $\xi$ is given by Eq (64). There is no singularity as long as either $\xi$ or $\bar{m}$ is nonzero. The quantity $\xi$ depends on $\delta$ only and Eq (64) is not coupled with Eqs (65) and (66). Thus $\xi$ can be solved separately from Eq (64) for a given $\delta$. When $\delta$ is small, to the leading order, Eq (64) can be solved with the result

$$ \xi \approx \frac{\pi}{4}\delta, \qquad (67) $$

which leads to an approximate relation between $x_0$ and $\bar{m}$ on the $\delta$-axis ($T \to 0$).

$$ x_0^2 \approx \left(\frac{\pi}{4}\delta\right)^2 + \bar{m}^2. \qquad (68) $$

A more complete solution of the system of equations (64)-(66) can be obtained numerically for given $\delta$, $\bar{p}$, and $\bar{q}$. Now we discuss two special points on the $\delta$-axis ($T \to 0$).

An interesting point on the $\delta$-axis is where $m$ vanishes. At this point, $\xi = |x_0|$ and the integrals in the self-consistent equations (65)-(66) can be carried out analytically. Use $\delta_0$ to denote the doping concentration at this point, which can be determined by the following equations

$$ \delta_0 = -\frac{2}{\pi} \left[ x_0 \sqrt{1 - x_0^2} + \arcsin(x_0) \right], \qquad (69) $$

$$ \Lambda = -\frac{4}{\pi z_1} (1 - x_0^2)^{3/2}, \qquad (70) $$

$$ b' = \frac{4}{\pi} a' \left[ \sqrt{1 - x_0^2} + \ln|x_0| - \ln\left(1 + \sqrt{1 - x_0^2}\right) \right]. \qquad (71) $$
To the leading order of $\delta_0$, $x_0$ can be solved from Eq (69) as

$$x_0 \approx -\frac{\pi}{4}\delta_0. \qquad (72)$$

The term $\ln|x_0|$ in Eq (71) shows that a meaningful $\delta_0$ must be finite. Solving $x_0$ and $\Lambda$ from Eqs (69) and (70), and then substituting their values into Eq (71), we obtain a single variable equation for $\delta_0$. Given $\bar{p}$ and $\bar{q}$, this equation can be solved to obtain $\delta_0(\bar{p}, \bar{q})$. In practice, the value of $\delta_0$ is known from experiments, so that $\delta_0(\bar{p}, \bar{q})$ specifies a relation between the dimensionless parameters $\bar{p}$ and $\bar{q}$. For example, in the case of $b > 0$, for $\delta_0 = 0.05$, $\bar{q} = 0.63$ corresponds to $r = \bar{p}/\bar{q} \approx 2.5$, while $\bar{q} = 1$ corresponds to $r = \bar{p}/\bar{q} \approx 3.73$.

The most interesting case is the ground state at half-filling. This is the origin of the $\delta-T$ plane. Eq (64) indicates that $\xi = 0$ and therefore $\tilde{x}_0^2 = \bar{m}^2$ at this point. Noticing that at $\delta = 0$, $z_1 = z$, and $b' = b$, the self-consistent equations become

$$\Lambda = -z \int_0^1 dx \rho(x) \frac{x^2}{\sqrt{\bar{m}^2 + x^2}}, \qquad (73)$$

$$b = -2za' \int_0^1 dx \rho(x) \frac{1}{\sqrt{\bar{m}^2 + x^2}}. \qquad (74)$$

Apparently, the integral in Eq (74) is divergent when $m = 0$, which indicates that $m = 0$ is not a possible solution of the system as $T \to 0$ at $\delta = 0$. That is, the system must be in an antiferromagnetic state with a nonzero $m$ at the origin of the $\delta-T$ plane. This is one of the most striking features of the model described by Eq (7). Since in this case $T \to 0$ and $\delta = 0$, the magnetization $m$ depends on the two parameters $\bar{q}$ and $\bar{p}$ (or $\bar{q}$ and $r = \bar{p}/\bar{q} = p/q$) only. Our results show that for a given value of $\bar{q}$, a nonzero $m$ starts with a threshold value $r \ge r_0$. Then $m$ increases with $r$, and finally saturates to a maximum value of $m_{max} \approx 0.6698$ when $r$ is sufficiently large. Generally, the greater the $\bar{q}$, the greater the threshold value $r_0$. The threshold value of $r$ for a given $\bar{q}$ is related to the condition of $a' < 0$ and $b' > 0$.

In the simplest case where $m = 0$, for the case of $b > 0$, the self-consistent equations become

$$n = 2z \int_0^1 dx \rho(x) \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} + \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \qquad (75)$$

$$\Lambda = -\frac{z^2}{z_1} \int_0^1 dx \rho(x) x \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} - \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \qquad (76)$$

$$b' = -2za' \int_0^1 dx \rho(x) \frac{1}{x} \left[ \frac{1}{e^{\beta\epsilon_1(x)} + 1} - \frac{1}{e^{\beta\epsilon_1(x)} + 1} \right], \qquad (77)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Antiferromagnetic phase boundary lines for typical values of $\bar{q} = q/t$ and $r$. The region below each line is the antiferromagnetic state.</figcaption></figure>

where

$$ \epsilon_1(x) = zb'(-x_0 - x), \qquad (78) $$

$$ \epsilon_2(x) = zb'(-x_0 + x). \qquad (79) $$

Now there are only two unknowns $x_0$ and $\Lambda$, which can be determined from Eqs (75) and (76) when $\delta$ and $T$ are given. The third equation (77) will give the Néel temperature at a doping concentration $\delta$, $T_N(\delta)$, in the $\delta$-$T$ plane, which is the boundary of the area where $m \neq 0$. Fig. 1 shows the antiferromagnetic phase diagrams for some typical values of $\bar{q}$ and $r$. As discussed above, Eq (77) is singular at the point of $T \to 0$ and $\delta = 0$. Therefore, the boundary of $m \neq 0$ must bypass the origin of the $\delta$-$T$ plane.
# 6 Superconductivity

In order to study the superconductivity of the model described by Eq (7), one needs to calculate the ensemble average $\langle c_{l\downarrow}c_{l'\uparrow} \rangle$ or its complex conjugate $\langle c_{l'\uparrow}^{\dagger}c_{l\downarrow}^{\dagger} \rangle$. These quantities provide the information about the state of the Cooper pairs, from which the properties of the superconducting and pseudogap state can be obtained. The ensemble average $\langle c_{l'\uparrow}^{\dagger}c_{l\downarrow}^{\dagger} \rangle$ can be calculated from the Green function $\ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg$. Further analysis shows that the Green function $\ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg$ is coupled with $\ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg$ even at the lowest order of decoupling approximation. Therefore, we need to deal with both $\ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg$ and $\ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg$ simultaneously in the study of the superconductivity. In $\omega$-space, the equations of motion for both Green functions $\ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg_{\omega}$ and $\ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg_{\omega}$ are

$$ \omega \ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg_{\omega} = \delta_{ll'} + \ll [c_{l\uparrow}, K] |c_{l'\uparrow}^{\dagger} \gg_{\omega}, \quad (80) $$

$$ \omega \ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg_{\omega} = \ll [c_{l\downarrow}^{\dagger}, K] |c_{l'\uparrow}^{\dagger} \gg_{\omega}, \quad (81) $$

where the Fermion anticommutation relations $\{c_{l\uparrow}, c_{l'\uparrow}^{\dagger}\} = \delta_{ll'}$ and $\{c_{l\downarrow}^{\dagger}, c_{l'\uparrow}^{\dagger}\} = 0$ have been taken into account. For translational invariant systems, $\ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg_{\omega}$ and $\ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg_{\omega}$ are functions of frequency $\omega$ and $R_l - R_{l'}$.

To the lowest order of decoupling approximation, we obtain the Green functions as follows

$$ \ll c_{l\uparrow}|c_{l'\uparrow}^{\dagger} \gg_{\omega} = \frac{1}{N} \sum_{\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{R}_l - \mathbf{R}_{l'})} \left[ \frac{A_{\mathbf{k}}}{\omega - E_{\mathbf{k}}} + \frac{B_{\mathbf{k}}}{\omega + E_{\mathbf{k}}} \right], \quad (82) $$

$$ \ll c_{l\downarrow}^{\dagger}|c_{l'\uparrow}^{\dagger} \gg_{\omega} = \frac{1}{N} \sum_{\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{R}_l - \mathbf{R}_{l'})} \left[ \frac{C_{\mathbf{k}}}{\omega - E_{\mathbf{k}}} + \frac{D_{\mathbf{k}}}{\omega + E_{\mathbf{k}}} \right], \quad (83) $$

where the partial-fraction decomposition coefficients are given by

$$ A_{\mathbf{k}} = \frac{1}{2} \left( 1 + \frac{-\mu + za + b'\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \right), \quad (84) $$

$$ B_{\mathbf{k}} = \frac{1}{2} \left( 1 - \frac{-\mu + za + b'\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \right), \quad (85) $$

$$ C_{\mathbf{k}} = \frac{1}{2} \frac{(zc + d'\gamma_{\mathbf{k}})\Delta^*}{E_{\mathbf{k}}} = -D_{\mathbf{k}}, \quad (86) $$

and the elementary excitation energy spectrum takes the form

$$ E_{\mathbf{k}} = \sqrt{(\mu - za - b'\gamma_{\mathbf{k}})^2 + (zc + d'\gamma_{\mathbf{k}})^2 |\Delta|^2}. \quad (87) $$
Here the dispersion relation $\gamma_k$ is given in Eq (23) and $b' = b(1 - \delta)/(1 + r)$, $d' = d(1 - \delta)r/(1 + r)$, where $r = p/q$, and the coefficients $a, b, c$, and $d$ are expressions of $\Lambda$ and $\Delta$

$$a = t(1 + \delta)\Lambda + (p + q)\delta\left[\frac{1}{4}(1 - \delta^2) + \Lambda^2 - |\Delta|^2\right], \qquad (88)$$

$$b = t\left[-\frac{1}{4}(1 + \delta)^2 + 3\Lambda^2 + |\Delta|^2\right] + \left[-\frac{3}{2}(p - q) + 2(p + q)\left(-\frac{1}{4}\delta^2 + \Lambda^2 + |\Delta|^2\right)\right]\Lambda, \qquad (89)$$

$$c = t(1 + \delta) + 2(p + q)\delta\Lambda, \qquad (90)$$

$$d = t2\Lambda - \frac{3}{2}(p - q) + 2(p + q)\left[\frac{1}{4}\delta^2 + \Lambda^2 + |\Delta|^2\right]. \qquad (91)$$

It is easy to check that the dimensionless coefficients $\bar{a} = a/t$, $\bar{b} = b/t$, $\bar{c} = c/t$, and $\bar{d} = d/t$ satisfy the following relations

$$\bar{a} = \bar{c}\Lambda + (\bar{p} + \bar{q})\delta\left[\frac{1}{4}(1 - \delta^2) - \Lambda^2 - |\Delta|^2\right], \qquad (92)$$

$$\bar{b} = \bar{d}\Lambda - \frac{1}{2}\bar{c}\delta - \frac{1}{4}(1 + \delta)^2 + \Lambda^2 + |\Delta|^2, \qquad (93)$$

where $\bar{p} = p/t$ and $\bar{q} = q/t$. Here $\Lambda$ and $\Delta$ are the values of the nonzero formal ensemble averages $\langle c_{l\sigma}^{\dagger}c_{l'\sigma} \rangle$ and $\langle c_{l\downarrow}c_{l'\uparrow} \rangle$ for the nearest-neighbor sites $l$ and $l'$, which will be calculated self-consistently [14]. We have assumed that $\Lambda$ and $\Delta$ do not depend on spin $\sigma$ explicitly in the discussion of superconductivity.

Using Eq (15), the ensemble averages $\langle c_{l'\uparrow}^{\dagger}c_{l\uparrow} \rangle$ and $\langle c_{l'\uparrow}^{\dagger}c_{l\downarrow}^{\dagger} \rangle$ can be calculated by the following equations

$$\langle c_{l'\uparrow}^{\dagger}c_{l\uparrow} \rangle = \frac{1}{2}\delta_{ll'} - \frac{1}{2N}\sum_{\mathbf{k}} e^{i\mathbf{k}\cdot(1-\mathbf{r})} \frac{-\mu + za + b'\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh\left(\frac{1}{2}\beta E_{\mathbf{k}}\right), \qquad (94)$$

$$\langle c_{l'\uparrow}^{\dagger}c_{l\downarrow}^{\dagger} \rangle = -\frac{\Delta^*}{2N} \sum_{\mathbf{k}} e^{i\mathbf{k}\cdot(1-\mathbf{r})} \frac{zc + d'\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh\left(\frac{1}{2}\beta E_{\mathbf{k}}\right). \qquad (95)$$

These are the general formulas of the ensemble averages for the study of superconductivity. Eqs (94) and (95) confirm that $\langle c_{l'\uparrow}^{\dagger}c_{l\uparrow} \rangle$, and $\langle c_{l'\uparrow}^{\dagger}c_{l\downarrow}^{\dagger} \rangle$ are functions of $R_l - R_{l'}$, which is a result of the translational invariance of the system.

Three quantities $\mu$, $\Lambda$, and $\Delta$ need to be determined from this formalism. The self-consistent equations for $\mu$, $\Lambda$, and $\Delta$ can be obtained from Eqs (94) and
(95) with the following form

$$ \delta = \frac{1}{N} \sum_{\mathbf{k}} \frac{-\mu + za + b'\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh \left( \frac{1}{2}\beta E_{\mathbf{k}} \right). \qquad (96) $$

$$ 2z_1\Lambda = -\frac{1}{N} \sum_{\mathbf{k}} \frac{(-\mu + za + b'\gamma_{\mathbf{k}})\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh \left(\frac{1}{2}\beta E_{\mathbf{k}}\right), \qquad (97) $$

$$ 2z_2 = -\frac{1}{N} \sum_{\mathbf{k}} \frac{(zc + d'\gamma_{\mathbf{k}})\gamma_{\mathbf{k}}}{E_{\mathbf{k}}} \tanh \left(\frac{1}{2}\beta E_{\mathbf{k}}\right), \qquad (98) $$

where $z_1 = z(1 - \delta)/(1 + r)$ and $z_2 = z(1 - \delta)r/(1 + r)$, the approximate average numbers of the nearest neighbors of the site $l$ with the same and opposite spin, respectively.

Furthermore, introducing a new variable $x \to \gamma_{\mathbf{k}}/z$, the self-consistent equations can be written in the integration form

$$ \delta = zb' \int_{-1}^{1} dx \rho(x) \frac{x - x_0}{E(x)} \tanh \left( \frac{1}{2} z \beta E(x) \right), \qquad (99) $$

$$ 2z_1\Lambda = -z^2b' \int_{-1}^{1} dx\rho(x) \frac{(x-x_0)x}{E(x)} \tanh \left(\frac{1}{2}z\beta E(x)\right), \qquad (100) $$

$$ 2z_2 = -z^2 d' \int_{-1}^{1} dx \rho(x) \frac{(x - x'_0)x}{E(x)} \tanh \left( \frac{1}{2} z \beta E(x) \right), \qquad (101) $$

where $x_0 = (\mu/z - a)/b', x'_0 = -c/d'$, the density of states $\rho(x)$ is given in Eq (28), and the elementary excitation energy can be written as

$$ E(x) = \sqrt{b'^2(x - x_0)^2 + d'^2(x - x'_0)^2 \Delta^2}. \qquad (102) $$

In the dimensionless system where the coefficients $a, b, c$, and $d$ are rescaled with $t$, there are only two independent parameters $\bar{p}$ and $\bar{q}$. Except for the expression of density of state (28), all equations are applicable to any dimensional systems.

Eqs (99)-(101) are the basic self-consistent equations for the superconductivity state of the model. There are three unknowns, $x_0$, $\Lambda$, and $\Delta$, to be determined from the system of equations. The most interesting quantity is the energy gap parameter $\Delta$. It is the order parameter of the superconducting state. The quantity $x_0$, the effective chemical potential of the system, is relevant to the features of the Fermi surface. In the case where $d \neq 0$ and $x_0 \neq x'_0$, a nonzero $\Delta$ causes an energy gap at the Fermi surface, as shown in Eq (102). This is the most remarkable feature of the superconducting state. The system of equations (99)-(101) can be solved numerically. The phase diagram of the superconductivity for some typical values of $\bar{q}$ and $r$ is shown in Fig. 2. Other detailed results, such as the $\delta$ and $T$ dependence of the energy gap parameter $\Delta$, will be presented in the forthcoming papers.
<figure><img src="image_3.png" /><figcaption>Figure 2: Phase diagram of superconductivity state for typical values of $\bar{q} = q/t$ and $r$. The region below each phase boundary line (within each dome) is the superconductivity state.</figcaption></figure>

# 7 Pseudogap State

We believe that the pseudogap state observed in experiments can be described by the special case of the system of equations (99)-(101) where $x_0 = x'_0$. It can be shown that $b = \Lambda d$ in this case, and therefore the self-consistent equations take the following form

$$ \delta = \frac{z r_1 \Lambda}{\kappa} \int_{-1}^{1} dx \rho(x) \tanh \left( \frac{1}{2} z \beta \kappa d(x - x_0) \right), \qquad (103) $$

$$ 2 = -z \frac{1}{\kappa} \int_{-1}^{1} dx \rho(x) x \tanh \left( \frac{1}{2} z \beta \kappa d(x - x_0) \right), \qquad (104) $$

$$ \Lambda^2 + \Delta^2 = \frac{1}{4}(1 + \delta)^2 + (\bar{p} + \bar{q})\delta^2\Lambda, \qquad (105) $$

where $\kappa$ is given by

$$ \kappa = \sqrt{(r_1 \Lambda)^2 + (r_2 \Delta)^2}. \qquad (106) $$

Note that there is no singularity in Eqs (103)-(104) at any doping level as $T \to 0$. At a given temperature $T$ and doping level $\delta$, the quantities $x_0$, $\Lambda$, and $\Delta$ are determined from Eqs (103)-(105).

In the pseudogap state, a new degree of freedom, $n^*$, which is the number of electrons that cannot ‘see’ the energy gap, emerges. This is the density of the
charge carriers contributing to the electrical conductivity in the pseudogap region. Explicitly, the condition of $x_0 = x'_0$ can be expressed as

$$ \bar{c} + x_0 r_2 \bar{d} = 0, \qquad (107) $$

where

$$ \bar{c} = 2 - n + 2(\bar{p} + \bar{q})(1 - n)\Lambda, \qquad (108) $$

$$ \bar{d} = 2\Lambda - \bar{p} + 2\bar{q} + 2(\bar{p} + \bar{q}) \left[ -\frac{1}{4}n(2-n) + \Lambda^2 + \Delta^2 \right]. \qquad (109) $$

When $n$ is the total number of the electrons, $n = 1 - \delta$, Eq (107) describes the line of $T^*(\delta)$ [15]. If we require the condition (107) to be satisfied for $T < T^*$, then the only possibility is that $n < 1 - \delta$. We assume that the new degree of freedom $n^*$ is determined by Eq (107) when $T < T^*$. That is, the condition $x_0 = x'_0$ determines the number of charge carriers $n^*$ in the pseudogap state.

<figure><img src="image_4.png" /><figcaption>Figure 3: Temperature dependence of the charge carriers density $n^*$ in pseudogap state for typical doping levels ($\bar{q} = 0.63$, $r = 2.5$).</figcaption></figure>

Replacing $n$ by $n^*$ in Eq (107) leads to a quadratic equation of $n^*$, which can
be solved with the result

$$
\begin{align}
n^* &= 1 - \delta^*, \nonumber \\
\delta^* &= \frac{-1}{(\bar{p} + \bar{q})r_2 x_0} \left[ 1 + 2(\bar{p} + \bar{q})\Lambda + \sqrt{D} \right], \nonumber \\
D &= \left[ 1 + 2(\bar{p} + \bar{q})\Lambda \right]^2 - 2(\bar{p} + \bar{q})r_2 x_0 \left[ 1 + 2 \left( -\frac{3}{4}(\bar{p} - \bar{q}) + \Lambda + (\bar{p} + \bar{q})(\Lambda^2 + \Delta^2) \right) r_2 x_0 \right]. \tag{110}
\end{align}
$$

The symbol $\delta^*$ should not be confused with the doping parameter. As Fig. 3 shows, at a doping level $\delta$, $n^*$ increases linearly with temperature $T$ for the range of higher $T$. This remarkable result indicates that the Hall coefficient $R_H = 1/en^*$ has a temperature dependence of $1/T$.

As expected, the condition of $n^* = 0$ should give the metal-insulator (MI) transition line. A semi-phenomenological analysis can show that the resistivity in the pseudogap region $\rho^* \sim -1/x_0$ at the Fermi surface. Detailed discussion of the resistivity and the phase diagram of the pseudogap state (the MI transition line and $T^*$ line) will be presented in forthcoming papers.

## 8 Conclusion

We have established a novel model Hamiltonian for the strongly correlated electron systems, which contains the Coulomb potential energy terms up to the order of two-site interactions. The magnetic properties (ferromagnetism and antiferromagnetism) of this model have been studied. Our results indicate that the system is not necessarily in the ferromagnetic state as $T \to 0$ at any doping level. Nevertheless, the system must be in an antiferromagnetic state with a nonzero magnetization at the state point of $T \to 0$ and $\delta = 0$. Moreover, the system exhibits a superconducting state in a doped region at sufficiently low temperatures. The phase diagrams for antiferromagnetism and superconductivity of the model are presented. This model also predicts the temperature dependence of $1/T$ of the Hall coefficient in the pseudogap state.

## References

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

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

[3] M. C. Gutzwiller, Phys. Rev. Lett. **10**, 159 (1963).
[4] https://arxiv.org/pdf/0706.4236v1.pdf J. Spalek, *t-J model then and now: A personal perspective from the pioneering times*.

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

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

[7] Gerald D. Mahan, *Many-Particle Physics* (Plenum Press, New York, 1990).

[8] 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.

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

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

[11] 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.

[12] D. N. Zubarev, Sov. Phys. Uspekhi **3**, 320 (1960).

[13] A. L. Fetter and J. D. Walecka, *Quantum Theory of Many-Particle Systems* (McGraw-Hill, San Francisco, 1971).

[14] The ensemble average can only be performed *formally* at this stage, because we do not know the quantity $\Lambda_{\bar{\sigma}}$, yet. In fact, $\Lambda_{\bar{\sigma}}$ must be calculated *self-consistently* from the formalism when the Green function is established. That is, $\Lambda_{\bar{\sigma}}$ both determines and is determined by the Green functions. See Chap. 13 of Ref. [13] for details. This understanding is applicable to the quantity $\Delta$, which appears in the discussion of superconductivity.

[15] $T^*$ line should be determined by either 1) all of the electrons in the system cannot ‘see’ the energy gap ($n^* = n = 1 - \delta$) or 2) the energy gap disappears ($\Delta = 0$, and therefore no electron can ‘see’ the energy gap), whichever occurs first. Our study shows that the condition 2), i.e., $\Delta = 0$, occurs first.