\[\].

Model and Method
================

Our model is the 2D Hubbard model defined by $$H =
-t \sum_{<jl>,\sigma} (c_{j\sigma}^{\dag}c_{l\sigma} + 
                               \mbox{H.c.})
+ U \sum_j c_{j\uparrow}^{\dag} c_{j\uparrow}
                                  c_{j\downarrow}^{\dag} c_{j\downarrow},
                             \label{eq:H}$$ where $c_{j\sigma}^{\dag}$ ($c_{j\sigma}$) is the creation (annihilation) operator of an electron with spin $\sigma$ at the $j$th site; the sites form a rectangular lattice; $t$ is the transfer energy between the nearest-neighbor (n.n.) sites; $t$ is our energy unit; $<jl>$ denotes summation over all the n.n. bonds. $U$ is the on-site Coulomb energy. In this report, we also study the effect of $t'$ between n.n.n. sites by including $$H_{nnn} = -t' \sum_{<<jl>>,\sigma} (c_{j\sigma}^{\dag}c_{l\sigma} + 
                               \mbox{H.c.}) \label{eq:H2}$$ in the Hamiltonian; in the above equation $<<jl>>$ means summation over the n.n.n. pairs.

Our trial wave function is a Gutzwiller-projected BCS-type wave function defined as \[\]: $$\begin{aligned}
\Psi_{\mbox{\scriptsize s}} &=& P_{N_e} P_{\mbox{\scriptsize G}}
                                      \psi_{\mbox{\scriptsize BCS}}, \\
\psi_{\mbox{\scriptsize BCS}} &=& \prod_{{\mib k}} (u_{\mib k} + v_{\mib k}
        c_{{\mib k}\uparrow}^{\dag}c_{-{\mib k}\downarrow}^{\dag})|0>,\end{aligned}$$ where $P_{\mbox{\scriptsize G}}$ is the Gutzwiller projection operator given by $$P_{\mbox{\scriptsize G}} = \prod_j [1 - (1-g)n_{j\uparrow}n_{j\downarrow}];$$ $g$ is a variational parameter in the range from 0 to unity and $j$ labels a site in the real space. $P_{N_e}$ is a projection operator which extracts 