format: text

text:
# Effect of Exclusion of Double Occupancies in 𝑡-𝐽 Model: Extension of Gutzwiller Approximation

## I Introduction

For understanding the basic physics of high-$`T_\mathrm{c}`$ superconductivity, the two-dimensional $`t`$-$`J`$ model has been extensively studied as one of the most promising and simplified models which describe charge and spin dynamics in the CuO2 plane. Mean-field theories and numerical calculations have shown that the d$`_{x^2y^2}`$ wave superconductivity (SC) takes place in a reasonable parameter region in the phase diagram ($`J/t=0.3`$ and the doping rate $`\delta <0.3`$). Therefore, as far as the d-wave SC is concerned, the experimentally observed phase diagram of high-$`T_\mathrm{c}`$ superconductors and the results in the $`t`$-$`J`$ model agree quite well.

Recently, however, there are some experiments which indicate the importance of the interplay between the antiferromagnetism and d-wave SC. The most interesting phenomenon is a realization of stripe state in some materials, which has charge ordering as well as incommensurate antiferromagnetic (AF) spin ordering. Other potential problems are the AF state induced around nonmagnetic impurities and a possibility of AF vortex cores.

In order to study these problems, it is necessary to develop a theory in which the d-wave SC and the AF correlation are treated in a reliable way. Usual mean-field theories generally overestimate the AF long-range order, so that they give unphysical results even in the uniform case. In this paper we develop a new type of analytic theory which gives a reliable estimate of the AF correlation as well as the d-wave SC. Our scheme is an extension of the Gutzwiller approximation.

In the $`t`$-$`J`$ model the double occupancy of up- and down-spin electrons at any site is prohibited. To study its ground state, it is natural to use a projected mean-field wave function in which the double occupancies are excluded. As the mean-field wave functions, a SDW mean-field solution at half filling, BCS wave functions, and a coexistent state of the AF and d-wave SC have been used. Variational energies are calculated in the variational Monte Carlo (VMC) simulations, which treat exactly the constraints of no double occupancies. Among the above trial states, the coexistent state has the lowerest variational energy in the doping region $`\delta <0.1`$.

On the other hand, the analytic theories in which the constraint is treated approximately give large discrepancies in the estimation of the AF correlation. For example, in the slave-boson mean-field theory which takes account of the constraint as an average, the AF order is too overestimated and it extends up to unphysical doping rates ($`\delta <0.2`$). Thus in the slave-boson theory, we are unable to discuss the stripe state, for example, because it is stabilized near $`\delta =0.125`$ doping where the slave-boson theory gives the AF state even for the bulk.

The Gutzwiller approximation is more advanced approximation than the simple mean-field theories in treating the constraint. Renormalizations of expectation values are introduced by comparing the statistical weighting factors in the wave functions with and without projection. As a result, the parameters $`t`$ and $`J`$ are renormalized to $`g_tt`$ and $`g_sJ`$. It has been shown that the Gutzwiller approximation gives a fairly reliable estimation of the variational energy for the pure d-wave SC state when it is compared with the VMC results. However it was shown that there is no region in the phase diagram where the AF state is stabilized. This contradicts with the VMC results.

In order to clarify the origin of these discrepancies, we investigated the VMC data and found that the Gutzwiller approximation has to be modified when the AF correlations are present. Based on these observations, we extend the Gutzwiller approximation in this paper and derive an analytic formalism for the renormalization factors which reproduces the VMC results. We show that it is important to take account of the longer-range correlations for the weighting factors, in contrast to the previous approximation where only the site-diagonal expectation values such as $`n_{i\sigma }`$ are considered. As a result, the renormalization factors for $`t`$ and $`J`$ have nonlinear dependences on the expectation values of Cooper pairs and AF moment. We think that this is the essence of the strong correlation in the sense that the renormalization appears solely from the projection operators.

The merit of the present scheme is that it can be easily applied to the inhomogeneous systems, where the VMC simulations are not so feasible. Thus a reliable analytic formulation can be given to the problems, such as the stripe state, vortex cores and states around impurities, where the interplay between the AF and d-wave SC plays an important role. Preliminary results of these problems have been described elsewhere.

The outline of this paper is as follows: In Sec. II we briefly review the formulation of the original Gutzwiller approximation in order to make the present paper self-contained. We also show our final results for the extended Gutzwiller approximation before going into the details of the derivation. In Sec. III, the original Gutzwiller approximation is extended in order to include longer-range correlations. We divide the whole system into cells and introduce the configurations in each cell for evaluating the weighting factors. Using the general formulation developed in Sec. III, we obtain the renormalization factors for the half-filled case in Sec. IV. We approximate the weighting factors in a perturbative way with respect to the nearest-neighbor correlations. In IV.D, the physical origin of the enhancement of the renormalization factor for the $`z`$-component of exchange interaction is discussed in a viewpoint of statistical weights of real-space spin configurations. In Sec. V, the results at half filling are generalized to the less-than-half-filled case. The present method is applied to the projected variational states in Sec. VI. The self-consistency equations are derived and solved numerically to show that they give reasonable estimate of AF long-range order near half filling. Finally section VII is devoted to a summary and discussions on related problems.

## II $`t`$-$`J`$ model and the original Gutzwiller approximation

We consider the Gutzwiller approximation for the two-dimensional $`t`$-$`J`$ model on a square lattice:

$$\widehat{}=t\underset{ij\sigma }{}P_G(c_{i\sigma }^{}c_{j\sigma }+h.c.)P_G+J\underset{ij}{}𝑺_i𝑺_j,$$
(1)
where $`ij`$ represents the sum over the nearest-neighbor sites and $`𝑺_i=c_{i\alpha }^{}(\frac{1}{2}𝝈)_{\alpha \beta }c_{i\beta }`$. The Gutzwiller’s projection operator $`P_G`$ is defined as $`P_G=_j(1\widehat{n}_j\widehat{n}_j)`$.

For this Hamiltonian, we assume the projected BCS-SDW mean-field wave function

$$|\psi =P_G|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V}),$$
(2)
where $`\mathrm{\Delta }_d^\mathrm{V}`$, $`\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V}`$ and $`\mu ^\mathrm{V}`$ are the variational parameters relating to d-wave SC, AF and chemical potential, respectively, and $`|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V})`$ is a Hartree-Fock type wave function with d-wave SC and AF orders. The wave function $`|\psi `$ is a natural generalization of the RVB state proposed by Anderson. It has been shown that the coexistent state of AF and d-wave SC has the best variational energy in the VMC simulations for the doping rate $`\delta <0.1`$.

In evaluating the variational energy, the projection operator makes difficulties for an analytic approach. The Gutzwiller approximation was developed for this purpose. In this method, the effect of the projection is taken into account by renormalizations of expectation values as follows:

$$c_{i\sigma }^{}c_{j\sigma }=g_tc_{i\sigma }^{}c_{j\sigma }_0,𝑺_i𝑺_j=g_s𝑺_i𝑺_j_0,$$
(3)
where $`\mathrm{}_0`$ represents the expectation value in terms of $`|\psi _0=|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V})`$, and $`\mathrm{}`$ represents the normalized expectation value in $`|\psi =P_G|\psi _0`$;

$`\widehat{𝒪}{\displaystyle \frac{\psi |\widehat{𝒪}|\psi }{\psi |\psi }}={\displaystyle \frac{\psi _0|P_G\widehat{𝒪}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}.`$ (4)
The coefficients $`g_t`$ and $`g_s`$ are the renormalization factors of the expectation values, which we call as Gutzwiller factors in the following. Using these notations, the variational energy $`E_{\mathrm{var}}=\widehat{}`$ is rewritten as

$$E_{\mathrm{var}}=\widehat{}_{\mathrm{eff}}_0,$$
(5)
where the parameters $`t`$ and $`J`$ in $`\widehat{}`$ are replaced with

$$t_{\mathrm{eff}}=g_tt,J_{\mathrm{eff}}=g_sJ.$$
(6)
A systematic estimation of $`g_t`$ and $`g_s`$ was developed by Ogawa et al. whose method is briefly reviewed in II.A. For the case with AF order parameter, it can be shown that

$$g_t=\frac{2\delta (1\delta )}{1\delta ^2+4m^2},g_s=\frac{4(1\delta )^2}{(1\delta ^2+4m^2)^2},$$
(7)
where $`\delta `$ is the hole concentration, $`\delta =1n`$, and $`m`$ is the expectation value of AF order parameter in terms of $`|\psi _0`$ defined as

$$m=\frac{(1)^j}{2}\left(\widehat{n}_j_0\widehat{n}_j_0\right).$$
(8)
These Gutzwiller factors, however, do not reproduce the VMC results, as mentioned in I. We find that this discrepancy is due to the fact that only the site-diagonal expectation values are taken into account in obtaining Eq. (7). In this paper we extend the method by Ogawa et al. systematically to include the longer range correlations. We show that it is important to include the effects of the nearest-neighbor expectation values, such as

$`\chi `$ $`=c_{i\sigma }^{}c_{j\sigma }_0,`$ (9)
$`\mathrm{\Delta }`$ $`=c_i^{}c_j^{}_0,`$ (10)
in the estimation of $`g_t`$ and $`g_s`$. This corresponds to taking account of the effect of surroundings of the corresponding bond.

Furthermore we find that the Gutzwiller factors $`g_s`$ for the exchange interaction have different values for the $`z`$ component (denoted as $`g_s^Z`$) and $`xy`$ component ($`g_x^{XY}`$); i.e.,

$`S_i^zS_j^z`$ $`=g_s^ZS_i^zS_j^z_0,`$ (11)
$`S_i^+S_j^{}`$ $`=g_s^{XY}S_i^+S_j^{}_0,`$ (12)
instead of (3) in the presence of $`m`$. Actually from the VMC calculations, we have estimated the behaviors of $`g_s^Z`$ and $`g_s^{XY}`$ as a function of $`m`$, using the relations

$$g_s^Z=\frac{S_i^zS_j^z}{S_i^zS_j^z_0},g_s^{XY}=\frac{S_i^+S_j^{}}{S_i^+S_j^{}_0},$$
(13)
where the numerators $`S_i^zS_j^z`$ and $`S_i^+S_j^{}`$ are evaluated numerically in VMC simulations. It was found that $`g_s^Z`$ is enhanced compared with $`g_s^{XY}`$ and that this enhancement is essential for the stabilization of the AF order. In the present extended Gutzwiller approximation, this result is reproduced. Physically the weighting factor of $`g_s^Z`$ for a specific configuration of a bond is enhanced due to the effect of the surrounding AF correlations.

Before going into details of the derivation, we show our final results. They are summarized as follows:

$`g_s^{XY}=\left({\displaystyle \frac{2(1\delta )}{1\delta ^2+4m^2}}\right)^2a^7,`$ (14)
and

$$g_s^Z=g_s^{XY}\frac{1}{4m^2+X_2}\left[X_2+4m^2\left\{1+\frac{6X_2(1\delta )^2}{1\delta ^2+4m^2}a^3\right\}^2\right],$$
(15)
where

$`a`$ $`=1+{\displaystyle \frac{4X}{(1\delta ^2+4m^2)^2}},`$ (16)
$`X`$ $`=2\delta ^2(\mathrm{\Delta }^2\chi ^2)+8m^2(\chi ^2+\mathrm{\Delta }^2)+4(\chi ^2+\mathrm{\Delta }^2)^2,`$ (17)
$`X_2`$ $`=2(\chi ^2+\mathrm{\Delta }^2).`$ (18)
The above expressions sufficiently reproduce the results obtained in the VMC simulations. However we find that the slight difference can be improved by using

$$X=2\delta ^2(\mathrm{\Delta }^2\chi ^2)+8m^2(\chi ^2+\mathrm{\Delta }^2)+2(\chi ^2+\mathrm{\Delta }^2)^2,$$
(19)
instead of Eq. (18). A typical $`m`$-dependence of $`g_s^Z`$ and $`g_s^{XY}`$ are shown in Fig. 1 for the half-filled case, $`\delta =0`$, and with $`\mathrm{\Delta }`$ being fixed at the two values, $`\mathrm{\Delta }=0.02`$ and $`0.18`$. The behaviors in Fig. 1 agree with those obtained in the VMC simulations using the relation in Eq. (13). The meanings of the rather complicated Gutzwiller factors become apparent in the following sections.

The Gutzwiller factor for the hopping term is given by

$`g_t={\displaystyle \frac{2\delta (1\delta )}{1\delta ^2+4m^2}}{\displaystyle \frac{(1+\delta )^24m^22X_2}{(1+\delta )^24m^2}}a.`$ (20)
Figure 2 shows the $`m`$-dependences of $`g_s^{XY},g_s^Z`$ and $`g_t`$ for $`\delta =0.12`$, which can be compared with the VMC results in Ref. . Note here that the slave-boson mean-field theory simply gives $`g_t=b_i^{}b_j=\delta `$ and $`g_s=1`$. The nonlinear dependences of $`g_s^{XY},g_s^Z`$ and $`g_t`$ on $`\mathrm{\Delta },\chi `$ and $`m`$ are beyond the slave boson theory. In the following subsection we review the original formulation of the Gutzwiller approximation which is useful for the generalization in the later sections.

### A Original formulation of Gutzwiller approximation

Let us start with the formulation of the Gutzwiller approximation developed by Ogawa et al  when it is applied to the $`t`$-$`J`$ model. In this subsection we do not consider the AF order. The ultimate aim is to evaluate the variational energy

$$E_{\mathrm{var}}=\widehat{}=\frac{\psi |\widehat{}|\psi }{\psi |\psi },$$
(21)
using $`|\psi =P_G|\psi _0`$. First, Ogawa et al rewrote the denominator as ($`P_G^2=P_G`$)

$`\psi |\psi `$ $`=\psi _0|P_GP_G|\psi _0`$ (25)
$`=\psi _0|{\displaystyle \underset{j}{}}(1\widehat{n}_j\widehat{n}_j)|\psi _0`$
$`=\psi _0|{\displaystyle \underset{j}{}}\{\widehat{n}_j(1\widehat{n}_j)+\widehat{n}_j(1\widehat{n}_j)`$
$`+(1\widehat{n}_j)(1\widehat{n}_j)\}|\psi _0.`$
The expansion of the product in (25) leads to various real-space configurations. Therefore $`\psi |\psi `$ can be rewritten as

$`{\displaystyle \underset{\mathrm{config}.}{}}\psi _0`$ $`|{\displaystyle \underset{j𝒜}{}}\widehat{n}_j(1\widehat{n}_j){\displaystyle \underset{j^{}}{}}\widehat{n}_j^{}(1\widehat{n}_j^{})`$ (27)
$`\times {\displaystyle \underset{j^{\prime \prime }}{}}(1\widehat{n}_{j^{\prime \prime }})(1\widehat{n}_{j^{\prime \prime }})|\psi _0,`$
where $`𝒜`$ ($``$) is a subset of the lattice sites which is singly occupied by an up-spin (down-spin) electron and $``$ is a subset which is vacant. The summation is over all the possible configurations, i.e., all the possible combinations of $`\{𝒜,,\}`$.

In principle, Eq. (27) can be calculated using Wick’s theorem since $`|\psi _0`$ is a mean-field wave function, although it is very complicated. The simplest Gutzwiller approximation is to estimate (27) by using only the site-diagonal expectation values, which gives

$`{\displaystyle \underset{\mathrm{config}.}{}}`$ $`{\displaystyle \underset{j𝒜}{}}\widehat{n}_j(1\widehat{n}_j)_0{\displaystyle \underset{j^{}}{}}\widehat{n}_j^{}(1\widehat{n}_j^{})_0`$ (28)
$`\times `$ $`{\displaystyle \underset{j^{\prime \prime }}{}}(1\widehat{n}_{j^{\prime \prime }})(1\widehat{n}_{j^{\prime \prime }})_0,`$ (29)
with $`\mathrm{}_0`$ meaning the expectation value in $`|\psi _0`$.

For example, the contribution from the site $`j𝒜`$ is evaluated as

$$\omega _A\widehat{n}_j(1\widehat{n}_j)_0=\frac{n}{2}\left(1\frac{n}{2}\right),$$
(30)
where $`n`$ is the average electron density, $`n=N_e/N`$, with $`N`$ ($`N_e`$) being the total number of sites (total electron number). We call $`\omega _A`$ as the weight for $`j`$ site belonging to the subset $`𝒜`$. Since the hole density is $`\delta =1n`$, we obtain

$`\omega _A`$ $`=\omega _B={\displaystyle \frac{n}{2}}\left(1{\displaystyle \frac{n}{2}}\right)={\displaystyle \frac{1\delta ^2}{4}},`$ (31)
$`\omega _E`$ $`=\left(1{\displaystyle \frac{n}{2}}\right)^2={\displaystyle \frac{(1+\delta )^2}{4}}.`$ (32)
Then the expectation value (29) is rewritten as

$$\underset{\mathrm{config}.}{}\omega _A^{N_A}\omega _B^{N_B}\omega _E^{N_E}=\frac{N!}{N_A!N_B!N_E!}\omega _A^{N_A}\omega _B^{N_B}\omega _E^{N_E},$$
(33)
where $`N_A`$ is the number of sites in the subset $`𝒜`$, and so on, ($`N_A+N_B+N_E=N`$). In this simple case, $`N_A`$ is equal to the number of up-spin electrons, so that

$$N_A=N_B=\frac{N_e}{2},N_E=NN_e.$$
(34)
A similar procedure can be carried out for the estimation of the numerator in $`E_{\mathrm{var}}`$. For the exchange term, we have

$`\psi _0`$ $`|P_GS_{\mathrm{}}^+S_m^{}P_G|\psi _0`$ (39)
$`=\psi _0|S_{\mathrm{}}^+S_m^{}{\displaystyle \underset{j\mathrm{},m}{}}\{\widehat{n}_j(1\widehat{n}_j)+\widehat{n}_j(1\widehat{n}_j)`$
$`+(1\widehat{n}_j)(1\widehat{n}_j)\}|\psi _0`$
$`={\displaystyle \underset{\mathrm{config}.}{}}\psi _0|S_{\mathrm{}}^+S_m^{}{\displaystyle \underset{j𝒜^{}}{}}\widehat{n}_j(1\widehat{n}_j){\displaystyle \underset{j^{}^{}}{}}\widehat{n}_j^{}(1\widehat{n}_j^{})`$
$`\times {\displaystyle \underset{j^{\prime \prime }^{}}{}}(1\widehat{n}_{j^{\prime \prime }})(1\widehat{n}_{j^{\prime \prime }})|\psi _0.`$
Here the summation is over all the possible configurations $`\{𝒜^{},^{},^{}\}`$ in which the two sites $`\mathrm{}`$ and $`m`$ are excluded. The evaluation of (39) with the site-diagonal expectation values for $`j\mathrm{},m`$ leads to

$$\frac{(N2)!}{N_A^{}!N_B^{}!N_E^{}!}\omega _A^{N_A^{}}\omega _B^{N_B^{}}\omega _E^{N_E^{}}S_{\mathrm{}}^+S_m^{}_0.$$
(40)
Since the sites $`\mathrm{}`$ and $`m`$ contain one up-spin electron and one down-spin electron, we have $`N_A^{}=\frac{N_e}{2}1`$, $`N_B^{}=\frac{N_e}{2}1`$, and $`N_E^{}=N_E=NN_e`$.

By combining the estimation of the denominator in Eq. (33), we obtain the original Gutzwiller approximation

$`S_{\mathrm{}}^+S_m^{}`$ $`={\displaystyle \frac{\psi _0|P_GS_{\mathrm{}}^+S_m^{}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}`$ (43)
$`={\displaystyle \frac{\frac{N_e}{2}\frac{N_e}{2}}{N(N1)}}\omega _A^1\omega _B^1S_{\mathrm{}}^+S_m^{}_0`$
$`={\displaystyle \frac{4}{(1+\delta )^2}}S_{\mathrm{}}^+S_m^{}_0.`$
This gives $`g_s^{XY}=4/(1+\delta )^2`$.

### B Antiferromagnetic case

Ogawa et al also considered the above formalism in the case with AF long-range order for the Hubbard model. The application to the $`t`$-$`J`$ model was carried out by Zhang et al. In this subsection we follow their methods to obtain the Gutzwiller approximation in the $`t`$-$`J`$ model.

In the AF case, the sublattices 1 and 2 are distinguished and their magnetizations are defined as

$$\frac{1}{2}\left(\widehat{n}_j_0\widehat{n}_j_0\right)=m,$$
(44)
for the sublattice 1 and $`m`$ for the sublattice 2. Thus we denote

$`\widehat{n}_j_0`$ $`={\displaystyle \frac{n}{2}}+mr,`$ (45)
$`\widehat{n}_j_0`$ $`={\displaystyle \frac{n}{2}}mw,`$ (46)
for the sublattice 1, where $`r`$ ($`w`$) means the average electron density with the right (wrong) spin direction on the sublattice. $`r`$ and $`w`$ are exchanged for the sublattice 2.

Using these notations, the denominator of $`E_{\mathrm{var}}`$ is evaluated as

$`\psi _0|P_GP_G|\psi _0=`$ $`{\displaystyle \underset{\mathrm{config}.}{}}\omega _{A1}^{N_{A1}}\omega _{B1}^{N_{B1}}\omega _{E1}^{N_{E1}}\omega _{A2}^{N_{A2}}\omega _{B2}^{N_{B2}}\omega _{E2}^{N_{E2}}`$ (47)
$`=`$ $`{\displaystyle \underset{N_{A1},N_{B1}}{}}{\displaystyle \frac{\left(\frac{N}{2}\right)!}{N_{A1}!N_{B1}!N_{E1}!}}{\displaystyle \frac{\left(\frac{N}{2}\right)!}{N_{A2}!N_{B2}!N_{E2}!}}`$ (49)
$`\times \omega _{A1}^{N_{A1}}\omega _{B1}^{N_{B1}}\omega _{E1}^{N_{E1}}\omega _{A2}^{N_{A2}}\omega _{B2}^{N_{B2}}\omega _{E2}^{N_{E2}},`$
where the configurations in the sublattice 1 are specified by $`(N_{A1},N_{B1},N_{E1})`$, and so on. Contrary to the previous subsection, Eq. (49) has a summation over possible values of $`N_{A1}`$ and $`N_{B1}`$ because there are only four constraints between six numbers, $`(N_{A1},N_{B1},N_{E1},N_{A2},N_{B2},N_{E2})`$:

$`N_{A1}+N_{B1}+N_{E1}`$ $`={\displaystyle \frac{N}{2}},N_{A2}+N_{B2}+N_{E2}={\displaystyle \frac{N}{2}},`$ (50)
$`N_{A1}+N_{A2}`$ $`={\displaystyle \frac{N_e}{2}},N_{B1}+N_{B2}={\displaystyle \frac{N_e}{2}}.`$ (51)
The weights in the site-diagonal expectation values are

$`\omega _{A1}`$ $`=\widehat{n}_j_01\widehat{n}_j_0=r(1w),`$ (53)
$`\omega _{B1}`$ $`=\widehat{n}_j_01\widehat{n}_j_0=w(1r),`$ (54)
$`\omega _{E1}`$ $`=1\widehat{n}_j_01\widehat{n}_j_0=(1r)(1w),`$ (55)
for the sublattice 1, and

$`\omega _{A2}`$ $`=\omega _{B1},`$ (56)
$`\omega _{B2}`$ $`=\omega _{A1},`$ (57)
$`\omega _{E2}`$ $`=\omega _{E1},`$ (58)
for the sublattice 2.

Since $`N_{A1},N_{B1}`$ and $`N`$ are large numbers, Eq. (49) is further approximated by the largest term in the summation. By taking the partial derivative of the logarithm of each term in (49) with respect to $`N_{A1}`$, we obtain

$$\mathrm{ln}\overline{N_{A1}}+\mathrm{ln}(\frac{N_e}{2}\overline{N_{A1}})+\mathrm{ln}\omega _{A1}\mathrm{ln}\omega _{A2}=0,$$
(59)
with $`\overline{N_{A1}}`$ being the value of $`N_{A1}`$ which gives the largest term in the summation. By solving (59), $`\overline{N_{A1}}`$ is given by

$$\overline{N_{A1}}=\frac{r(1w)}{r(1w)+w(1r)}\frac{N_e}{2}=\frac{r(1w)}{n2rw}\frac{N_e}{2}.$$
(60)
Similarly we obtain

$`\overline{N_{B1}}={\displaystyle \frac{w(1r)}{n2rw}}{\displaystyle \frac{N_e}{2}}.`$ (61)
The summation in (49) is approximated by the term with $`\overline{N_{A1}}`$ and $`\overline{N_{B1}}`$.

Again the numerator such as $`\psi _0|P_GS_{\mathrm{}}^+S_m^{}P_G|\psi _0`$ is approximated in the similar way. As a result, the Gutzwiller approximation in the presence of AF order, $`m`$, becomes

$`S_{\mathrm{}}^+S_m^{}`$ $`={\displaystyle \frac{\psi _0|P_GS_{\mathrm{}}^+S_m^{}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}`$ (65)
$`={\displaystyle \frac{\overline{N_{A1}}\overline{N_{B1}}}{\frac{N}{2}\frac{N}{2}}}\omega _{A1}^1\omega _{B1}^1S_{\mathrm{}}^+S_m^{}_0`$
$`={\displaystyle \frac{n^2}{(n2rw)^2}}S_{\mathrm{}}^+S_m^{}_0`$
$`={\displaystyle \frac{4(1\delta )^2}{(1\delta ^2+4m^2)^2}}S_{\mathrm{}}^+S_m^{}_0,`$
where a useful relation

$$\frac{n}{n2rw}=\frac{2(1\delta )}{1\delta ^2+4m^2}$$
(66)
has been used. Equation (65) gives

$$g_s^{XY}=\frac{4(1\delta )^2}{(1\delta ^2+4m^2)^2},$$
(67)
which is simply Eq. (7).

It is straightforward to see

$$g_s^Z=g_s^{XY},$$
(68)
and

$$g_t=\frac{2\delta (1\delta )}{1\delta ^2+4m^2}.$$
(69)
## III Extension of the Gutzwiller approximation: Formulation

The approximation in the preceding section is restricted to the site-diagonal (on-site) expectation values. It was shown that this approximation is not good enough to reproduce the AF order obtained in the VMC simulation. A straightforward extension of the previous Gutzwiller approximation is to consider longer-range expectation values. The nearest-neighbor expectation was considered before. In this section, we develop a general formalism to take account of the longer-range effects systematically.

In evaluating the denominator of $`E_{\mathrm{var}}`$, i.e., $`\psi _0|P_GP_G|\psi _0`$, various configurations have been classified depending on the state ($`,`$ or a hole) on each site as $`\{𝒜,,\}`$ in the previous section. Instead of this we divide the whole system into cells consisting of $`N_c`$ sites. (Later we take $`N_c`$ as a large enough number.) For each cell, there are $`3^{N_c}`$ configurations because each site has $`,`$ or a hole. We denote these configurations as states of a cell. Using these states, the denominator of $`E_{\mathrm{var}}`$ can be approximated as

$`\psi _0|P_GP_G|\psi _0`$ $`={\displaystyle \underset{\mathrm{all}\mathrm{the}\mathrm{possible}\mathrm{states}}{}}{\displaystyle \underset{i=1}{\overset{K}{}}}\omega _i^{N_i}`$ (71)
$`={\displaystyle \underset{\{N_i\}}{}}{\displaystyle \frac{\left(\frac{N}{N_c}\right)!}{_{i=1}^K(N_i)!}}{\displaystyle \underset{i=1}{\overset{K}{}}}\omega _i^{N_i},`$
with $`K=3^{N_c}`$. Here $`N_i`$ is the number of cells in the $`i`$-th state, and $`\omega _i`$ is the weight of the $`i`$-th state cell in the analogy with Eq. (49). Explicitly the weight $`\omega _i`$ is defined as

$`\omega _i=`$ $`\psi _0|{\displaystyle \underset{j}{}}\widehat{n}_j(1\widehat{n}_j){\displaystyle \underset{j^{}}{}}\widehat{n}_j^{}(1\widehat{n}_j^{})`$ (73)
$`\times {\displaystyle \underset{j^{\prime \prime }}{}}(1\widehat{n}_{j^{\prime \prime }})(1\widehat{n}_{j^{\prime \prime }})|\psi _0,`$
with $`j,j^{}`$ and $`j^{\prime \prime }`$ being the sites in the cell. To calculate $`\omega _i`$ is the most important part of the present theory which we carry out later.

The last summation in (71) is over all the possible values of $`N_i`$ under the following constraints:

$`{\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i={\displaystyle \frac{N_e}{2}},{\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i={\displaystyle \frac{N_e}{2}},`$ (74)
$`{\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}N_i=NN_e,`$ (75)
where $`n_i`$ ($`n_i`$) represents the number of up-spins (down-spins) in the $`i`$-th state cell, and $`n_{\mathrm{h}i}`$ represents the number of holes, respectively. Since a cell has $`N_c`$ sites, a relation

$`n_i+n_i+n_{\mathrm{h}i}=N_c,`$ (76)
holds.

Under the above constraints, we look for the largest term in Eq. (71) as was done in the original Gutzwiller approximation (see Eq. (59)). Since the constraints (75) are rather complicated, it is useful to introduce the Lagrange multipliers ($`\mu _{},\mu _{},\lambda `$). Taking the partial derivative of

$`\mathrm{ln}\left[{\displaystyle \frac{\left(\frac{N}{N_c}\right)!}{_{i=1}^K(N_i)!}}{\displaystyle \underset{i=1}{\overset{K}{}}}\omega _i^{N_i}\right]\mu _{}\left({\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i{\displaystyle \frac{N_e}{2}}\right)`$ (77)
$`\mu _{}\left({\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i{\displaystyle \frac{N_e}{2}}\right)\lambda \left({\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}N_iN+N_e\right),`$ (78)
we obtain

$`\overline{N_i}=\omega _i\mathrm{exp}(\mu _{}n_i\mu _{}n_i\lambda n_{\mathrm{h}i}).`$ (79)
From symmetry we find $`\mu _{}=\mu _{}=\mu `$, and using the relation (76) we get

$`\overline{N_i}=\omega _i\mathrm{e}^{\mu N_c}\mathrm{e}^{(\mu \lambda )n_{\mathrm{h}i}}.`$ (80)
Furthermore we introduce new variables $`W`$ and $`p`$ as

$`\mathrm{e}^{\mu N_c}`$ $`{\displaystyle \frac{N}{N_c}}{\displaystyle \frac{1}{W}},`$ (81)
$`\mathrm{e}^{\mu \lambda }`$ $`p.`$ (82)
Using these variables, $`\overline{N_i}`$ is rewritten as

$`\overline{N_i}={\displaystyle \frac{N}{N_c}}{\displaystyle \frac{\omega _i}{W}}p^{n_{\mathrm{h}i}}.`$ (83)
The variables $`W`$ and $`p`$ are to be determined from the constraints (75) which are equivalent to

$`{\displaystyle \underset{i=1}{\overset{K}{}}}\overline{N_i}={\displaystyle \frac{N}{N_c}},{\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}\overline{N_i}=NN_e.`$ (84)
These conditions become

$`{\displaystyle \underset{i=1}{\overset{K}{}}}{\displaystyle \frac{\omega _i}{W}}p^{n_{\mathrm{h}i}}`$ $`=1,`$ (85)
$`{\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}{\displaystyle \frac{\omega _i}{W}}p^{n_{\mathrm{h}i}}`$ $`=\delta N_c.`$ (86)
It is convenient to classify all the possible states of a cell into subgroups which contain $`j`$ holes ($`j=0,1,2,\mathrm{},N_c`$). (We call these subgroups as $`j`$-hole sectors.) Then the quantity

$`W_j={\displaystyle \underset{i\mathrm{with}j\mathrm{holes}}{}}\omega _i,`$ (87)
represents the total weight of the states in this subgroup. Using $`W_j`$, we can rewrite the constraints as

$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}{\displaystyle \frac{W_j}{W}}p^j`$ $`=1,`$ (88)
$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}j{\displaystyle \frac{W_j}{W}}p^j`$ $`=\delta N_c.`$ (89)
The first equality is rewritten as $`W=_jW_jp^j`$ so that $`W`$ represents the total weight. In the following sections, we evaluate $`W_j`$ and then determine $`W`$ and $`p`$ from Eq. (89).

The numerator in $`E_{\mathrm{var}}`$ is evaluated in a similar way. Let us consider the expectation value of an operator $`\widehat{𝒪}`$ (such as $`S_{\mathrm{}}^+S_m^{}`$) which operates on a certain part inside a cell of the $`i_0`$-th state. (We call this cell as the central cell.) Then all the configurations in the other cells are classified by $`\{N_i^{}\}`$ which is the number of cells of the $`i`$-th state. Then the expectation value becomes

$$\psi _0|P_G\widehat{𝒪}P_G|\psi _0=\underset{i_0}{}\underset{\{N_i^{}\}}{}\frac{\left(\frac{N}{N_c}1\right)!}{_{i=1}^K(N_i^{})!}\underset{i=1}{\overset{K}{}}\omega _i^{N_i^{}}\widehat{𝒪}_{i_0},$$
(90)
in the analogy with Eq. (40). Here $`\widehat{𝒪}_{i_0}`$ indicates the expectation value of $`\widehat{𝒪}`$ together with terms such as $`_i\widehat{n}_i(1\widehat{n}_i)`$ inside the central cell of the $`i_0`$-th state.

Determining the largest terms on the right-hand side and taking the ratio to the denominator, we finally obtain

$`\widehat{𝒪}{\displaystyle \frac{\psi _0|P_G\widehat{𝒪}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}`$ $`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{1}{\left(\frac{N}{N_c}\right)}}N_{i_0}\omega _{i_0}^1\widehat{𝒪}_{i_0}`$ (92)
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{p^{n_{\mathrm{h}i_0}}}{W}}\widehat{𝒪}_{i_0}.`$
The evaluation of the largest term and the derivation of the above result is shown in Appendix A.

If we consider only the site-diagonal expectation values to estimate $`\omega _i`$, Eq. (92) reproduces the original Gutzwiller approximation shown in the previous section. This is summarized in Appendix B.

## IV Half-Filled Case

To obtain the Gutzwiller factors formulated in the last section, it is necessary to evaluate the weights $`\omega _i`$ for each state. Although these weights are generally complicated, we estimate them assuming that the corrections to the original Gutzwiller approximation (reproduced in Appendix B) are small. To be more specific, we evaluate $`W_j`$ in the lowest orders with respect to the intersite correlations $`\chi `$ and $`\mathrm{\Delta }`$ defined in Eq. (10). Typical values for $`\chi `$ and $`\mathrm{\Delta }`$ are less than $`0.2`$, so that the perturbation with respect to them will be justified.

In this section we calculate the Gutzwiller approximation for the half filling, since it is less complicated than the doped case. All the cells do not contain holes so that only the weight, $`W_0`$, is to be calculated. We call this subset of states as 0-hole sector. After the half-filled case, it is rather straightforward to extend the results to the doped case.

### A Evaluation of $`W_0`$

First we calculate the weight $`\omega _i`$ of the $`i`$-th state which has $`n_i`$ up-spin electrons, $`n_i`$ down-spin electrons:

$$\omega _i=\underset{j𝒜}{}\widehat{n}_j(1\widehat{n}_j)\underset{j^{}}{}\widehat{n}_j^{}(1\widehat{n}_j^{})_0.$$
(93)
In the lowest order of $`\chi `$ and $`\mathrm{\Delta }`$ or in the zeroth order, we have to take only the site-diagonal expectation values. This gives

$`\omega _i^{(0)}=[r(1w)]^{n_{\mathrm{right}}}[w(1r)]^{n_{\mathrm{wrong}}},`$ (94)
where $`n_{\mathrm{right}}`$ ($`n_{\mathrm{wrong}}`$) means the number of sites on which the right (wrong) spices of spin direction is located depending on sublattice 1 and 2. The summation over all the possible configurations gives

$`W_0^{(0)}=[r(1w)+w(1r)]^{N_c}=(n2rw)^{N_c}.`$ (95)
In the next order with respect to $`\chi `$ and $`\mathrm{\Delta }`$, we have to consider contributions from expectation values of bonds in the cell as shown in Fig. 3. For example, consider a bond connecting the sites $`i`$ and $`j`$ (see Fig. 3(a)) and a state in which two up-spin electrons occupy both $`i`$ and $`j`$ sites. The expectation value in this bond gives

$`P_{}`$ $`=(1\widehat{n}_i)\widehat{n}_i\widehat{n}_j(1\widehat{n}_j)_0`$ (98)
$`=rw(1r)(1w)rw\chi ^2(1r)(1w)\chi ^2`$
$`r(1r)\mathrm{\Delta }^2w(1w)\mathrm{\Delta }^2+(\chi ^2+\mathrm{\Delta }^2)^2,`$
where we have used the Wick’s theorem. Generally $`\chi `$ and $`\mathrm{\Delta }`$ can be complex numbers and $`\chi ^2,\mathrm{\Delta }^2`$ in the above expression mean $`|\chi |^2,|\mathrm{\Delta }|^2`$ implicitly. In the similar way we calculate $`P_{},P_{},P_{}`$. Writing that the left-spin in the subscript indicates the spin on the sublattice 1, they are calculated as

$`P_{}`$ $`=P_{},`$ (99)
$`P_{}`$ $`=r^2(1w)^2+2r(1w)\chi ^2`$ (101)
$`+r^2\mathrm{\Delta }^2+(1w)^2\mathrm{\Delta }^2+(\chi ^2+\mathrm{\Delta }^2)^2,`$
$`P_{}`$ $`=w^2(1r)^2+2w(1r)\chi ^2`$ (103)
$`+w^2\mathrm{\Delta }^2+(1r)^2\mathrm{\Delta }^2+(\chi ^2+\mathrm{\Delta }^2)^2.`$
In the estimation of $`W_0^{(1)}`$ we need their summation

$`P_{}+P_{}+P_{}+P_{}=(n2rw)^2+X,`$ (104)
with $`X`$ defined in Eq. (18), where we have used the relations $`r=n/2+m`$ and $`w=n/2m`$. The first term $`(n2rw)^2`$ on the r.h.s. is the zeroth order contribution which has been included in $`W_0^{(0)}`$. In this sense, $`X`$ is a kind of connected (or cummulant) contribution coming from the real-space diagram in Fig. 3(a). The actual values of $`X`$ is roughly $`1/20`$ so that our perturbation scheme can be justified. Denoting the number of bonds in a cell as $`N_b`$, we have a contribution to $`W_0`$ as

$`W_0^{(1)}=N_bX(n2rw)^{N_c2},`$ (105)
where the factor $`(n2rw)^{N_c2}`$ comes from the contributions of the other $`N_c2`$ sites in the cell except for $`i`$ and $`j`$.

In the similar way, we evaluate the higher order contributions to $`W_0`$. By considering two bonds in the cell we approximate their contribution as

$`W_0^{(2)}=_{N_b}C_2X^2(n2rw)^{N_c4},`$ (106)
where $`{}_{N_b}{}^{}C_{2}^{}=N_b!/(N_b2)!2!`$ is the number of choices for the positions of the two bonds in a cell. The summation of these series leads to

$`W_0`$ $`={\displaystyle \underset{j=0}{\overset{N_b}{}}}{}_{N_b}{}^{}C_{j}^{}X^j(n2rw)^{N_c2j}`$ (108)
$`=(n2rw)^{N_c}\left(1+{\displaystyle \frac{X}{(n2rw)^2}}\right)^{N_b}.`$
Let us discuss the neglected contributions to $`W_0`$. Figures 3(b) and 3(c) show real-space diagrams whose connected expectation values can contribute to $`W_0`$. Although their contributions are not included, it can be shown that they are smaller than the terms in (108). For example, the contribution from Fig. 3(b) is in the order of $`N_b\delta ^2\chi ^4`$ and $`N_bm^2\chi ^4`$ etc., so that it is smaller than $`W_0^{(1)}`$. The contribution from Fig. 3(c) is in the order of $`X^2N_b`$ which is smaller than (106) by a factor $`1/N_b`$.

There is another effect neglected in $`W_0`$: The higher order terms in (108) become less and less correct because the number of the bonds giving $`X`$ is not exactly $`{}_{N_b}{}^{}C_{j}^{}`$ due to their overlapping. However the error is again in the order smaller than $`{}_{N_b}{}^{}C_{j}^{}`$. In this sense, Eq. (108) is a kind of a summation of most dominant terms in the order of $`X^jN_b^j`$.

### B Evaluation of $`g_s^{XY}`$

In order to obtain the Gutzwiller factor for $`S_{\mathrm{}}^+S_m^{}`$ at half filling we need to calculate

$`S_{\mathrm{}}^+S_m^{}={\displaystyle \underset{i_0}{}}{\displaystyle \frac{1}{W_0}}S_{\mathrm{}}^+S_m^{}_{i_0},`$ (109)
according to the general formulation (92). Note that $`W=W_0`$ and the central cell (which includes sites $`\mathrm{}`$ and $`m`$) does not contain holes, i.e., $`n_{\mathrm{h}i_0}=0`$. The average in a central cell of the $`i_0`$-th state, $`S_{\mathrm{}}^+S_m^{}_{i_0}`$, is evaluated in a similar perturbation scheme as in $`W_0`$. In the lowest order of $`\chi `$ and $`\mathrm{\Delta }`$, we have (Fig. 4(a))

$`S_{\mathrm{}}^+S_m^{}_0(n2rw)^{N_c2},`$ (110)
where the factor $`(n2rw)^{N_c2}`$ comes from the contributions from the sites in the central cell other than the sites $`\mathrm{}`$ and $`m`$.

In the next order of $`\chi `$ and $`\mathrm{\Delta }`$, there are contributions of $`X`$ from the bonds in the cell. This leads to

$`\stackrel{~}{N}_bXS_{\mathrm{}}^+S_m^{}_0(n2rw)^{N_c4},`$ (111)
where $`\stackrel{~}{N}_b`$ is the number of bonds which are not connected to the sites $`\mathrm{}`$ and $`m`$ directly (Fig. 4(b)). In general $`\stackrel{~}{N}_b`$ depends on the shape of the cell. However, if we choose the large enough cell, we have

$`\stackrel{~}{N}_b=N_b7,`$ (112)
as is evident from Fig. 4(b).

Although the dotted bonds in Fig. 4(b) do not give $`X`$, they may still give different contributions to $`S_{\mathrm{}}^+S_m^{}_{i_0}`$, such as

$`S_{\mathrm{}}^+S_m^{}\widehat{n}_m^{}(1\widehat{n}_m^{})_c,`$ (113)
where $`m^{}`$ is a nearest neighbor site of $`m`$ and $`\mathrm{}_c`$ means a connected expectation value excluding the disconnected terms such as $`S_{\mathrm{}}^+S_m^{}_0\widehat{n}_m^{}(1\widehat{n}_m^{})_0.`$ (The disconnected expectation values have been already taken into account in Eq. (110).) By calculating (113), however, we find that there are no connected conributions: Actually the two electron operators in $`S_{\mathrm{}}^+=c_{\mathrm{}}^{}c_{\mathrm{}}`$ have to make contractions with the electron operators on the site $`m`$ in the Wick’s expansion, which gives just the disconnected expectation value.

By summing up the series of corrections such as (111) we obtain

$`{\displaystyle \underset{i_0}{}}S_{\mathrm{}}^+S_m^{}_{i_0}`$ (114)
$`=S_{\mathrm{}}^+S_m^{}_0\{(n2rw)^{N_c2}`$ (115)
$`+\stackrel{~}{N}_bX(n2rw)^{N_c4}+\mathrm{}\}`$ (116)
$`=S_{\mathrm{}}^+S_m^{}_0(n2rw)^{N_c2}\left(1+{\displaystyle \frac{X}{(n2rw)^2}}\right)^{\stackrel{~}{N}_b}.`$ (117)
Combining with $`W_0`$, the final result for the extended Gutzwiller approximation for $`S_{\mathrm{}}^+S_m^{}`$ becomes

$`S_{\mathrm{}}^+S_m^{}={\displaystyle \frac{S_{\mathrm{}}^+S_m^{}_0}{(n2rw)^2}}\left(1+{\displaystyle \frac{X}{(n2rw)^2}}\right)^{(N_b\stackrel{~}{N}_b)}.`$ (118)
Since the Gutzwiller factor is defined as the ratio between $`S_{\mathrm{}}^+S_m^{}`$ and $`S_{\mathrm{}}^+S_m^{}_0`$, we have

$`g_s^{XY}`$ $`={\displaystyle \frac{1}{(n2rw)^2}}\left(1+{\displaystyle \frac{X}{(n2rw)^2}}\right)^{(N_b\stackrel{~}{N}_b)}`$ (120)
$`={\displaystyle \frac{a^{(N_b\stackrel{~}{N}_b)}}{(n2rw)^2}},`$
where we have put

$`a1+{\displaystyle \frac{X}{(n2rw)^2}},`$ (121)
which appears frequently in the following discussion. If we choose a large enough cell like $`N_c>>2`$, we obtain

$`g_s^{XY}={\displaystyle \frac{a^7}{(n2rw)^2}},`$ (122)
according to Eq. (112). Physically the factor $`a^7`$ with $`a>1`$ represents the exclusion effect of the bonds as shown in Fig. 4(b).

### C Evaluation of $`g_s^Z`$

The Gutzwiller factor for the $`z`$-component of $`𝑺_{\mathrm{}}𝑺_m`$ is calculated in the similar way. The lowest orders with respect to $`\chi `$ and $`\mathrm{\Delta }`$ give the similar results to (110) and (111), in which $`S_{\mathrm{}}^+S_m^{}_0`$ is replaced with $`S_{\mathrm{}}^zS_m^z_0`$. However, for $`S_{\mathrm{}}^zS_m^z_{i_0}`$, an additional contribution appears from the diagram in Fig. 5(a) which did not give contributions in $`S_{\mathrm{}}^+S_m^{}_{i_0}`$. This is the reason why the anisotropy for the Gutzwiller factor between $`xy`$ component and $`z`$ component appears, which is the most important feature of our extended Gutzwiller approximation.

Assuming that the site $`\mathrm{}`$ and $`m^{}`$ in Fig. 5(a) are in the sublattice 1, we have a contribution

$`S_{\mathrm{}}^zS_m^z\left\{\widehat{n}_m^{}(1\widehat{n}_m^{})+\widehat{n}_m^{}(1\widehat{n}_m^{})\right\}_c`$ (123)
$`=S_{\mathrm{}}^z_0S_m^z\widehat{n}_m^{}(1\widehat{n}_m^{})+S_m^z\widehat{n}_m^{}(1\widehat{n}_m^{})_c`$ (124)
$`={\displaystyle \frac{m}{2}}(P_{}P_{}+P_{}P_{})`$ (125)
$`=m^2X_2,`$ (126)
with $`X_2`$ defined in Eq. (18). Since the number of bonds connected to the site $`m`$ is

$`N_2{\displaystyle \frac{N_b\stackrel{~}{N}_b1}{2}},`$ (127)
the contribution to $`_{i_0}S_{\mathrm{}}^zS_m^z_{i_0}`$ turns out to be

$`2N_2m^2X_2(n2rw)^{N_c3},`$ (128)
where the contributions from the bonds connected to the site $`\mathrm{}`$ are also included.

Similarly the diagram in Fig. 5(b) gives a contribution

$`S_{\mathrm{}}^z\widehat{n}_{\mathrm{}^{}}(1\widehat{n}_{\mathrm{}^{}})+S_{\mathrm{}}^z\widehat{n}_{\mathrm{}^{}}(1\widehat{n}_{\mathrm{}^{}})_c`$ (129)
$`\times S_m^z\widehat{n}_m^{}(1\widehat{n}_m^{})+S_m^z\widehat{n}_m^{}(1\widehat{n}_m^{})_c`$ (130)
$`={\displaystyle \frac{1}{4}}(P_{}P_{}+P_{}P_{})^2`$ (131)
$`=m^2X_2^2.`$ (132)
By counting the number of possible combination of the bonds, we have

$`N_2^2m^2X_2^2(n2rw)^{N_c4}.`$ (133)
The higher order terms of $`X^j`$ are calculated as before, which become

$`S_{\mathrm{}}^zS_m^z={\displaystyle \frac{1}{W_0}}\{`$ $`S_{\mathrm{}}^zS_m^z_0(n2rw)^{N_c2}a^{\stackrel{~}{N}_b}`$ (134)
$``$ $`2N_2m^2X_2(n2rw)^{N_c3}a^{\stackrel{~}{N}_{b}^{}{}_{}{}^{}}`$ (135)
$``$ $`N_2^2m^2X_2^2(n2rw)^{N_c4}a^{\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime }}\},`$ (136)
with $`\stackrel{~}{N}_{b}^{}{}_{}{}^{}`$ and $`\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime }`$ being the numbers of the bonds which are not connected directly to the diagrams in Figs. 5(a) and 5(b), respectively. They are shown in Figs. 6(a) and (b). In the large enough cell, we have

$`\stackrel{~}{N}_{b}^{}{}_{}{}^{}`$ $`=N_b10,`$ (137)
$`\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime }`$ $`=N_b13.`$ (138)
Combining with $`W_0`$, we obtain

$`S_{\mathrm{}}^zS_m^z=`$ $`{\displaystyle \frac{a^{(N_b\stackrel{~}{N}_b)}}{(n2rw)^2}}\{S_{\mathrm{}}^zS_m^z_0`$ (142)
$`{\displaystyle \frac{2N_2}{n2rw}}m^2X_2a^{(\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{})}`$
$`{\displaystyle \frac{N_2^2}{(n2rw)^2}}m^2X_2^2a^{(\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime })}\}.`$
In this case, $`S_{\mathrm{}}^zS_m^z`$ is not proportional to $`S_{\mathrm{}}^zS_m^z_0`$. Since the Gutzwiller factor is the ratio between the two, we have some nontrivial contribution from the second and third terms on the r.h.s of (142). Using

$`S_{\mathrm{}}^zS_m^z_0=m^2{\displaystyle \frac{X_2}{4}},`$ (143)
the final expression for the Gutzwiller factor turns out to be

$`g_s^Z`$ $`=g_s^{XY}{\displaystyle \frac{1}{4m^2+X_2}}`$ (145)
$`\times \left[X_2+4m^2\left\{1+{\displaystyle \frac{N_2X_2}{n2rw}}a^{(\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{})}\right\}^2\right],`$
where we have assumed

$`\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime }=2(\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{}).`$ (146)
For a large enough cell we have $`N_2=3`$ and $`\stackrel{~}{N}_b\stackrel{~}{N}_{b}^{}{}_{}{}^{}=3`$, so that

$$g_s^Z=g_s^{XY}\frac{1}{4m^2+X_2}\left[X_2+4m^2\left\{1+\frac{3X_2}{n2rw}a^3\right\}^2\right].$$
(147)
The above formula is one of the most important results in this paper. Typical $`m`$-dependences have been shown in Fig. 1. Let us here check some limiting cases. When $`m=0`$, we have

$$g_s^{XY}=g_s^Z=4a^7,$$
(148)
which reproduces the original Gutzwiller approximation. For small $`m`$, we obtain

$$g_s^Z=g_s^{XY}\left\{1+4m^2\left(\frac{6}{n2rw}a^3+\frac{9X_2}{(n2rw)^2}a^6\right)\right\}.$$
(149)
It is apparent that $`g_s^Z`$ has an enhancement as a function of $`m`$, which gives the reasonable AF long-range order at half filling.

### D Physical meaning of the enhancement of $`g_s^Z`$

In the usual interpretation of the Gutzwiller approximation, a comparison is made between the probabilities of spin configurations in the wave functions with and without the projection operators. When we consider only two sites $`\mathrm{}`$ and $`m`$ for $`𝑺_{\mathrm{}}𝑺_m`$, the summation of the probabilities of spin configurations $`|,|,|`$ and $`|`$ in the wave function without the projection is given by $`P_{}+P_{}+P_{}+P_{}`$ which was calculated as $`(n2rw)^2+X`$ in Eq. (104). On the other hand, in the wave function with the projection, we have always one of the above four states at half filling, and thus the probability is equal to 1. The ratio of these probabilities gives the Gutzwiller factor

$`g_s^{XY}=g_s^Z={\displaystyle \frac{1}{(n2rw)^2+X}}.`$ (150)
However this Gutzwiller approximation does not give a reasonable answer compared with the VMC results as mentioned before.

For the physical understanding of the results in the previous subsections, it is necessary to go beyond the two sites, $`\mathrm{}`$ and $`m`$, because the enhancement of $`g_s^Z`$ comes from the diagrams in Fig. 5. For this purpose we consider the spin configurations on three sites as shown in Fig. 7. To calculate the probabilities of these configurations is straightforward. For example for Fig. 7(c) we have

$`\left(\begin{array}{ccc}& & \\ & & \end{array}\right)=2r(1w)P_{}2r^3(1w)^3,`$ (151)
where the second term on the r.h.s. appears to avoid the double counting of the zeroth order with respect to $`\chi `$ and $`\mathrm{\Delta }`$. We can see that the presence of the third site $`m^{}`$ enhances the probability of this configuration, $`\left(\begin{array}{ccc}& & \\ & & \end{array}\right)`$, because $`P_{}>P_{}=P_{}>P_{}`$ in the presence of AF order. As a result, the weight of the configuration $``$ on $`\mathrm{}`$ and $`m`$ sites is increased in the wave function. This causes the enhancement of $`g_s^Z`$.

To understand this effect more quantitatively, we calculate $`S_{\mathrm{}}^zS_m^z`$ directly from the configurations in Fig. 7 as follows:

$`4S_{\mathrm{}}^zS_m^z`$ (152)
$`={\displaystyle \frac{(\mathrm{a})+(\mathrm{b})(\mathrm{c})(\mathrm{d})(\mathrm{e})(\mathrm{f})+(\mathrm{g})+(\mathrm{h})}{(\mathrm{a})+(\mathrm{b})+(\mathrm{c})+(\mathrm{d})+(\mathrm{e})+(\mathrm{f})+(\mathrm{g})+(\mathrm{h})}}.`$ (153)
The denominator becomes

$`(n2rw)\{(n2rw)^2+2X\},`$ (154)
while the numerator becomes

$`(n2rw)(P_{}P_{}P_{}+P_{})`$ (155)
$`+2m(P_{}P_{}+P_{}P_{})+4m^2(n2rw)`$ (156)
$`=4(n2rw)S_{\mathrm{}}^zS_m^z_04m^2X_2.`$ (157)
Although the denominator is larger than the two-site case in Eq. (150), the second term in the numerator enhances $`S_{\mathrm{}}^zS_m^z`$, which is mainly from the contribution of the configuration 7(c). Using $`4S_{\mathrm{}}^zS_m^z_0=(4m^2+X_2)`$ and taking the ratio, we obtain

$`{\displaystyle \frac{S_{\mathrm{}}^zS_m^z}{S_{\mathrm{}}^zS_m^z_0}}`$ $`={\displaystyle \frac{1}{(n2rw)^2+2X}}{\displaystyle \frac{1}{4m^2+X_2}}`$ (159)
$`\times \left\{X_2+4m^2(1+{\displaystyle \frac{X_2}{n2rw}})\right\}.`$
This is essentially the result obtained in the previous subsection for $`g_s^Z`$. Since we have considered only the site $`m^{}`$ as the third term, $`N_2`$ in Eq. (145) is replaced with $`1`$ in this case. The first factor $`1/\{(n2rw)^2+2X\}`$ corresponds the exclusion effect which was represented by $`a^{(N_b\stackrel{~}{N}_b)}`$ in $`g_s^{XY}`$. From this analysis it is apparent that the enhancement of $`g_s^Z`$ as a function of $`m`$ is due to the increase of the probability of $``$ in the presence of the AF circumstances.

## V Less-than-half-filled case

### A Evaluation of the total weight $`W`$

In the presence of holes, the Gutzwiller factors are calculated similarly as in the half-filled case. In this subsection we obtain the total weight, $`W`$. Since the weight for the zero-hole sector $`W_0`$ has been calculated, we have to evaluate

$`W_1={\displaystyle \underset{i\mathrm{with}1\mathrm{hole}}{}}\omega _i,`$ (160)
and $`W_2`$ and so on.

In the lowest order of $`\chi `$ and $`\mathrm{\Delta }`$, we have

$`W_1^{(0)}=N_c(1r)(1w)(n2rw)^{N_c1},`$ (161)
where the factor $`N_c`$ comes from the possible position of the hole in a cell and the factor $`(1r)(1w)`$ comes from the expectation value of the hole position, $`(1\widehat{n}_i)(1\widehat{n}_i)_0`$. In the next order, by counting the number of bonds which contribute $`X`$, we obtain

$`N_c(1r)(1w)N_{1b}X(n2rw)^{N_c3},`$ (162)
where $`N_{1b}`$ means the number of bonds which are not connected the hole site as shown in Fig. 8(a). In the large enough cell we have

$`N_{1b}=N_b4.`$ (163)
Although the dotted bonds in Fig. 8(a) do not give $`X`$, they give alternative contributions to $`W_1`$ in the same order. This contribution is calculated as

$`(1\widehat{n}_i)(1\widehat{n}_i)\{\widehat{n}_j(1\widehat{n}_j)+\widehat{n}_j(1\widehat{n}_j)\}_0`$ (164)
$`=(n2rw)(1r)(1w)+Y,`$ (165)
with

$$Y=\delta (1+\delta )(\chi ^2\mathrm{\Delta }^2)4m^2(\chi ^2+\mathrm{\Delta }^2)2(\chi ^2+\mathrm{\Delta }^2)^2.$$
(166)
Here $`i`$ represents the site of the hole and $`j`$ is one of the nearest neighbor sites. Since the number of the dotted bonds in Fig. 8(a) is $`(N_bN_{1b})`$, the additional contribution to $`W_1`$ becomes

$`N_c(N_bN_{1b})Y(n2rw)^{N_c2}.`$ (167)
The higher order terms with respect to $`X`$ can be taken into account in the similar way as before and we obtain

$`W_1`$ $`=N_c(1r)(1w)(n2rw)^{N_c1}a^{N_{1b}}`$ (169)
$`+N_c(N_bN_{1b})Y(n2rw)^{N_c2}a^{\stackrel{~}{N}_b},`$
where $`\stackrel{~}{N}_b`$ appears because the diagram in Fig. 8(a) excludes the same number of bonds as in Fig. 5(a). We assume that

$`\stackrel{~}{N}_b=N_{1b}3,`$ (170)
which leads to

$`W_1`$ $`=N_c\left\{(1r)(1w)+{\displaystyle \frac{N_bN_{1b}}{n2rw}}Ya^3\right\}`$ (172)
$`\times (n2rw)^{N_c1}a^{N_{1b}}.`$
Strictly speaking, when the hole is located on the boundary of the cell, the weight should be different. However we make an approximation of large enough cell.

Continuing the similar arguments, we obtain

$`W_2`$ $`=_{N_c}C_2\left\{(1r)(1w)+{\displaystyle \frac{N_bN_{1b}}{n2rw}}Ya^3\right\}^2`$ (174)
$`\times (n2rw)^{N_c2}a^{N_{2b}},`$
for the two-hole sector, where $`N_{2b}`$ is the number of bonds which are not connected to the two sites of holes (Fig. 8(b)). We have neglected the case in which two holes come to the nearest neighbor sites. From Eq. (174) it is apparent that the $`j`$-hole sector, $`W_j`$, can be generally approximated as

$`W_j=_{N_c}C_jz^j(n2rw)^{N_cj}a^{N_{j,b}},`$ (175)
with

$`z=(1r)(1w)+{\displaystyle \frac{N_bN_{1b}}{n2rw}}Ya^3.`$ (176)
Finally approximating that $`N_{j,b}=N_b4j`$, the constraints which determine $`W`$ and $`p`$ (Eq. (89)) become

$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}{\displaystyle \frac{W_j}{W}}p^j`$ $`={\displaystyle \frac{1}{W}}(n2rw+pza^4)^{N_c}a^{N_b}=1,`$ (177)
$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}j{\displaystyle \frac{W_j}{W}}p^j`$ $`={\displaystyle \frac{N_c}{W}}pza^4(n2rw+pza^4)^{N_c1}a^{N_b}`$ (179)
$`=\delta N_c.`$
From these constraints we obtain

$`p`$ $`={\displaystyle \frac{\delta (n2rw)}{nz}}a^4,`$ (180)
$`W`$ $`=\left({\displaystyle \frac{n2rw}{n}}\right)^{N_c}a^{N_b},`$ (181)
for less-than-half-filling.

### B Evaluation of $`g_s^{XY}`$ and $`g_s^Z`$

In the similar way we estimate $`g_s^{XY}`$. In the case when there are $`j`$ holes in the central cell, we have

$`{\displaystyle \underset{i_0\mathrm{with}j\mathrm{holes}}{}}S_{\mathrm{}}^+S_m^{}_{i_0}`$ (182)
$`=_{N_c2}C_jz^j(n2rw)^{N_c2j}a^{\stackrel{~}{N}_b4j}S_{\mathrm{}}^+S_m^{}_0.`$ (183)
By summing up all the sectors with different number of holes, we have (according to the general formula (92)),

$`S_{\mathrm{}}^+S_m^{}`$ $`={\displaystyle \underset{j=0}{\overset{N_c2}{}}}{\displaystyle \frac{p^j}{W}}{\displaystyle \underset{i_0\mathrm{with}j\mathrm{holes}}{}}S_{\mathrm{}}^+S_m^{}_{i_0}`$ (186)
$`={\displaystyle \frac{1}{W}}(n2rw+pza^4)^{N_c2}a^{\stackrel{~}{N}_b}S_{\mathrm{}}^+S_m^{}_0`$
$`=\left({\displaystyle \frac{n}{n2rw}}\right)^2a^{(N_b\stackrel{~}{N}_b)}S_{\mathrm{}}^+S_m^{}_0.`$
From this we obtain the Gutzwiller factor

$`g_s^{XY}=\left({\displaystyle \frac{n}{n2rw}}\right)^2a^7.`$ (187)
This is a simple generalization of the results at half filling obtained in Eq. (122).

Although the estimation of $`g_s^Z`$ is essentially the same as in the half-filled case, there are some additional diagrams to be taken into account which are shown in Fig. 9. The detailed calculations are summarized in Appendix C. The final result is

$`g_s^Z`$ $`=g_s^{XY}{\displaystyle \frac{1}{4m^2+X_2}}`$ (189)
$`\times \left[X_2+4m^2\left\{1+{\displaystyle \frac{N_2X_2n}{n2rw}}\left(1{\displaystyle \frac{p}{2}}\right)a^3\right\}^2\right].`$
The quantity $`p`$ is given in Eq. (181). For $`g_s^Z`$, we use the lowest order approximation

$$p\frac{\delta (n2rw)}{nz}2\delta ,$$
(190)
since $`n2rw=1/2`$ and $`z=1/4`$ in the lowest order with respect to $`\delta ,\chi `$ and $`\mathrm{\Delta }`$. This gives our final expression for $`g_s^Z`$ given in Eq. (15).

### C Evaluation of $`g_t`$

As for the Gutzwiller factor for the hopping term, $`g_t`$, we have to calculate

$`(1\widehat{n}_i)c_i^{}c_j(1\widehat{n}_j)_{i_0}.`$ (191)
Since the hopping term needs at least one hole, the lowest sector is the one-hole sector. For this sector, we have

$`{\displaystyle \underset{i_0\mathrm{with}1\mathrm{hole}}{}}(1\widehat{n}_i)c_i^{}c_j(1\widehat{n}_j)_{i_0}`$ (192)
$`=(1\widehat{n}_i)c_i^{}c_j(1\widehat{n}_j)_0a^{\stackrel{~}{N}_b}(n2rw)^{N_c2}`$ (193)
$`=Ta^{\stackrel{~}{N}_b}(n2rw)^{N_c2},`$ (194)
where the expectation value of the correlated hopping, $`T`$, is calculated as

$`T`$ $`(1\widehat{n}_i)c_i^{}c_j(1\widehat{n}_j)_0`$ (197)
$`=(1r)(1w)\chi \chi ^3\chi \mathrm{\Delta }^2`$
$`=\chi \{(1r)(1w){\displaystyle \frac{X_2}{2}}\}.`$
It is easy to show that for the $`j+1`$-hole sector

$`{}_{N_c2}{}^{}C_{j}^{}Ta^{\stackrel{~}{N}_{j,b}}(n2rw)^{N_c2j}z^j.`$ (198)
Therefore summing up all the sectors, we obtain

$`c_i^{}c_j`$ $`={\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c2}{}}}{}_{N_c2}{}^{}C_{j}^{}Ta^{\stackrel{~}{N}_{j,b}}(n2rw)^{N_c2j}z^jp^{j+1}`$ (200)
$`=\left({\displaystyle \frac{n}{n2rw}}\right)^2Ta^{(N_b\stackrel{~}{N}_b)}p.`$
Since $`c_i^{}c_j_0=\chi `$, the Gutzwiller factor $`g_t=c_i^{}c_j/c_i^{}c_j_0`$ becomes

$`g_t`$ $`=\left({\displaystyle \frac{n}{n2rw}}\right)^2\{(1r)(1w){\displaystyle \frac{X_2}{2}}\}a^{(N_b\stackrel{~}{N}_b)}p`$ (203)
$`={\displaystyle \frac{n}{n2rw}}{\displaystyle \frac{\delta \{(1r)(1w)\frac{X_2}{2}\}}{z}}a^{(N_b\stackrel{~}{N}_b)+4},`$
where we have substituted the value of $`p`$ obtained in Eq. (181). The diagrams in Fig. 5 give contributions of the order of $`\delta ^2\chi ^2`$ or $`\delta m^2\chi ^2`$ so that they are neglected.

Finally comparing the definitions of $`X`$ and $`Y`$, we use an approximation

$`Y{\displaystyle \frac{X}{2}},`$ (204)
for the small $`\delta `$ case. Using this approximation, the quantity $`z`$ can be rewritten as

$`z`$ $`=(1r)(1w)\left\{1+{\displaystyle \frac{N_bN_{1b}}{(1r)(1w)(n2rw)}}Ya^3\right\}`$ (208)
$`=(1r)(1w)\left\{14(N_bN_{1b})X\right\}`$
$`=(1r)(1w)a^{(N_bN_{1b})}`$
$`=(1r)(1w)a^4.`$
Here we are considering $`\delta ,m`$ and $`\chi `$ as small quantities. Substituting this approximate $`z`$ into $`g_t`$ we have

$`g_t`$ $`={\displaystyle \frac{\delta n}{n2rw}}{\displaystyle \frac{(1r)(1w)\frac{X_2}{2}}{(1r)(1w)}}a^{(N_b\stackrel{~}{N}_b)+8}`$ (210)
$`={\displaystyle \frac{2\delta (1\delta )}{1\delta ^2+4m^2}}{\displaystyle \frac{(1+\delta )^24m^22X_2}{(1+\delta )^24m^2}}a.`$
This is our final expression for $`g_t`$ given in Eq. (20).

## VI Optimized variational state

In this section we calculate the energy of the variational state using the Gutzwiller factors obtained in the previous sections. We will show that the variational energies and the magnitudes of order parameters agree with the results in VMC simulations.

In the variational state $`|\psi =P_G|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V})`$ in Eq. (2), $`|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V})`$ is a Hartree-Fock type wave function with the $`d`$-wave SC and AF orders. It is expressed as

$`|\psi _0(\mathrm{\Delta }_d^\mathrm{V},\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V},\mu ^\mathrm{V})={\displaystyle \underset{k,s(=\pm )}{}}(u_k^{(s)}+v_k^{(s)}d_k^{(s)}d_k^{(s)})|0`$ (211)
$`={\displaystyle \underset{k,s}{}}u_k^{(s)}\mathrm{exp}\left[{\displaystyle \underset{k,s}{}}{\displaystyle \frac{v_k^{(s)}}{u_k^{(s)}}}d_k^{(s)}d_k^{(s)}\right]|0,`$ (212)
where

$`{\displaystyle \frac{v_k^{(\pm )}}{u_k^{(\pm )}}}`$ $`=`$ $`{\displaystyle \frac{\pm \mathrm{\Delta }_d^\mathrm{V}\eta _k}{\left(\pm E_k\mu ^\mathrm{V}\right)+\sqrt{\left(\pm E_k\mu ^\mathrm{V}\right)^2+(\mathrm{\Delta }_d^\mathrm{V}\eta _k)^2}}},`$ (213)
$`E_k`$ $`=`$ $`\sqrt{ϵ_k^2+\mathrm{\Delta }_{\mathrm{af}}^{\mathrm{V2}}},`$ (214)
$`ϵ_k=t\gamma _k`$, $`\gamma _k=2(\mathrm{cos}k_x+\mathrm{cos}k_y)`$ and $`\eta _k=2(\mathrm{cos}k_x\mathrm{cos}k_y)`$. The annihilation operators $`d_{k\sigma }^{(s)}`$ are related to the electron operators through the following unitary transformation,

$$\left(\begin{array}{c}d_{k\sigma }^{(+)}\\ d_{k\sigma }^{()}\end{array}\right)=\left(\begin{array}{cc}\alpha _{k\sigma }& \beta _{k\sigma }\\ \beta _{k\sigma }& \alpha _{k\sigma }\end{array}\right)\left(\begin{array}{c}c_{Ak\sigma }\\ c_{Bk\sigma }\end{array}\right),$$
(215)
with

$$\{\begin{array}{c}\alpha _{k\sigma }=\sqrt{\frac{1}{2}\left(1\frac{\sigma \mathrm{\Delta }_{\mathrm{af}}^\mathrm{V}}{E_k}\right)}\\ \beta _{k\sigma }=\sqrt{\frac{1}{2}\left(1+\frac{\sigma \mathrm{\Delta }_{\mathrm{af}}^\mathrm{V}}{E_k}\right)}\end{array}.$$
(216)
Here $`c_{Ak\sigma }(c_{Bk\sigma })`$ are annihilation operators of an electron on the A(B)-sublattice and $`\sigma `$ represent $``$(+1) and $``$(-1). The wave vector $`𝒌`$ is limited to half of the Brillouin zone where $`ϵ_k<0`$. We can confirm that $`|\psi _0`$ is a vacuum of the annihilation operators which diagonalize

$`{\displaystyle \underset{k}{}}`$ $`[{\displaystyle \underset{\sigma }{}}\{ϵ_k(c_{Ak\sigma }^{}c_{Bk\sigma }+h.c.)`$ (219)
$`(\mu ^\mathrm{V}+\sigma \mathrm{\Delta }_{\mathrm{af}}^\mathrm{V})c_{Ak\sigma }^{}c_{Ak\sigma }(\mu ^\mathrm{V}\sigma \mathrm{\Delta }_{\mathrm{af}}^\mathrm{V})c_{Bk\sigma }^{}c_{Bk\sigma }\}`$
$`\mathrm{\Delta }_d^\mathrm{V}\eta _k(c_{Ak}c_{Bk}+c_{Bk}c_{Ak}+h.c.)].`$
In order to clarify the correspondence to the mean-field theory, let us consider the effective Hamiltonian, $`\widehat{}_{\mathrm{eff}}`$ in Eq. (5). In $`\widehat{}_{\mathrm{eff}}`$ the parameter $`t`$ in $``$ is replaced with

$$t_{\mathrm{eff}}=g_tt,$$
(220)
and the exchange term is replaced with

$`J𝑺_i𝑺_j=g_s^{XY}J\left(S_i^xS_j^x+S_i^yS_j^y\right)+g_s^ZJS_i^zS_j^z,`$ (221)
where $`g_t`$, $`g_s^{XY}`$ and $`g_s^Z`$ are the Gutzwiller factors obtained in the previous sections. When we apply the mean-field theory to $`\widehat{}_{\mathrm{eff}}`$, we obtain the similar Hamiltonian as in Eq. (219) but with the replacements

$`t`$ $`t_{\mathrm{eff}}+J_{\mathrm{eff}}\chi ^\mathrm{V},`$ (222)
$`\mathrm{\Delta }_d^\mathrm{V}`$ $`J_{\mathrm{eff}}\mathrm{\Delta }^\mathrm{V},`$ (223)
$`\mathrm{\Delta }_{\mathrm{af}}^\mathrm{V}`$ $`2J_{\mathrm{eff}}^Zm^\mathrm{V},`$ (224)
with

$$J_{\mathrm{eff}}=\frac{1}{2}g_s^{XY}J+\frac{1}{4}g_s^ZJ,$$
(225)
$$J_{\mathrm{eff}}^Z=g_s^ZJ.$$
(226)
Note that in the usual mean-field theory the self-consistency equations give $`\chi ^\mathrm{V}=\chi ,\mathrm{\Delta }^\mathrm{V}=\mathrm{\Delta }`$ and $`m^\mathrm{V}=m`$, where $`\chi ,\mathrm{\Delta }`$ and $`m`$ are the expectation values $`\mathrm{\Delta }=c_i^{}c_j^{}_0`$, $`\chi =c_{i\sigma }^{}c_{j\sigma }_0`$, and $`m=\frac{1}{2}(1)^i\left(\widehat{n_i}_0\widehat{n_i}_0\right)`$ used in the previous sections. We will see shortly that $`\chi ^\mathrm{V}`$ and $`\chi `$ etc. are slightly different due to the dependence of the Gutzwiller factors on $`\chi ,\mathrm{\Delta }`$ and $`m`$.

Using the wave function $`|\psi _0`$, the expectation values become

$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{1}{8N}}{\displaystyle \underset{k,\pm }{}}{\displaystyle \frac{F_k\eta _k}{\sqrt{(\pm E_k\mu ^\mathrm{V})^2+F_k^2}}},`$ (227)
$`\chi `$ $`=`$ $`{\displaystyle \frac{1}{8N}}{\displaystyle \underset{k,\pm }{}}{\displaystyle \frac{(t_{\mathrm{eff}}+J_{\mathrm{eff}}\chi ^\mathrm{V})\gamma _k^2(\pm E_k\mu ^\mathrm{V})}{\pm E_k\sqrt{(\pm E_k\mu ^\mathrm{V})^2+F_k^2}}},`$ (228)
$`m`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{k,\pm }{}}{\displaystyle \frac{2J_{\mathrm{eff}}^Zm^\mathrm{V}(\pm E_k\mu ^\mathrm{V})}{\pm E_k\sqrt{(\pm E_k\mu ^\mathrm{V})^2+F_k^2}}},`$ (229)
$`n`$ $`=`$ $`1{\displaystyle \frac{1}{N}}{\displaystyle \underset{k,\pm }{}}{\displaystyle \frac{(\pm E_k\mu ^\mathrm{V})}{\sqrt{(\pm E_k\mu ^\mathrm{V})^2+F_k^2}}},`$ (230)
with

$`F_k`$ $`=`$ $`J_{\mathrm{eff}}\mathrm{\Delta }^\mathrm{V}\eta _k,`$ (231)
$`E_k`$ $`=`$ $`\sqrt{(t_{\mathrm{eff}}+J_{\mathrm{eff}}\chi ^\mathrm{V})^2\gamma _k^2+(2J_{\mathrm{eff}}^Zm^\mathrm{V})^2},`$ (232)
where $`N`$ is the number of sites and the summation over $`𝒌`$ is limited to half of the Brillouin zone.

Using these expectation values we obtain

$`E_{\mathrm{var}}`$ $`=`$ $`\widehat{}_{\mathrm{eff}}_0`$ (233)
$`=`$ $`8Nt_{\mathrm{eff}}\chi _{}4NJ_{\mathrm{eff}}\left(\mathrm{\Delta }^2+\chi ^2\right)`$ (235)
$`2NJ_{\mathrm{eff}}^Zm^2\mu (NnN_\mathrm{e}).`$
By taking the derivatives with respect to the variational parameters, $`\mathrm{\Delta }^\mathrm{V},\chi ^\mathrm{V}`$, and $`m^\mathrm{V}`$, we find the following self-consistency equations:

$`\mathrm{\Delta }^\mathrm{V}`$ $`=`$ $`\mathrm{\Delta }{\displaystyle \frac{1}{8NJ_{\mathrm{eff}}}}{\displaystyle \frac{H_{\mathrm{eff}}}{\mathrm{\Delta }}}_0,`$ (236)
$`\chi ^\mathrm{V}`$ $`=`$ $`\chi {\displaystyle \frac{1}{8NJ_{\mathrm{eff}}}}{\displaystyle \frac{H_{\mathrm{eff}}}{\chi }}_0,`$ (237)
$`m^\mathrm{V}`$ $`=`$ $`m{\displaystyle \frac{1}{4NJ_{\mathrm{eff}}^Z}}{\displaystyle \frac{H_{\mathrm{eff}}}{m}}_0,`$ (238)
$`\mu ^\mathrm{V}`$ $`=`$ $`\mu {\displaystyle \frac{1}{N}}{\displaystyle \frac{H_{\mathrm{eff}}}{n}}_0,`$ (239)
where the partial derivative of $`H_{\mathrm{eff}}`$ is applied to the Gutzwiller factors, $`g_s^{XY},g_s^Z`$ and $`g_t`$.

Figure 10 shows the self-consistent parameters $`\mathrm{\Delta },\chi `$ and $`m`$ satisfying (239) as a function of the doping $`\delta =1n`$ for $`J/t=0.3`$. Note that these are the expectation values in the wave function $`|\psi _0`$ without the projection. The dashed line in Fig. 10 represents the results when the AF order is suppressed, i.e., $`m`$ is fixed to zero. It is apparent that the presence of the AF order has a small effect on the expectation values $`\mathrm{\Delta }`$ and $`\chi `$. The coexistent state between the d-wave SC and AF order parameters is stabilize up to the doping rate $`\delta =0.1`$. This is consistent with the VMC simulations.

The actual expectation values in the wave function $`P_G|\psi _0`$ with projection are different from $`\mathrm{\Delta },\chi `$ and $`m`$. For example, the expectation value of $`c_{i\sigma }^{}c_{j\sigma }`$ is

$$\chi _{\mathrm{exp}}c_{i\sigma }^{}c_{j\sigma }=g_tc_{i\sigma }^{}c_{j\sigma }_0=g_t\chi .$$
(240)
In the similar way, the actual expectation values in $`P_G|\psi _0`$ are obtained using corresponding Gutzwiller factors. For $`m_{\mathrm{exp}}`$ we repeat the similar arguments as for $`g_s^Z`$ to obtain

$`m_{\mathrm{exp}}{\displaystyle \frac{1}{2}}n_in_i=g_mm,`$ (241)
$`g_m={\displaystyle \frac{n}{n2rw}}a^4\left[1+{\displaystyle \frac{(N_bN_{1b})X_2n}{n2rw}}\left(1{\displaystyle \frac{p}{2}}\right)a^3\right].`$ (242)
(243)
For the expectation value of $`c_i^{}c_j^{}`$, we have to take care of the fact that the number of holes changes in the cell. Here we use the average number of holes as an approximation. In the same level of approximations as for $`g_t`$ we obtain

$`\mathrm{\Delta }_{\mathrm{exp}}`$ $`c_i^{}c_j^{}=g_\mathrm{\Delta }\mathrm{\Delta },`$ (244)
$`g_\mathrm{\Delta }`$ $`={\displaystyle \frac{\delta n}{n2rw}}{\displaystyle \frac{(1r)^2+(1w)^2+X_2}{2(1r)(1w)}}a^{(N_b\stackrel{~}{N}_b)+8}`$ (246)
$`={\displaystyle \frac{2\delta (1\delta )}{1\delta ^2+4m^2}}{\displaystyle \frac{(1+\delta )^2+4m^2+2X_2}{(1+\delta )^24m^2}}a.`$
Figure 11 shows the actual expectation values, $`\mathrm{\Delta }_{\mathrm{exp}}`$ and $`m_{\mathrm{exp}}`$ as a function of the doping $`\delta =1n`$ for $`J/t=0.3`$.

## VII Summary and Discussion

In this paper we have developed a new approach for studying the effect of strong correlation or the Gutzwiller’s projection in the two-dimensional $`t`$-$`J`$ model. It is based on the extended Gutzwiller approximation, in which the effects of longer-range correlations are taken into account. These correlations play important roles for the interplay between the AF and d-wave SC. Let us summarize our main results and discuss related problems.

(1) Generally the expectation values with respect to the projected wave function are strongly renormalized due to the exclusion of double occupancies. In the slave-boson mean-field theory, this effect is taken into account by assuming the replacement $`t\delta t`$ and $`JJ`$. However we have shown that the renormalization is not so simple. First of all the renormalization factor for the exchange interaction is anisotropic, i.e., $`J_{\mathrm{eff}}^{XY}=g_s^{XY}JJ_{\mathrm{eff}}^Z=g_s^ZJ`$ in the presence of AF moment. Furthermore $`g_t,g_s^{XY}`$ and $`g_s^Z`$ have nonlinear dependences on the expectation values $`\mathrm{\Delta },\chi `$ and $`m`$. In this sense the extended Gutzwiller approximation is beyond the simple slave-boson mean-field theory. We think that essence of the strong correlation is contained in these Gutzwiller factors, because they stem solely from the projection.

The physical meanings of the Gutzwiller factors are clarified:

(i) The factor $`a^7`$ appearing in $`g_s^{XY}`$ and $`g_s^Z`$ represents the exclusion effect shown in Fig. 4(b). When we calculate $`𝑺_i𝑺_j`$, the surrounding six bonds cannot make contribution of nearest neighbor correlation $`X`$. As is evident in Fig. 1, the factor $`a^7`$ reduces the value of $`J_{\mathrm{eff}}`$ even for the case $`m=0`$.

(ii) The enhancement of $`g_s^Z`$ as a function of $`m`$ is the most important feature. This was observed in the VMC simulations. Here we have identified its origin as the effect of surrounding AF correlations by discussing the probabilities of spin configurations in Fig. 7. The enhancement of $`g_s^Z`$ is due to the increase of $``$ configuration caused by the AF circumstances.

(2) We have studied the projected variational states as an application of the new Gutzwiller approximation. It is shown that our scheme reproduces the results in the VMC simulations. The enhancement of $`g_s^Z`$ as a function of $`m`$ is essential for reproducing the coexistent state of the AF and d-wave SC orders found in VMC.

(3) At half filling ($`\delta =0`$), there are several interesting features. In the original Gutzwiller approximation, the AF state at half filling was not obtained which was unphysical. However in the present Gutzwiller approximation, the AF state with d-wave SC order parameter becomes the most optimized state as obtained in VMC. Note that the expectation value of the SC order parameter, $`\mathrm{\Delta }_{\mathrm{exp}}`$, is zero because of the Gutzwiller factor $`g_\mathrm{\Delta }=0`$ at half filling, although $`\mathrm{\Delta }`$ and the variational parameter $`\mathrm{\Delta }^V`$ are nonzero. The self-consistent values which we obtain are (Fig. 10)

$`\mu ^V`$ $`=0,`$ (247)
$`\mathrm{\Delta }^V`$ $`=\chi ^V=0.080,`$ (248)
$`\mathrm{\Delta }`$ $`=\chi =0.163,`$ (249)
$`m^V`$ $`=0.017,`$ (250)
$`m`$ $`=0.147.`$ (251)
The actual expectation value with the projection, $`m_{\mathrm{exp}}`$, is (see Fig. 11)

$$m_{\mathrm{exp}}=g_mm=0.373.$$
(252)
This is a reasonable magnitude and close to the Monte Carlo result for the Heisenberg model ($`m_{\mathrm{exp}}=0.31\pm 0.02`$).

The relations $`\mathrm{\Delta }^V=\chi ^V`$ and $`\mathrm{\Delta }=\chi `$ are the manifestation of the SU(2) symmetry at half filling. Actually the self-consistency equations for $`\chi `$ and $`\mathrm{\Delta }`$ become the same at $`\delta =0`$. It is worth while noting here that, due to the SU(2) symmetry, the coexistent state of the AF and d-wave SC is equivalent to the $`\pi `$-flux state with AF long range order discussed by Hsu.

The self-consistent value, $`m=0.147`$ in Eq. (251), is smaller than the result $`m=1/2`$ in the simple mean-field theory. Instead it is slightly larger than the value at which the Gutzwiller factor $`g_s^Z`$ has a maximum, as shown in Fig. 1. If $`m`$ is increased from the self-consistent value, $`g_s^Z`$ decreases. This indicates that $`m`$ is determined so as to minimize the exchange energy $`2NJ_{\mathrm{eff}}^Zm^2=2Ng_s^ZJm^2`$ by optimizing the energy gain due to the long-range order and the energy loss due to the reduction of $`g_s^Z`$. This situation is, in some sense, similar to the quantum fluctuations of the Heisenberg spin system discussed in the spin-wave theory: The spin fluctuation, which reduces $`m`$ from the mean-field value $`1/2`$, leads to the enhancement of the Gutzwiller’s renormalization factors to gain the energy.

(4) For less-than-half-filled case, the difference between the variational parameter $`\mathrm{\Delta }^V`$ and the expectation value $`\mathrm{\Delta }_{\mathrm{exp}}`$ is remarkable: $`\mathrm{\Delta }^V`$ is finite and moreover it increases as the doping rate $`\delta `$ decreases, while $`\mathrm{\Delta }_{\mathrm{exp}}`$ is proportional to $`\delta `$ near half filling. This is because the Gutzwiller factor, $`g_\mathrm{\Delta }`$, is proportional to $`\delta `$ (see section VI). We interpret that $`\mathrm{\Delta }^V`$ is the BCS-type energy gap observed in scanning tunnel spectroscopy or in break junctions, because it is the parameter embedded in the wave function even at half filling. The excitation spectra will have a large energy gap corresponding to $`\mathrm{\Delta }^V`$, while the true long-range order, $`\mathrm{\Delta }_{\mathrm{exp}}`$, is reduced due to the projection or the strong correlation. The increase of $`\mathrm{\Delta }^V`$ in decreasing $`\delta `$ is consistent with the dependence observed experimentally.

(5) We discuss here the relation to the SO(5) theory. In our formulation, the combination $`m^2+\mathrm{\Delta }^2`$ appears frequently. If the SO(5) symmetry is exact in the $`t`$-$`J`$ model, the free energy will have a systematic dependence such as

$$F(\mathrm{\Delta }^2+m^2)\mathrm{or}F(\mathrm{\Delta }_{\mathrm{exp}}^2+m_{\mathrm{exp}}^2).$$
(253)
For example, the numerator in $`a`$ has a combination

$$16(m^2+\mathrm{\Delta }^2+\chi ^2).$$
(254)
Also the factor giving the enhancement of $`g_s^Z`$ contains

$$4m^2+X_2=4m^2+2\mathrm{\Delta }^2+2\chi ^2.$$
(255)
Although these combinations remind us of the SO(5) symmetry, our Gutzwiller approximation does not show exact symmetry.

Nevertheless a tendency similar to the SO(5) prediction can be seen from the Gutzwiller factors in Fig. 1. Comparing $`g_s^{XY}`$ and $`g_s^Z`$ for the cases with $`\mathrm{\Delta }=0.02`$ and $`0.18`$, we find that $`g_s^{XY}`$ and $`g_s^Z`$ are larger for smaller $`\mathrm{\Delta }`$. This is mainly from the exclusion effect of $`a`$, because even at $`m=0`$ this effect is observed. Therefore if we consider a situation where the d-wave SC order parameter is suppressed, then the Gutzwiller factor is enhanced causing the increase of AF moment. This gives the similar phenomena predicted in the SO(5) theory.

(6) One advantage of the present theory is that it is easily applied to the inhomogeneous cases, such as the stripe state, vortex cores, and magnetic states around nonmagnetic or magnetic impurities, where the interplay between the AF and d-wave SC correlations plays an important role. Since our Gutzwiller approximation gives a reasonable estimate of the variational energies in the presence of AF correlations, a reliable analytic formulation can be given.

The simplest way of applying the present scheme to inhomogeneous problems is to assume that $`g_s^{XY},g_s^Z`$ and $`g_t`$ for each bond $`i,j`$ are determined locally from the expectation values $`\mathrm{\Delta }_{ij},\chi _{ij}`$ and $`m_i,m_j`$ for the bond. In this case, if the d-wave order parameter is reduced around vortex cores, impurities or stripes, then the Gutzwiller factors $`g_s^{XY}`$ and $`g_s^Z`$ are enhanced locally as expected from Fig. 1. This effect causes the local development of AF correlations, which can be observed experimentally. From these viewpoint, preliminary calculations for the vortex cores and stripe states have been published elsewhere.

Of course the VMC simulations give more accurate evaluation of the variational energies, if they are used for the inhomogeneous systems like stripe states. However there are some difficulties in applying the VMC. Firstly we have to treat a fairly large unit cell to study the slowly varying order parameters especially near half filling. For example, the incommensurability in the stripe state is close to $`(\pi ,\pi )`$ so that the period of the stripe pattern becomes fairly long. In this case the VMC simulations become difficult. Furthermore the choice of the functional form of the trial state in VMC is restricted, since only the small number of variational parameters can be used practically. On the contrary, our scheme based on the mean-field-type Gutzwiller approximation can be used in a fairly large system sizes. Moreover the order parameters on all the bonds, $`\mathrm{\Delta }_{ij},\chi _{ij}`$ and $`m_i`$ can be optimized in the similar sense to unrestricted Hartree-Fock theory. Therefore we can search for the microscopically optimized variational states in our scheme.

The authors wish to thank T. M. Rice, H. Fukuyama, Y. Tanaka, H. Tsuchiura and M. Sigrist for useful discussions. They also thank C.-M. Ho for critical reading of the manuscript. This work is supported in part by a Grant-in-Aid of of the Ministry of Education, Science, Sports and Culture.

## A Derivation of the expectation value of an operator

The evaluation of Eq. (90) is slightly complicated. The constraints for $`\{N_i^{}\}`$ are

$`{\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i^{}`$ $`={\displaystyle \frac{N_e}{2}}n_{i_0},{\displaystyle \underset{i=1}{\overset{K}{}}}n_iN_i^{}={\displaystyle \frac{N_e}{2}}n_{i_0},`$ (A1)
$`{\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}N_i^{}`$ $`=NN_en_{\mathrm{h}i_0},`$ (A2)
instead of (75) because $`N_i^{}`$ represents the number of cells of the $`i`$-th state except for the central cell.

In the same way as in the denominator, the largest term is given by

$`\overline{N_i^{}}={\displaystyle \frac{N}{N_c}}{\displaystyle \frac{\omega _i}{W^{}}}(p^{})^{n_{\mathrm{h}i}},`$ (A3)
with slightly different values of $`W^{}`$ and $`p^{}`$ because of the difference of the constraints. From the constraints (A1), we can see that $`\overline{N_i^{}}`$ satisfies the relations

$$\underset{i=1}{\overset{K}{}}\overline{N_i^{}}=\frac{N}{N_c}1,\underset{i=1}{\overset{K}{}}n_{\mathrm{h}i}\overline{N_i^{}}=NN_en_{\mathrm{h}i_0}.$$
(A4)
Therefore if we define the difference $`\mathrm{\Delta }\overline{N_i}=\overline{N_i^{}}\overline{N_i}`$, we have important relations

$`{\displaystyle \underset{i=1}{\overset{K}{}}}\mathrm{\Delta }\overline{N_i}=1,{\displaystyle \underset{i=1}{\overset{K}{}}}n_{\mathrm{h}i}\mathrm{\Delta }\overline{N_i}=n_{\mathrm{h}i_0}.`$ (A5)
By use of these, the ratio between the numerator and the denominator in Eq. (92) is calculated as follows:

$`\widehat{𝒪}`$ $`={\displaystyle \frac{\psi _0|P_G\widehat{𝒪}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}`$ (A11)
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{1}{\left(\frac{N}{N_c}\right)}}{\displaystyle \underset{i=1}{\overset{K}{}}}{\displaystyle \frac{\overline{N_i}!}{\overline{N_i^{}}!}}{\displaystyle \underset{i=1}{\overset{K}{}}}\omega _i^{\overline{N_i^{}}\overline{N_i}}\widehat{𝒪}_{i_0}`$
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{N_c}{N}}{\displaystyle \underset{i=1}{\overset{K}{}}}\left({\displaystyle \frac{\overline{N_i}}{\omega _i}}\right)^{\mathrm{\Delta }\overline{N_i}}\widehat{𝒪}_{i_0}`$
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{N_c}{N}}{\displaystyle \underset{i=1}{\overset{K}{}}}\left({\displaystyle \frac{N}{N_c}}{\displaystyle \frac{p^{n_{\mathrm{h}i}}}{W}}\right)^{\mathrm{\Delta }\overline{N_i}}\widehat{𝒪}_{i_0}`$
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{N_c}{N}}\left({\displaystyle \frac{N}{N_c}}{\displaystyle \frac{1}{W}}\right)^{{\scriptscriptstyle \mathrm{\Delta }\overline{N_i}}}\times p^{{\scriptscriptstyle n_{\mathrm{h}i}\mathrm{\Delta }\overline{N_i}}}\widehat{𝒪}_{i_0}`$
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{p^{n_{\mathrm{h}i_0}}}{W}}\widehat{𝒪}_{i_0}.`$
This gives Eq. (92).

## B Reproduction of the original Gutzwiller approximation

In this appendix we show that the original Gutzwiller approximation described in II is reproduced if we assume only the site-diagonal expectation values in the generalized formulation in section III.

First we calculate the weight $`\omega _i`$ of the $`i`$-th state which has $`n_i`$ up-spin electrons, $`n_i`$ down-spin electrons and $`n_{\mathrm{h}i}`$ holes. If we take only the site-diagonal expectation values, we have

$$\omega _i=[r(1w)]^{n_{\mathrm{right}}}[w(1r)]^{n_{\mathrm{wrong}}}[(1r)(1w)]^{n_{\mathrm{h}i}},$$
(B1)
where $`n_{\mathrm{right}}`$ ($`n_{\mathrm{wrong}}`$) means the number of sites where the right (wrong) spices of spin direction is located depending on the sublattice 1 and 2. Then the total weight in the subgroup with $`j`$ holes (Eq. (87)) can be calculated exactly as

$`W_j`$ $`={\displaystyle \underset{i\mathrm{with}j\mathrm{holes}}{}}\omega _i`$ (B4)
$`=_{N_c}C_j[r(1w)+w(1r)]^{N_cj}\left[(1r)(1w)\right]^j`$
$`=_{N_c}C_j(n2rw)^{N_cj}(1r)^j(1w)^j,`$
where $`{}_{N_c}{}^{}C_{j}^{}`$ is the number of choices of the positions of $`j`$ holes.

Using these $`W_j`$, the constraints (89) have simple forms as

$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}{\displaystyle \frac{W_j}{W}}p^j={\displaystyle \frac{1}{W}}\{n2rw+(1r)(1w)p\}^{N_c}=1,`$ (B5)
$`{\displaystyle \underset{j=0}{\overset{N_c}{}}}j{\displaystyle \frac{W_j}{W}}p^j={\displaystyle \frac{N_c}{W}}(1r)(1w)p`$ (B6)
$`\times \{n2rw+(1r)(1w)p\}^{N_c1}=\delta N_c.`$ (B7)
These can be solved easily to give

$`p={\displaystyle \frac{\delta (n2rw)}{n(1r)(1w)}},`$ (B8)
$`W=\left({\displaystyle \frac{n2rw}{n}}\right)^{N_c}.`$ (B9)
Finally the expectation value for $`S_{\mathrm{}}^+S_m^{}`$, for example, becomes

$`{\displaystyle \frac{\psi _0|P_GS_{\mathrm{}}^+S_m^{}P_G|\psi _0}{\psi _0|P_GP_G|\psi _0}}`$ (B10)
$`={\displaystyle \underset{i_0}{}}{\displaystyle \frac{p^{n_{\mathrm{h}i_0}}}{W}}S_{\mathrm{}}^+S_m^{}_{i_0}`$ (B11)
$`={\displaystyle \underset{j=0}{\overset{N_c2}{}}}{\displaystyle \frac{p^j}{W}}_{N_c2}C_j(n2rw)^{N_c2j}(1r)^j(1w)^jS_{\mathrm{}}^+S_m^{}_0`$ (B12)
$`={\displaystyle \frac{1}{W}}\{n2rw+(1r)(1w)p\}^{N_c2}S_{\mathrm{}}^+S_m^{}_0`$ (B13)
$`={\displaystyle \frac{n^2}{(n2rw)^2}}S_{\mathrm{}}^+S_m^{}_0,`$ (B14)
which is exactly the same as the results in the original Gutzwiller approximation in section II.B, Eq. (65).

## C Evaluation of $`g_s^Z`$ for the less-than-half-filled case

In the similar way to $`g_s^{XY}`$, we evaluate $`g_s^Z`$. In the zero-hole sector, we have

$`{\displaystyle \underset{i_0\mathrm{with}0\mathrm{holes}}{}}S_{\mathrm{}}^zS_m^z_{i_0}`$ $`=S_{\mathrm{}}^zS_m^z_0(n2rw)^{N_c2}a^{\stackrel{~}{N}_b}`$ (C1)
$``$ $`2N_2m^2X_2(n2rw)^{N_c3}a^{\stackrel{~}{N}_b3}`$ (C2)
$``$ $`N_2^2m^2X_2^2(n2rw)^{N_c4}a^{\stackrel{~}{N}_b6}.`$ (C4)
For the one-hole sector, there are five contributions from the diagrams in Fig. 9. The first three diagrams are simple extensions of the zero-hole sector. However Figs. 9(d) and 9(e) are new-type contributions due to the presence of a hole. For the diagrams in Figs. 9(d) and 9(e) we calculate

$`\left(\begin{array}{ccc}& & \\ & & \end{array}\right)`$ $`=S_{\mathrm{}}^zS_m^z(1\widehat{n}_m^{})(1\widehat{n}_m^{})_c`$ (C6)
$`={\displaystyle \frac{1}{2}}m^2X_2,`$
$`\left(\begin{array}{ccc}& & \\ & & \\ & & \end{array}\right)`$ $`=\{\widehat{n}_{\mathrm{}^{}}(1\widehat{n}_{\mathrm{}^{}})+\widehat{n}_{\mathrm{}^{}}(1\widehat{n}_{\mathrm{}^{}})\}S_{\mathrm{}}^zS_m^z`$ (C9)
$`\times (1\widehat{n}_m^{})(1\widehat{n}_m^{})_c`$
$`={\displaystyle \frac{1}{2}}m^2X_2^2.`$
Using these expectation values, we obtain

$`{\displaystyle \underset{i_0\mathrm{with}1\mathrm{holes}}{}}S_{\mathrm{}}^zS_m^z_{i_0}`$ (C10)
$`=_{N_c2}C_1S_{\mathrm{}}^zS_m^z_0z(n2rw)^{N_c3}a^{\stackrel{~}{N}_{1b}}`$ (C11)
$`_{N_c3}C_12N_2m^2X_2z(n2rw)^{N_c4}a^{\stackrel{~}{N}_{1b}3}`$ (C12)
$`_{N_c4}C_1N_2^2m^2X_2^2z(n2rw)^{N_c5}a^{\stackrel{~}{N}_{1b}6}`$ (C13)
$`+N_2m^2X_2(n2rw)^{N_c3}a^{\stackrel{~}{N}_b3}`$ (C14)
$`+N_2^2m^2X_2^2(n2rw)^{N_c4}a^{\stackrel{~}{N}_b6}.`$ (C15)
In the two hole sector, we need to calculate a diagram in Fig. 9(f) which is $`\frac{1}{4}m^2X_2^2`$. Taking account of these diagrams and counting the higher order terms with respect to $`X`$, we approximate as

$`{\displaystyle \underset{i_0\mathrm{with}2\mathrm{holes}}{}}S_{\mathrm{}}^zS_m^z_{i_0}`$ (C16)
$`=_{N_c2}C_2S_{\mathrm{}}^zS_m^z_0z^2(n2rw)^{N_c4}a^{\stackrel{~}{N}_{2b}}`$ (C17)
$`_{N_c3}C_22N_2m^2X_2z^2(n2rw)^{N_c5}a^{\stackrel{~}{N}_{2b}3}`$ (C18)
$`_{N_c4}C_2N_2^2m^2X_2^2z^2(n2rw)^{N_c6}a^{\stackrel{~}{N}_{2b}6}`$ (C19)
$`+_{N_c3}C_1N_2m^2X_2z(n2rw)^{N_c4}a^{\stackrel{~}{N}_{1b}3}`$ (C20)
$`+_{N_c4}C_1N_2^2m^2X_2^2z(n2rw)^{N_c5}a^{\stackrel{~}{N}_{1b}6}`$ (C21)
$`{\displaystyle \frac{N_2^2}{4}}m^2X_2^2(n2rw)^{N_c4}a^{\stackrel{~}{N}_b6}.`$ (C22)
Finally generalization to higher order contributions leads to

$`S_{\mathrm{}}^zS_m^z={\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c}{}}}p^j{\displaystyle \underset{i_0\mathrm{with}j\mathrm{holes}}{}}S_{\mathrm{}}^zS_m^z_{i_0}`$ (C23)
$`={\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c2}{}}}p_{N_c2}^jC_jS_{\mathrm{}}^zS_m^z_0z^j(n2rw)^{N_c2j}a^{\stackrel{~}{N}_{j,b}}`$ (C24)
$`{\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c3}{}}}p_{N_c3}^jC_j2N_2m^2X_2z^j(n2rw)^{N_c3j}a^{\stackrel{~}{N}_{j,b}3}`$ (C25)
$`{\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c4}{}}}p_{N_c4}^jC_jN_2^2m^2X_2^2z^j(n2rw)^{N_c4j}a^{\stackrel{~}{N}_{j,b}6}`$ (C26)
$`+{\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c3}{}}}p_{N_c3}^{j+1}C_jN_2m^2X_2z^j(n2rw)^{N_c3j}a^{\stackrel{~}{N}_{j,b}3}`$ (C27)
$`+{\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c4}{}}}p_{N_c4}^{j+1}C_jN_2^2m^2X_2^2z^j(n2rw)^{N_c4j}a^{\stackrel{~}{N}_{j,b}6}`$ (C28)
$`{\displaystyle \frac{1}{W}}{\displaystyle \underset{j=0}{\overset{N_c4}{}}}p_{N_c4}^{j+2}C_j{\displaystyle \frac{N_2^2}{4}}m^2X_2^2z^j(n2rw)^{N_c4j}a^{\stackrel{~}{N}_{j,b}6}`$ (C29)
$`=\left({\displaystyle \frac{n}{n2rw}}\right)^2S_{\mathrm{}}^zS_m^z_0a^{(N_b\stackrel{~}{N}_b)}`$ (C30)
$`\left({\displaystyle \frac{n}{n2rw}}\right)^32N_2m^2X_2\left(1{\displaystyle \frac{p}{2}}\right)a^{(N_b\stackrel{~}{N}_{b}^{}{}_{}{}^{})}`$ (C31)
$`\left({\displaystyle \frac{n}{n2rw}}\right)^4N_2^2m^2X_2^2\left(1{\displaystyle \frac{p}{2}}\right)^2a^{(N_b\stackrel{~}{N}_{b}^{}{}_{}{}^{\prime \prime })}`$ (C32)
$`=\left({\displaystyle \frac{n}{n2rw}}\right)^2a^{(N_b\stackrel{~}{N}_b)}`$ (C33)
$`\times \left[{\displaystyle \frac{X_2}{4}}m^2\left\{1+{\displaystyle \frac{N_2X_2n}{n2rw}}\left(1{\displaystyle \frac{p}{2}}\right)a^3\right\}^2\right].`$ (C34)
Thus the Gutzwiller factor is

$$g_s^Z=g_s^{XY}\frac{1}{4m^2+X_2}[X_2+4m^2\{1+\frac{N_2X_2n}{n2rw}(1\frac{p}{2})a^3)^2].$$
(C35)