Correlated-Electron Systems and High-Temperature Superconductivity ()
1. Introduction
The effect of the strong correlation between electrons is important for many quantum critical phenomena, such as unconventional superconductivity (SC) and the metalinsulator transition. Typical correlated electron systems are high-temperature superconductors [1-5], heavy fermions [6-9] and organic conductors [10]. The phase diagram for the typical high-Tc cuprates is shown in Figure 1 [9]. It has a characteristics that the region of antiferromagnetic order exists at low carrier concentrations and the superconducting phase is adjacent to the antiferromagnetism.
In the low-carrier region shown in Figure 2, there is the anomalous metallic region where the susceptibility and 
show a peak above Tc suggesting an existence of the pseudogap. To clarify an origin of the anomalous metallic behaviors is also a subject attracting many physicists as a challenging problem.
It has been established that the Cooper pairs of high-Tc cuprates have the 
-wave symmetry in the hole-doped materials [11,12]. Several evidences of 
-wave pairing
Figure 1. Phase diagram delineating the regions of superconductivity and antiferromagnetic ordering of the Cu2+ ions for the hole-doped La2–xSrxCuO4 and electron-doped Nd2–xCexCuO4–y systems.
Figure 2. Phase diagram showing the regions of non-Fermi liquid and pseudogap metal for the hole-doped case. The boundaries indicated in the figure are not confirmed yet.
symmetry were provided for the electron-doped cuprates Nd2–xCexCuO4 [13-15]. Thus it is expected that the superconductivity of electronic origin is a candidate for the high-Tc superconductivity. We can also expect that the origin of 
-wave superconductivity lies in the onsite Coulomb interaction of the two-dimensional Hubbard model.
The antiferromagnetism should also be examined in correlated electron systems. The undoped oxide compounds exhibit rich structures of antiferromagnetic correlations over a wide range of temperature that are described by the two-dimensional quantum antiferromagnetism [16-18]. A small number of holes introduced by doping are responsible for the disappearance of long-range antiferromagnetic order [19-24].
Recent neutron scattering experiments have suggested an existence of incommensurate ground states with modulation vectors given by 
and 
(or 
and 
) where 
denotes the hole-doping ratio [25]. We can expect that the incommensurate correlations are induced by holes doped into the Cu-O plane in the underdoped region. A checkerboard-like charge-density modulation with a roughly 
period has also been observed by scanning tunneling microscopy experiments in Bi2212 and Na-CCOC compounds.
Recently the mechanism of superconductivity in hightemperature superconductors has been extensively studied using various two-dimensional (2D) models of electronic interactions. Among them the 2D Hubbard model [26] is the simplest and most fundamental model. This model has been studied intensively using numerical tools, such as the Quantum Monte Carlo method [27-42], and the variational Monte Carlo method [24,43-50]. The two-leg ladder Hubbard model was also investigated with respect to the mechanism of high-temperature superconductivity [51-59].
Since the discovery of cuprate high-temperature superconductors, many researchers have tried to explain the occurrence of superconductivity of these materials in terms of the two-dimensional (2D) Hubbard model. However, it remains matter of considerable controversial as to whether the 2D Hubbard model accounts for the properties of high-temperature cuprate superconductors. This is because the membership of the the two-dimensional Hubbard model in the category of strongly correlated systems is a considerable barrier to progress on this problem. The quest for the existence of a superconducting transition in the 2D Hubbard model is a long-standing problem in correlated-electron physics, and has been the subject of intensive study [35,36,38,46, 60,61]. In particular, the results of quantum Monte Carlo methods, which are believed to be exact unbiased methods, have failed to show the existence of superconductivity in this model [38,61].
In the weak coupling limit we can answer this question. We can obtain the superconducting order parameter of the Hubbard model in the limit of small
, that is given by [62-66]
(1)
where 
is the strength of the on-site Coulomb interaction and the exponent 
is determined by solving the gap equation. Thus the existence of the superconducting gap is guaranteed by the weak coupling theory although 
is extremely small because of the exponential behavior given above. 
indicates the strength of superconductivity. In the intermediate or large coupling region, we must study it beyond the perturbation theory.
We investigate the ground state of the Hubbard model by employing the variational Monte Carlo method. In the region
, the finite superconducting gap is obtained by using the quantum variational Monte Carlo method. The superconducting condensation energy obtained by adopting the Gutzwiller ansatz is in reasonable agreement with the condensation energy derived for YBa2Cu3O7. We have further investigated the stability of striped and checkerboard states in the under-doped region. Holes doped in a half-filled square lattice lead to an incommensurate spin and charge density wave. The relationship of the hole density x and incommensurability
, 
, is satisfied in the lower doping region. This is consistent with the results by neutron scattering measurements. To examine the stability of a 
checkerboard state, we have performed a variational Monte Carlo simulation on a two-dimensional 
Hubbard model with a Bi-2212 type band structure. We found that the 
period checkerboard checkerboard spin modulation that is characterized by multi 
vectors is stabilized.
Further investigation has been performed by using the quantum Monte Carlo method which is a numerical method that can be used to simulate the behavior of correlated electron systems. This method is believed to be an exact unbiased method. We compute pair correlation functions to examine a possibility of superconductivity.
The Quantum Monte Carlo (QMC) method is a numerical method employed to simulate the behavior of correlated electron systems. It is well known, however, that there are significant issues associated with the application to the QMC. First, the standard Metropolis (or heat bath) algorithm is associated with the negative sign problem. Second, the convergence of the trial wave function is sometimes not monotonic, and further, is sometimes slow. In past studies, workers have investigated the possibility of eliminating the negative sign problem [37,38,40-42]. We present the results obtained by a method, quantum Monte Carlo diagonalization, without the negative sign problem.
2. Hubbard Hamiltonian
The Hubbard Hamiltonian is
(2)
where 
and 
denote the creation and annihilation operators of electrons, respectively, and 
is the number operator. The second term represents the on-site Coulomb interaction which acts when the two electrons occupy the same site. The numbers of lattice sites and electrons are denoted as 
and
, respectively. The electron density is
.
In the non-interacting limit
, the Hamiltonian is easily diagonalized in terms of the Fourier transformation. In the ground state each energy level is occupied by electrons up to the Fermi energy. In the other limit
, each site is occupied by the upor down-spin electron, or is empty. The non-zero 
induces the movement of electrons that leads to a metallic state id
. The ground state is probably insulating at half-filling 
if 
is sufficiently large.
If 
are non-zero only for the nearest-neighbor pairs, the Hubbard Hamiltonian is transformed to the following effective Hamiltonian for large 
[67]:
(3)
where 
and 
and 
indicate the nearest-neighbor sites in the 
and 
directions, respectively. The second term contains the three-site terms when
. If we neglect the three-site terms, this effective Hamiltonian is reduced to the t-J model:
where
.
The Hubbard model has a long history in describing the magnetism of materials since the early works by Hubbard [26], Gutzwiller [68] and Kanamori [69]. Onedimensional Hubbard model has been well understood by means of the Bethe ansatz [70-72] and conformal field theory [73-75]. The solutions established a novel concept of the Tomonaga-Luttinger liquid [76] which is described by the scalar bosons corresponding to charge and spin sectors, respectively. The correlated electrons in twoand three-dimensional space are still far from a complete understanding in spite of the success for the one-dimensional Hubbard model. A possibility of superconductivity is a hot topic as well as the magnetism and metal-insulator transition for the twoand three-dimensional Hubbard model.
The three-band Hubbard model that contains 
and 
orbitals has also been investigated intensively with respect to high temperature superconductors [24,64, 77-88]. This model is also called the d-p model. The 2D three-band Hubbard model is the more realistic and relevant model for two-dimensional CuO2 planes which are contained usually in the crystal structures of high-
superconductors. The network of CuO2 layer is shown in Figure 3. The parameters of the three-band Hubbard model are given by the Coulomb repulsion
, energy levels of 
electrons 
and 
electron
, and transfer between 
orbitals given by
. Typical parameter values for the three-band (
-
) Hubbard model are shown in Table 1. The Hamiltonian of the
Figure 3. The lattice of the three-band Hubbard model on the CuO2 plane. Small circles denote Cu sites and large ones denote O sites.
Table 1. Typical parameter values for the three-band Hubbard model. Energies are measured in eV.
three-band Hubbard model is written as [24,47,48,80] (see Equation (4)) where 
and 
represent unit vectors along x and y directions, respectively. 
and 
denote the operators for the 
electrons at site
. Similarly 
and 
are defined. 
denotes the strength of Coulomb interaction between 
electrons. For simplicity we neglect the Coulomb interaction among 
electrons in this paper. Other notations are standard and energies are measured in 
units. The number of cells is denoted as 
for the three-band Hubbard model. In the non-interacting case 
the Hamiltonian in the k-space is written as (see Equation (5)).
where 
, 
and 
are operators for 
-, 
- and 
-electron of the momentum 
and spin
, respectively.
In the limit
, 
, and
, the 
-
model is mapped to the t-J model with
(6)
where
. 
gives the antiferromagnetic coupling between the neighboring 
and 
electrons. In real materials 
is not so large. Thus it seems that the mapping to the t-J model is not necessarily justified.
3. Variational Monte Carlo Studies
In this Section we present studies on the two-dimensional Hubbard model by using the variational Monte Carlo method.
3.1. Variational Monte Carlo Method
Let us start by describing the method based on the 2D Hubbard model. The Hamiltonian is given by
(7)
where 
denotes summation over all the nearestneighbor bonds and 
means summation over the next nearest-neighbor pairs. 
is our energy unit. The dispersion is given by
(8)
Our trial wave function is the Gutzwiller-projected wave functions defined as
(9)
(10)
|  (4)  (5) |
where
(11)
(12)
is the Gutzwiller projection operator given by
(13)
is a variational parameter in the range from 0 to unity and 
labels a site in the real space. 
is a projection operator which extracts only the sites with a fixed total electron number
. Coefficients 
and 
in 
appear in the ratio defined by
(14)
where 
and 
is a k-dependent gap function. 
is a variational parameter working like the chemical potential. 
is the Fourier transform of
. The wave functions 
and 
are expressed by the Slater determinants for which the expectations values are evaluated using a Monte Carlo procedure [43,44,92]. 
is written as
(15)
where
(16)
Then 
is written using the Slater determinants as
(17)
where 
is the Slater determinant defined by
(18)
In the process of Monte Carlo procedure the values of cofactors of the matrix in Equation (18) are stored and corrected at each time when the electron distribution is modified. We optimized the ground state energy
(19)
with respect to
, 
and 
for 
for
. For 
the variational parameter is only
. We can employ the correlated measurements method [93] in the process of searching optimal parameter values minimizing
.
A Monte Carlo algorithm developed in the auxiliary field quantum Monte Carlo calculations can also be employed in evaluating the expectation values for the wave functions shown above [94-96]. Note that the Gutzwiller projection operator is written as
(20)
where
. Then using the discrete HubbardStratonovich transformation, the Gutzwiller operator is the bilinear form:
(21)
where
. The Hubbard-Stratonovich auxiliary field 
takes the values of
. The norm 
is written as
(22)
where 
is a diagonal 
matrix corresponding to the potential
(23)
is written as
(24)
where 
denotes a diagonal matrix with elements given by the arguments
. The elements of 
are given by linear combinations of plane waves. For example,
(25)
Then we can apply the standard Monte Carlo sampling method to evaluate the expectation values [94,95]. This method is used to consider an off-diagonal Jastrow correlation factor of exp(–S)-type. The results for the improved wave functions are discussed in Section 3.10.
3.2. Superconducting Condensation Energy
We study the cases of the 
-, extended 
- (
-) and 
-wave gap functions in the following:
(26)
(27)
(28)
In Figure 4 calculated energies per site with 
on the 
lattice are shown for the case of 
and 
[46]. 
is plotted as a function of 
for three types of gap functions shown above. We impose the periodic and the antiperiodic boundary conditions for 
- and 
-direction, respectively. This set of boundary conditions is chosen so that 
does not vanish for any k-points occupied by electrons. 
was obtained as the average of the results of several Monte Carlo calculations each with 
steps. 
has minimum at a finite value of 
in the case of the d-wave gap function.
The energy gain 
in the superconducting state is called the SC condensation energy in this paper. 
is plotted as a function of 
in Figure 5 in order to examine the size dependence of the SC energy gain [97]. Lattice sizes treated are from 
to
. The electron density 
is in the range of
. Other parameters are 
and 
in 
units. Bulk limit 
of SC condensation energy 
was obtained by plotting as a function of
. The linear fitting line indicates very
Figure 5. Energy gain per site in the SC state with reference to the normal state for the 2D Hubbard model is plotted as a function of
. L is the length of the edge of the square lattice. YBCO attached to the vertical axis indicates the experimental value of the SC condensation energy for YBa2Cu3O4 [97].
clearly that the bulk limit remains finite when 
and
. When
, 
and
, the bulk-limit 
is
, where t = 0.51 eV is used [98]. Thus the superconductivity is a real bulk property, not a spurious size effect. The value is remarkably close to experimental values 0.17 ~ 0.26 meV/site estimated from specific heat data [99,100] and 0.26 meV/site from the critical magnetic field 
[101] for optimally doped YBa2Cu3O4 (YBCO). This good agreement strongly indicates that the 2D Hubbard model includes essential ingredients for the superconductivity in the cuprates.
We just point out that the t-J model gives 
at 
for J = 4 t2/U = 0.5 and 
[102]. This value is 50 times larger than the experimental values indicating a serious quantitative problem with this model. This means that the t-J model made from the leading two terms in the expansion in terms of 
of the canonical transformation of the Hubbard model should be treated with the higher-order terms in order to give a realistic SC condensation energy.
Here we show the SC condensation energy as a function of 
in Figure 6. The condensation energy 
is increased as 
is increased as far as
. In the strong coupling region
, we obtain the large condensation energy.
3.3. Fermi Surface and Condensation Energy
Now let us consider the relationship between the Fermi surface structure and the strength of superconductivity. The experimental SC condensation energy for (La,Sr)2CuO4 (LSCO) is estimated at 0.029 meV/(Cu site) or 0.00008 in units of t which is much smaller than that
Figure 6. Energy gain per site in the SC state with reference to the normal state for the 2D Hubbard model as a function of the Coulomb repulsion U. The system is 10 × 10 with the electron number 
and
.
for YBCO. [103] The band parameter values of LSCO were estimated as 
and 
[104]. This set corresponds roughly to
. The latter value is much larger than the above-mentioned experimental value for LSCO. However, the stripe-type SDW state coexists with superconductivity [105,106] and the SC part of the whole 
is much reduced. Therefore, such a coexistence allows us to qualitatively understand the SC 
in LSCO.
On the other hand, Tl2201 
and Hg1201 
band calculations by Singh and Pickett [107] give very much deformed Fermi surfaces that can be fitted by large 
such as
. For Tl2201, an Angular Magnetoresistance Oscillations (AMRO) work [108] gives information of the Fermi surface, which allows to get 
and
. There is also an Angle-Resolved Photoemission Study (ARPES) [109], which provides similar values. In the case of Hg1201, there is an ARPES work [110], form which we obtain by fitting 
and
. For such a deformed Fermi surface, 
in the bulk limit is reduced considerably. [111,112] Therefore, the SC 
calculated by VMC indicates that the Fermi surface of LSCO-type is more favorable for high
. The lower 
in LSCO may be attributed to the coexistence with antiferromagnetism of stripe type.
3.4. Ladder Hubbard Model
The SC condensation energy in the bulk limit for the ladder Hubbard model has also been evaluated using the variational Monte Carlo method [54]. The Hamiltonian is given by the 1D two-chain Hubbard model: [51,52,55, 56,85,113-116]
(29)
where 
is the creation (annihilation) operator of an electron with spin 
at the 
th site along the 
th chain
. 
is the intrachain nearest-neighbor transfer and 
is the interchain nearest-neighbor transfer energy. The energy is measured in 
units. The energy minimum was given when the components of the gap function 
take finite values plotted in Figure 7 for the lattice of 
sites with 34 electrons imposing the periodic boundary condition [54]. Each component of 
was optimized for 
and
. There are two characteristic features; one is that the components of the bonding and antibonding bands have opposite signs each other and second is that the absolute values of 
of the antibonding band 
is much larger than that of the bonding band
. In order to reduce the computation cpu time,
of each band was forced to take a fixed value specific to each band, i.e. 
for the bonding band and 
for the antibonding band. This drastically reduces the number of the variational parameters but still allows us to get a substantial value of the condensation energy. 
and 
take opposite sign, which is similar to that of the 
gap function.
The energy gain 
remains finite in the bulk limit when
. The SC condensation energy per site in the bulk limit is plotted as a function of 
in Figure 8 [54]. The SC region derived from the SC condensation energy in the bulk limit is consistent with the results obtained from the density-matrix renormalization group [55,56] and the exact-diagonalization method [51,52,115]. The maximum value of 
is 0.0008 which is of the same order of magnitude as the maximum condensation energy obtained for the 2D Hubbard model [46].
Figure 8. 
dependence of the SC condensation energy 
for the two-chain Hubbard model in the bulk limit [54].
3.5. Condensation Energy in the d-p Model
The SC energy gain for the d-p model, namely, threeband Hubbard model in Equation (4) has also been evaluated using the variational Monte Carlo method. For the three-band model the wave functions are written as
(30)
(31)
where 
is the linear combination of
, 
and 
constructed to express the operator for the lowest band (in the hole picture) or the highest band (in the electron picture) of the non-interacting Hamiltonian. The numerical calculations have been done in the hole picture. The Gutzwiller parameter
, effective level difference
, chemical potential 
and superconducting order parameter 
are the variational parameters.
The similar results to the single-band Hubbard model were obtained as shown in Figure 9 for
, 
and 
in 
units where the calculations were performed in the hole picture [24]. The SC condensation energy for the three-band model is 
meV per site in the optimally doped region. We set 
as estimated in Table 1. There is a tendency that 
increases as 
increases which is plotted in Figure 10. This tendency is not, however, in accordance with NQR (nuclear quadrupole resonance) study on cuprates. [117] We think that the NQR experiments indicate an importance of the Coulomb interaction on oxygen sites. This will be discussed in Section 3.11.
3.6. Antiferromagnetic State
When the density of doped holes is zero or small, the 2D single-band or three-band Hubbard model takes an antiferromganetic state as its ground state. The magnetic
order is destroyed and superconductivity appears with the increase of doped hole density. The transition between the d-wave SC and the uniform SDW states has been investigated by computing the energy of the SDW state by using the variational Monte Carlo method. The trial SDW wave function is written as
(32)
(33)
(34)
(35)
(36)
Summation over 
and 
in Equation (33) is performed over the filled k-points, as in the calculation of the normal state energy. 
is the SDW wave vector given by 
and 
is the SDW potential amplitude.
As shown in Figure 11, the energy gain per site in the SDW state rises very sharply from 
for 
[46]. At 
it is slightly larger than that in the SC state, and at 
there is no more stable SDW state. Thus the boundary between the SDW and the SC states is given at
. The results of the bulk limit calculations indicate that the energy gain in the SC state at 
takes the extremely small value and the value at 
vanishes as
. Hence the pure d-wave SC state possibly exists near the boundary at
, but the region of pure SC state is very restricted.
Let us turn to the three-band model. We show the antiferromagnetic-paramagnetic boundary for 
and 
in the plane of 
and the hole density in Figure 12 where AF denotes the antiferromagnetic region [47]. The value 
is taken from the estimations by cluster calculations [89-91]. The phase boundary in the region of small 
is drawn from an extrapolation. For the intermediate values of
,
Figure 12. Antiferromagnetic region in the plane of U and the hole density for 
and
.
the antiferromagnetic long-range ordering exists up to about 20 percent doping. Thus the similar features are observed compared to the single-band Hubbard model.
Since the three-band Hubbard model contains several parameters, we must examine the parameter dependence of the energy of SDW state. The energy gain 
in the SDW state is shown in Figure 13 as a function of doping ratio for several values of
. 
increases as 
increases as expected. In Figure 14 
- and 
-dependencies of 
are presented. The SDW phase extends up to 30 percent doping when 
is large. It follows from the calculations that the SDW region will be reduced if 
and 
decrease.
From the calculations for the SDW wave functions, we should set 
and 
small so that the SDW phase does not occupy a huge region near half-filling. In Figures 15 and 16, we show energy gains for both the SDW and SC states for
, 
, 
and
, where the right hand side and left hand side indicate the hole-doped and electron-doped case, respectively. Solid symbols indicate the results for 
and open symbols for
. For this set of parameters the SDW region extends up to 20 percent doping and the pure d-wave phase exists outside of the SDW phase. The d-wave phase may be possibly identified with superconducting phase in the overdoped region in the high-
superconductors.
3.7. Stripes and Its Coexistence with Superconductivity
Incommensurate magnetic and charge peaks have been observed from the elastic neutron-scattering experiments in the underdoped region of the Nd-doped La2–x–yNdySrxCuO4 [118] (Figure 17). Recent neutron experiments have also revealed the incommensurate spin structures [119-123]. Rapid decrease of the Hall resistivity in this region suggests that the electric con-
duction is approximately one dimensional [124]. The angle-resolved photo-emission spectroscopy measurements also suggested a formation of two sets of onedimensional Fermi surface [125]. Then it has been proposed that these results might be understood in the framework of the stripe state where holes are doped in the domain wall between the undoped spin-density-wave
domains. This state is a kind of incommensurate SDW state. It was also shown that the incommensurability is proportional to the hole density in the low-doping region in which the hole number per stripe is half of the site number along one stripe [118,120]. A static magnetically ordered phase has been observed by 
SR over a wide range of SC phase for 
in La2–xSrxCuO4 (LSCO) [126]. Thus the possibility of superconductivity that occurs in the stripe state is a subject of great interest [127-130]. The incommensurate magnetic scattering spots around 
were observed in the SC phase in the range of 
in the elastic and inelastic neutron-scattering experiments with LSCO [127,128,130]. The hole dependence of the incommensurability and the configuration of the spots around the Bragg spot in the SC phase indicated the vertical stripe. The neutronscattering experiments have also revealed that a diagonal spin modulation occurs across the insulating spin-glass phase in La2–xSrxCuO4 for
, where a one-dimensional modulation is rotated by 45 degrees from the stripe in the SC phase. The incommensurability 
versus hole density is shown in Figure 18 schematically [129,130]. The diagonal stripe changes into the vertical stripe across the boundary between the insulating and SC phase.
Let us investigate the doped system from the point of modulated spin structures [131-141]. The stripe SDW state has been studied theoretically by using the meanfiled theory [132-136]. They found that the stripe state appears when an incommensurate nesting becomes
Figure 17. Charge and spin density as a function of the distance for a striped state [50].
Figure 18. Schematic illustration of the incommensurability versus hole density.
favorable in the hole-doped 2D Hubbard model. When the electron correlation correlation is strong or intermediate, it was shown that the stripe state is more stable than the commensurate spin-density-wave state with the wave vector 
in the ground state of the 2D Hubbard model by using the variational Monte Carlo method [131]. It has also been confirmed by the same means that the stripe states are stabilized in the d-p model [48]. The purpose of this section is to examine whether the superconductivity can coexist with static stripes in the 2D Hubbard model in a wider doping region and investigate the doping dependence of the coexisting state.
We consider the 2D Hubbard model on a square lattice. We calculate the variational energy in the coexistent state that is defined by
(37)
where 
is a mean-field wave function. The effective mean-field Hamiltonian for the coexisting state is [105] represented by
(38)
where the diagonal terms describe the incommensurate spin-density wave state:
(39)
where 
is the chemical potential. The vertical stripe state is represented by the charge density 
and the spin density 
that are spatially modulated as
(40)
(41)
where 
denotes the position of vertical stripes. The amplitude 
is fixed by
. The off-diagonal terms in Equation (38) are defined in terms of the d-wave SC gap as
(42)
where
, 
(unit vectors). We consider two types of the spatially inhomogeneous superconductivity: anti-phase and in-phase defined as
(43)
(44)
and
(45)
(46)
respectively. Here, 
and 
is a incommensurability given by the stripe’s periodicity in the y direction with regard to the spin. The anti-phase (inphase) means that the sign if the superconducting gap is (is not) changed between nearest domain walls.
The wave function 
is made from the solution of the Bogoliubov-de Gennes equation represented by
(47)
(48)
The Bogoliubov quasiparticle operators are written in the form
(49)
(50)
Then the coexistence wave function is written as [105,142]
(51)
for constants 
and
. The calculations are performed for the wave function
. The variational parameters are
, 
, 
, 
and
. The system parameters were chosen as 
and 
suitable for cuprate superconductors. It has been shown that the “anti-phase” configuration is more stable than the “in-phase” one [105].
Here, the system parameters are 
and
. We use the periodic boundary condition in the x-direction and anti-periodic one in the y-direction. In Figure 19, we show the total energy of the coexistent state, 
, as a function of the SC gap 
for the cases of anti-phase and in-phase. The SC condensation energy 
is estimated as 0.0008 t per site at the hole density 0.125 on the 
lattice with periodic boundary condition in x-direction and antiperiodic one in y-direction. 
in the coexistence state is defined as the decrease of energy due to finite
. If we use
, this is evaluated as
. The SC condensation energy per site is shown as a function of hole density in Figure 20. One finds that 
in the stripe state decreases as the hole density decreases. This tendency is reasonable since the SC order is weakened in the domain of the incommensurate SDW because of the vanishingly small carrier concentration contributing the superconductivity in this domain. This behavior is consistent with the SC condensation energy estimated from the specific heat measurements [143].
There is a large renormalization of the Fermi surface due to the correlation effect in the striped state [144]. We considered the next-nearest transfer 
in the trial function as a variational parameter
. In Figure 21,
Figure 20. Superconducting condensation energy per site in the coexistence state as a function of the hole density
, 0.10 and 0.125. The model is the single-band Hubbard Hamiltonian with
. The stripe interval is preserved constant. The inset shows the hole dependence of the incommensurability in the coexistent state [105].
optimized values of 
for the striped state are shown as a function of
. The 
increases as 
increases. We also mention that the optimized 
almost vanishes. The renormalized Fermi surface of
, 
and 
are plotted in Figure 22. The system is a 
lattice with 
and the electron density 0.875. As 
is incresed, the Fermi surface is more deformed. We show the the gradient of the momentum distribution function, 
, calculated in the optimized stripe state in Figure 23. The brighter areas coincide with the renormalized Fermi surface with 
and 
for
.
The calculations for the three-band Hubbard model has also been done taking into account the coexistence of stripes and SC [15,106]. The energy of antiferromagnetic state would be lowered further if we consider the incommensurate spin correlation in the wave function. The phase diagram in Figure 24 presents the region of stable AF phase in the plane of 
and
. For large
, we have the region of the AF state with an eight-lattice periodicity in accordance with the results of neutron-scattering measurements [118,123]. The energy at 
is shown in Figure 25 where the 4-lattice stripe state has higher energy than that for 8-lattice stripe for all the values of
.
The Bogoliubov-de Gennes equation is extended to the case of three orbitals
, 
and
. 
and 
are now 
matrices. The energy of the state with double order parameters 
and 
is shown in Figure 26 on the 
lattice at the doping rate 1/8. The SC condensation energy per site is evaluated as 
for
, 
and
. If we use 
[89-91], we obtain 
which is slightly smaller than and close to the value obtained for the single-band Hubbard model. We show the size dependence of the SC condensation
Figure 23. Contour plot of 
measured for the projected stripe state on 24 × 24 lattice with
. The electron density is 0.875.
Figure 24. Phase diagram of stable antiferromagnetic state in the plane of 
and 
obtained for 16 × 4 lattice.
Figure 25. Energy as a function of 
for 16 × 16 square lattice at
. Circles, triangles and squares denote the energy for 4-lattice stripes, 8-lattice stripes, and commensurate SDW, respectively, where n-lattice stripe is the incommensurate state with one stripe per n ladders. The boundary conditions are antiperiodic in x-direction and periodic in y-direction [47,106].
energy for
, 0.125 0.08333 and 0.0625 in Figure 27. We set the parameters as 
and 
in 
units, which is reasonable from the viewpoint of the density of states and is remarkably in accordance with cluster estimations [89-91], and also in the region of eight-lattice periodicity at
. We have carried out the Monte Carlo calculations up to 
sites (768 atoms in total). In the overdoped region in the range of
, we have the uniform d-wave pairing state as the ground state. The periodicity of spatial variation increases as the doping rate 
decreases proportional to
. In the figure we have the 12-lattice periodicity at 
and 16-lattice periodicity at
. For
, 0.125 and 0.08333, the results strongly suggest a finite condensation energy in the bulk limit. The SC condensation energy obtained on the basis of specific heat measurements agrees well with the variational Monte Carlo computations [99]. In general, the Monte Carlo statistical errors are much larger than those for the single-band Hubbard model. The large number of Monte Carlo steps (more than
) is required to get convergent expectation values for each set of parameters.
In Figure 28 the order parameters 
and 
Figure 28. Phase diagram of the d-p model based on the Gutzwiller wave function [106].
were evaluated using the formula 
where 
is the density of states. Here we have set 
since 
is estimated as 
to 3 
for optimally doped YBCO using
[100]. The phase diagram is consistent with the recently reported phase diagram for layered cuprates [145].
3.8. Diagonal Stripe States in the Light-Doping Region
Here we examine whether the relationship 
holds in the lower doping region or not, and whether the diagonal stripe state is obtained in the further lower doping region [50]. The elastic neutron scattering experiments of LSCO in the light-doping region, 
, revealed that the position of incommensurate magnetic peaks changed from 
to 
as 
becomes less than 0.06 [129,130]. This means that the stripe direction rotates by 45 degrees to become diagonal at this transition. In the diagonal stripe states, the magnetic peaks are observed to keep a relation 
that holds in the vertical stripe state in the low doping region.
In Figure 29, we show the incommensurability of the most stable stripe state as a function of
. Open squares and triangles are values for diagonal and vertical incommensurate SDW’s obtained in the elastic neutron scattering experiments on LSCO, respectively. Solid squares and triangles show our results for the diagonal and vertical stripes, respectively. These results are in a good agreement with experimental data. We also found that the phase boundary 
between the diagonal and vertical stripe states lies at or above 0.0625 in the
case of 
and
. The following factors may give rise to slight changes of the phase boundary
: the diagonal stripe state may be stabilized in the low-temperature-orthorhombic (LTO) phase in LSCO. The diagonal stripe state is probably stabilized further by forming a line along larger next-nearest hopping direction due to the anisotropic 
generated by the Cu-O buckling in the LTO phase.
3.9. Checkerboard States
A checkerboard-like density modulation with a roughly 
period (
is a lattice constant) has also been observed by scanning tunneling microscopy (STM) experiments in Bi2Sr2CaCu2O8+d, [146] Bi2Sr2–xLaxCuO6+d, [147] and Ca2–xNaxCuO2Cl2 (Na-CCOC) [148]. It is important to clarify whether these inhomogeneous states can be understood within the framework of strongly correlated electrons.
Possible several electronic checkerboard states have been proposed theoretically [134,149,150]. The charge density 
and spin density 
are spatially modulated as
(52)
(53)
where 
and 
are variational parameters. The striped incommensurate spin-density wave state (ISDW) is defined by a single Q vector. On the other hand, the checkerboard ISDW state is described by the double-Q set; for example, vertical wave vectors 
and 
describe a spin vertical checkerboard state, where two diagonal domain walls are orthogonal. Diagonal wave vectors 
and 
lead to a spin diagonal checkerboard state with a 
-period. The hole density forms the charge vertical checkerboard pattern with vertical wave vectors 
and 
in which the hole density is maximal on the crossing point of magnetic domain walls in the spin diagonal checkerboard state. If 
is assumed, the charge modulation pattern is consistent with the 
charge structure observed in STM experiments.
We found that the coexistent state of bond-centered four-period diagonal and vertical spin-checkerboard structure characterized by a multi-Q set is stabilized and composed of 
period checkerboard spin modulation. [151] In Figure 30(a), we show the condensation energies of some heterogeneous states, 
, with fixing the transfer energies 
and 
suitable for Bi-2212. The system is a 
lattice with the electron-filling
. The energy of the normal state 
is calculated by adopting
. In our calculations, the condensation energies of both bondcentered stripe and checkerboard states are always larger than those of site-centered stripe and checkerboard states. The vertical stripe state is not suitable in this parameter set since this state is only stabilized with the LSCO-type band. The four-period spin-diagonal checkerboard (DC) state is significantly more stable than the eight-period spin-DC state. We found that the coexistent state of the bond-centered four-period spin-DC and four-period spin-vertical checkerboard (VC) with 
is most stable, as shown in Figure 30(a). The measured expectation value of the spin density is shown in Figure 30(b).
3.10. Improved Gutzwiller Function
We have presented the results based on the Gutzwiller functions for the normal state, SDW state and BCS state. We must consider a method to go beyond the Gutzwiller function-based Monte Carlo method. One method to achieve this purpose is to multiply the Jastrow correlation operators to take into account the intersite correlations. The simplest possible candidate is an introduction of the diagonal intersite correlation factor [152]:
(54)
for the variational parameter
. We have investigated the 2D Hubbard model by using the JastrowGutzwiller function [111]. The ground-state energy is lowered considerably by considering the intersite correlations such as nearest neighbor and next nearest neighbor spin and charge correlations.
Here we consider the method to improve the wave functions by an off-diagonal Jastrow correlation operators [94,95,153]. The off-diagonal correlation factors are more effective to lower the ground state energy in 2D systems. Let us consider the wave functions 
defined in the following way: [95]
(55)
(56)
(57)
(58)
and so on, where 
denotes the kinetic part of the Hamiltonian
(59)
and 
denotes the on-site Coulomb interaction.
, 
, 
, 
, 
and 
are variational parameters. It is considered that 
approaches the true ground state wave function as 
grows larger since the ground state wave function is written as
(60)
for large 
and small 
. If we can extrapolate the expectation values from the data for
, we can evaluate correct expectation values.
The computations are performed by using the discrete Hubbard-Stratonovich transformation as described in Section 3.1. In the evaluation of the expectation values we generate the Monte Carlo samples by the importance sampling [95] with the weight function 
where
(61)
Since the Monte Carlo samplings are generated with the weight
, the expectation values are calculated with the sign of 
in the summation over the generated samples. In our calculations the negative sign problem has become less serious due to the variational treatment, although we encounter the inevitable negative sign problem in the standard projector Monte Carlo approaches [154].
In Figure 31 the energy is shown as a function of 
where the SDW and normal states are chosen as the initial state
. The extrapolated values for different initial states coincide with each other within Monte Carlo statistical errors. The energy expectation values as a function of 
for 
square lattice are presented in Figure 32 for
, 
,
. The extrapolated curve is shown by the solid curve and the results obtained by the quantum Monte Carlo simulation (QMC) [28] are also shown as a reference. A good agreement between two calculations support the method although the QMC gives slightly higher energy for
.
One can formulate an approach to consider the BCS function with correlation operators. [96] For this end the electron-hole transformation is introduced for the down spin and the up-spin electrons are unaltered [155]. We show the energy versus 
in Figure 33 for 
and
. From an extrapolation to the limit
, both formulations predict the same limiting value for the energy. The energy is lowered considerably due to the
correlation operators compared to that for the Gutzwiller function. The energy in the d-wave SC state is always lower than that in the normal state for each
. The energy gain in the SC state remains the same order after the multiplication of correlation factors.
3.11. 
and 
Relationships between 
and structural features in cuprate high-temperature superconductors are very interesting. Torrance and Metzger found the first such relationship between 
and the Madelung potential difference 
[156]. Here 
is the potential difference between Cu and O sites in the CuO
plane. 
was found to increase with decreasing
. There is an interesting tendency of increasing 
with increasing relative ratio of hole density at oxygen site against that at copper site [117].
Here we show the results obtained by using the perturbation theory [62-66], There have been many similar works by making some kind of approximation such as random phase approximation (RPA) [157-159], fluctuation-exchange approximation (FLEX) [160-163], effective spin-fluctuation method [4,164,165], and perturbation theory in terms of U [166-168]. An application was made for Sr2RuO4 where we need to consider th emulti-band structure α and β orbitals [169], and also to the three-dimensional d-p model [170]. In our formulation the gap function is written as
(62)
The exponent 
indicates the strength of superconductivity. The results are in Figures 34 and 35 [171]. As shown in the figure, for positive
, with increase of 
the exponent 
increases monotonously. This means the increase of superconducting gap and so of
, and is consistent with the wide-range tendency of the variational Monte Carlo calculation [24,172]. This tendency can be understood in terms of
Figure 34. The exponent x (superconductivity strength) as a function of
, where the level difference 
is positive.
Figure 35. The exponent x as a function of
, where the level difference 
is negative.
(63)
where 
is the weight ofd electrons. This clearly indicates that increase of 
leads to the increase of 
and subsequently of
. In the case of
, we take account of finite Coulomb repulsion 
on oxygen sites. The effective interaction coming from 
is similarly given by the susceptibility with the weight of 
electrons. The results of 
with 
indicates that all four types of even parity (
, 
, 
and
) SC strength values increase, so that 
is raised, as the absolute value 
increases in this region. This result shows that 
also plays an important role as well.
Let us give a discussion on this result. Increase of 
in the region of 
means decrease of 
since
, where 
is the second electron affinity of oxygen atom and 
is the third electron ionicity of copper atom and 
is the charge of electron. Therefore, this relation is consistent with the systematics reported in [156]. With increase of the distance of the apex oxygen away from the CuO2 plane, cuprate superconductors are known to increase 
[173]. The accompanying raise of 
should tend to increase
.
The Coulomb interaction between p electrons on oxygen atom will raise the level of p electrons effectively. This leads to the lowering of p hole level 
or the raise of 
relatively. This indicates that 
will be increased by the Coulomb interaction between p electrons.
4. Quantum Monte Carlo Studies
4.1. Quantum Monte Carlo Method
The Quantum Monte Carlo (QMC) method is a numerical method that is employed to simulate the behavior of correlated electron systems. We outline the QMC method in this section. The Hamiltonian is the Hubbard model that contains the on-site Coulomb repulsion and is written as
(64)
where 
is the creation (annihilation) operator of an electron with spin 
at the 
-th site and
. 
is the transfer energy between the sites 
and
. 
for the nearest-neighbor bonds. For all other cases
. 
is the on-site Coulomb energy. The number of sites is 
and the linear dimension of the system is denoted as
. The energy unit is given by 
and the number of electrons is denoted as
.
In a Quantum Monte Carlo simulation, the ground state wave function is
(65)
where 
is the initial one-particle state represented by a Slater determinant. For large
, 
will project out the ground state from
. We write the Hamiltonian as 
where K and V are the kinetic and interaction terms of the Hamiltonian in Equation (64), respectively. The wave function in Equation (65) is written as
(66)
for
. Using the Hubbard-Stratonovich transformation [27,94], we have
(67)
for 
or
. The wave function is expressed as a summation of the one-particle Slater determinants over all the configurations of the auxiliary fields
. The exponential operator is expressed as
(68)
where we have defined
(69)
for
(70)
(71)
The ground-state wave function is
(72)
where 
is a Slater determinant corresponding to a configuration 
of the auxiliary fields:
(73)
The coefficients 
are constant real numbers:
. The initial state 
is a one-particle state. If electrons occupy the wave numbers 
for each spin
, 
is given by the product 
where 
is the matrix represented as [31]
(74)
is the number of electrons for spin
. In actual calculations we can use a real representation where the matrix elements are 
or
. In the real-space representation, the matrix of 
is a diagonal matrix given as
(75)
The matrix elements of 
are
(76)
is an 
matrix given by the product of the matrices
, 
and
. The inner product is thereby calculated as a determinant [38],
(77)
The expectation value of the quantity 
is evaluated as
(78)
If 
is a bilinear operator 
for spin
, we have
(79)
The expectation value with respect to the Slater determinants 
is evaluated using the singleparticle Green’s function [31,38],
(80)
In the above expression,
(81)
can be regarded as the weighting factor to obtain the Monte Carlo samples. Since this quantity is not necessarily positive definite, the weighting factor should be
; the resulting relationship is,
(82)
where 
and
(83)
This relation can be evaluated using a Monte Carlo procedure if an appropriate algorithm, such as the Metropolis or heat bath method, is employed [94]. The summation can be evaluated using appropriately defined Monte Carlo samples,
(84)
where 
is the number of samples. The sign problem is an issue if the summation of 
vanishes within statistical errors. In this case it is indeed impossible to obtain definite expectation values.
4.2. Quantum Monte Carlo Diagonalization
4.2.1. Basic Method and Optimization
Quantum Monte Carlo diagonalization (QMD) is a method for the evaluation of 
without the negative sign problem [41]. A bosonic version of this method was developed before in Ref.[174]. The configuration space of the probability 
in Equation (84) is generally very strongly peaked. The sign problem lies in the distribution of 
in the configuration space. It is important to note that the distribution of the basis functions 
is uniform since 
are constant numbers:
. In the subspace
, selected from all configurations of auxiliary fields, the right-hand side of Equation (78) can be determined. However, the large number of basis states required to obtain accurate expectation values is beyond the current storage capacity of computers. Thus we use the variational principle to obtain the expectation values.
From the variational principle,
(85)
where 
are variational parameters. In order to minimize the energy
(86)
the equation 
is solved for,
(87)
If we set
(88)
(89)
the eigen-equation is
(90)
for
. Since 
are not necessarily orthogonal, 
is not a diagonal matrix. We diagonalize the Hamiltonian
, and then calculate the expectation values of correlation functions with the ground state eigenvector.
In Quantum Monte Carlo simulations an extrapolation is performed to obtain the expectation values for the ground-state wave function. If 
is large enough, the wave function in Equation (72) will approach the exact ground-state wave function, 
, as the number of basis functions, 
, is increased. If the number of basis functions is large enough, the wave function will approach, 
, as 
is increased. In either case the method employed for the reliable extrapolation of the wave function is a key issue in calculating the expectation values. The variance method was recently proposed in variational and Quantum Monte Carlo simulations, where the extrapolation is performed as a function of the energy variance. We can expect linearity in some cases [175]:
(91)
where 
denotes the variance defined as
(92)
and 
is the expected exact value of the quantity
.
The simplest procedure for optimizing the groundstate wave function is to increase the number of basis states 
by random sampling. First, we set 
and
, for example, 
, and
. We denote the number of basis functions as
. We start with 
and then increase up to 10,000. This procedure can be outlined as follows:
1) Generate the auxiliary fields 
in 
randomly for 
for 
, and generate 
basis wave function
.
2) Evaluate the matrices 
and
, and diagonalize the matrix 
to obtain
. Then calculate the expectation values and the energy variance.
3) Repeat the procedure from 1) after increasing the number of basis functions.
For small systems this random method produces reliable energy results. The diagonalization plays an importance producing fast convergence. In order to lower the ground-state energy efficiently, we can employ a genetic algorithm [176] to generate the basis set from the initial basis set. One idea is to replace some parts of 
in 
that has the large weight 
to generate a new basis function
. The new basis function 
obtained in this way is expected to also have a large weight and contribute to
. The details of the method are shown in Ref.[41].
4.2.2. Ground State Energy and Correlation Functions
The energy as a function of the variance is presented in Figures 36-38 for
, 
and
, respectively. To obtain these results the genetic algorithm was employed to produce the basis functions except the open symbols in Figure 4. The 
where 
in Figure 2 is the energy for the closed shell case up to 2000 basis states. The other two figures are for open shell cases, where evaluations were performed up to 3000 states. We show the results for the
, 
and 
systems in Table 2.
The Figure 39 is the momentum distribution function
,
(93)
for 
sites where the results for the Gutzwiller VMC and the QMD are indicated. The Gutzwiller function gives the results that 
increases as 
approaches 
from above the Fermi surface. This is clearly unphysical. This flaw of the Gutzwiller function near the Fermi surface is not observed for the QMD result.
4.2.3. Spin Gap in the Hubbard Ladder
Here we show the results for one-dimensional models. The ground state of the 1D Hubbard model is no longer Fermi liquid for
. The ground state is insulating at half-filling and metallic for less than half-filling. The Figure 40 is the spin and charge correlation functions, 
and
, as a function of the wave number, for
Table 2. Ground state energy per site from the Hubbard model. The boundary conditions are periodic in both directions. The current results are presented under the column labeled QMD. The constrained path Monte Carlo (CPMC) results are from Ref.[38]. The column VMC is the results obtained for the optimized variational wave function 
except for 6 × 2 for which 
is employed. The QMC results are from Ref.[35]. Exact results are obtained using diagonalization [177].
Figure 37. Energy as a function of the variance for 6 × 2 
and
. The square is the exact value obtained using exact diagonalization.
Figure 38. Energy as a function of the variance 
for 6 × 6. with the periodic boundary conditions. Solid circles and crosses are data obtained from the QMD method for two different initial configurations of the auxiliary fields. Gray open circles show results obtained from the 
-renormalization method with 300 basis wave functions.
Figure 39. Momentum distribution function for the 14 × 14 lattice. Parameters are 
and
. The boundary conditions are periodic in both directions. The results for the Gutzwiller function (open circle) are also provided.
the 1D Hubbard model where
. The 
singularity can be clearly identified where the dotted line is for
. The spin correlation is enhanced and the charge correlation function is suppressed slightly because of the Coulomb interaction.
The spin correlation function 
for the Hubbard ladder is presented in Figure 41, where 
and
. 
is defined as
(94)
where 
denotes the site 
. We use the convention that 
where 
and 
indicate the lower band and upper band, respectively. There are four singularities at
, 
, 
, and 
for the Hubbard ladder, where 
and 
are the Fermi wave numbers of the lower and upper band, respectively.
Figure 40. Spin (solid circle) and charge (open circle) correlation functions obtained from the QMD method for the one-dimensional Hubbard model with 80 sites. The number of electrons is 66. We set 
and use the periodic boundary condition.
It has been expected that the charge gap opens up as 
turns on at half-filling for the Hubbard ladder model. In Figure 42 the charge gap at half-filling is shown as a function of
. The charge gap is defined as
(95)
where 
is the ground state energy for the 
electrons. The charge gap in Figure 42 was estimated using the extrapolation to the infinite system from the data for the
, 
, and 
systems. The data suggest the exponentially small charge gap for small 
or the existence of the critical value 
in the range of
, below which the charge gap vanishes.
4.2.4. Magnetization in 2D Hubbard Model
The ground state of the 2D Hubbard model at half-filling is antiferromagnetic for 
because of the nesting due to the commensurate vector
. The Gutzwiller function predicts that the magnetization
(96)
increases rapidly as 
increases and approaches 
for large
. In Figure 43, the QMD results are presented for 
as a function of
. The previous results obtained using the QMC method are plotted as open circles. The gray circles are for the 
-function VMC method and squares are the Gutzwiller VMC data. Clearly, the magnetization is reduced considerably because of the fluctuations, and is smaller than the Gutzwiller VMC method by about 50 percent.
4.3. Pair Correlation Function
The pair correlation function 
is defined by
Figure 42. Charge gap as a function of 
for 
(circles). The DMRG results (squares) are provided for comparison [58].
(97)
where 
, 
, denote the annihilation operators of the singlet electron pairs for the nearestneighbor sites:
(98)
Here 
is a unit vector in the 
-direction. We consider the correlation function of d-wave pairing:
(99)
where
(100)
and 
denote sites on the lattice.
We show how the pair correlation function is evaluated in quantum Monte Carlo methods. We show the pair correlation functions 
and 
on the lattice 
in Figures 44 and 45. The boundary condition is open in the 4-site direction and is periodic in the other direction. An extrapolation is performed as a function of 
in the QMC method with Metropolis algorithm and as a function of the energy variance 
in the QMD method with diagonalization. We keep 
a small constant 
and and increase
, where 
is the division number. In the Metropolis QMC method, we calculated averages over 
Monte Carlo steps. The exact values were obtained by using the exact diagonalization method. Two methods give consistent results as shown in figures. All
the 
and 
are suppressed on 
as 
is increased. In general, the pair correlation functions are suppressed in small systems. In Figures 46 and 47, we show the inter-chain pair correlation function for the ladder model
. We use the open boundary condition. The number of electrons is
, and the strength of the Coulomb interaction is
. 
indicates the electron pair along the rung, and 
is the expectation value of the parallel movement of the pair along the ladder. The results obtained by two methods are in good agreement except 
(nearest-neighbor correlation).
We first consider the half-filled case with
; in this case the antiferromagnetic correlation is dominant over the superconductive pairing correlation and thus the pairing correlation function is suppressed as the Coulomb repulsion 
is increased. The Figure 48 exhibits the d-wave pairing correlation function 
on 
lattice as a function of the distance. The 
is suppressed due to the on-site Coulomb interaction, as expected. Its reduction is, however, not so considerably large compared to previous QMC studies [39] where the pairing correlation is almost annihilated for
. We then turn to the case of less than half-filling. We show the results on 
with electron number
. We show 
as a function of the distance in Figure 49
. In the scale of this figure, 
for 
is almost the same as that of the non-interacting case, and is enhanced slightly for large
. Our results indicate that the pairing correlation is not suppressed and is indeed enhanced by the Coulomb interaction
, and its enhancement is very small.
The Figure 50 shows 
on 
lattice. This also indicates that the pairing correlation function is
enhanced for
. There is a tendency that 
is easily suppressed as the system size becomes small. We estimated the enhancement ratio compared to the non-interacting case 
at 
for 
as shown in Figure 51. This ratio increases as the system size is increased. To compute the enhancement, we picked the sites, for example on 
lattice, 
, (4,0), (4,1), (3,3), (4,2), (4,3), (5,0), (5,1) with 
and evaluate the mean value. In our computations, the ratio increases almost linearly indicating a possibility of superconductivity. This indicates 
for
.
Because
, we obtain
for
. This indicates that the exponent of the power law is 2. When
, the
enhancement is small and is almost independent of
. In the low density case, the enhancement is also suppressed being equal to 1. In Figure 52, the enhancement ratio is shown as a function of the electron density 
for
. A dome structure emerges even in small systems. The square in Figure 52 indicates the result for the half-filled case with 
on 
lattice. This is the open shell case and causes a difficulty in computations as a result of the degeneracy due to partially occupied electrons. The inclusion of 
enhances 
compared to the case with 
on 
lattice. 
is, however, not enhanced over the non-interacting case at half-filling. This also holds for
lattice where the enhancement ratio
. This indicates the absence of superconductivity at half-filling.
4.4. Spin Susceptibility
We propose a method to compute the magnetic susceptibility at absolute zero 
[178]. We add the source term 
to the Hamiltonian as follows
(101)
where 
is a small real number of the order 
or
. We calculate 
in the ground state, which is, as shown by the linear response theory, the magnetic susceptibility
(102)
in the limit of small
, where 
is the retarded Green function and 
is the dynamical susceptibility,
(103)
Indeed, the above formula gives the correct spin susceptibility 
on the finite lattice for the noninteracting case, which is given by 
with the Fermi distribution function
. We calculate 
by using the quantum Monte Carlo method to obtain
.
We examine the results obtained for the susceptibilities. Figure 53 shows the spin susceptibility 
for 
on a 
lattice as a function of 
or the energy variance
. The number of electrons is 10. The expectation values agree well with exact values given by the exact diagonalization method.
We now compute the staggered susceptibility 
by adding the source term
to the Hamiltonianwhere
. Here we set the periodic and antiperiodic boundary conditions in the 
and 
directions, respectively, to avoid a numerical difficulty caused by the degeneracy between states 
and 
where
. It has been shown that a long-range spin correlation exists in the ground state of the half-filled Hubbard model with 
for 
[40, 41,179,180]. In the case
, 
exhibits a double logarithmic behavior
. 
is shown as a function of 
in Figure 54 for
, 3, and 4. The obtained values are well fitted by 
and 
diverges in the limit of a large system size
:
(104)
This result is consistent with the existence of the long-range spin correlation for 
[179,180]. The degree of divergence of 
is beyond the criterion of the Kosterlitz-Thouless transition, and thus the long-range order represented by 
belongs to a different category. The 
behavior of 
is consistent with the predictions of perturbation theory in the 2D non-linear sigma model at low temperatures [181].
4.5. Pair Susceptibility
In this section we consider a method to evaluate the pair susceptibility 
at
. In order to compute the pair susceptibility, we use an electron-hole transformation for the down spin, 
, whereas the up-spin electrons are unaltered,
. For the on-site s-wave pairing, the source term is given by the following expression
(105)
For the anisotropic d-wave pairing, we add
(106)
where 
and
. The s-wave and d-wave pair susceptibility are respectively:
(107)
Using the Fourier transformation, the source term for the pair potential is written as follows 
for 
or 
with the 
-dependence factor
. If we define
, then for a small value of
, we have the following
(108)
where
(109)
for 
and
. On the basis of analytic continuation, using the thermal Green function, 
is written as
(110)
In the noninteracting system, this formula exhibits logarithmic divergence on the finite lattice
: 
with constants 
and
, which can be confirmed by numerical estimations on finite systems.
In the Kosterlitz-Thouless theory, the susceptibility is scaled as follows [182,183]:
, where 
is the coherence length. 
is of order 
on a lattice 
if long-range coherence exists. The exponent 
is expected to be 0 at absolute zero. Thus 
scales as 
in the ground state if the Kosterlitz-Thouless transition occurs at some temperature.
First, we investigate 
for the attractive Coulomb interaction
. For this model, the existence of a Kosterlitz-Thouless transition has been predicted on the basis of quantum Monte Carlo methods [34,183]. The results in Figure 55 show that the size dependence for 
and 
is
(111)
which is consistent with previous studies, and shows the existence of a Kosterlitz-Thouless transition for the attractive interaction. At near half-filling, 
is more enhanced than that at
. Second, let us investigate the d-wave pair susceptibility 
for the repulsive Coulomb interaction. Pair susceptibilities are
sensitively dependent on the band structure, particularly the energy of the van Hove singularity, as a characteristic of two-dimensional systems. We compute 
at an electron density
, a value near that of optimally doped high-temperature cuprates. We set
. Figure 56 shows 
as a function of 
for
, 3, 4, and 5 with 
and
. This shows that
(112)
if 
is moderately large. This result shows that a d-wave superconducting Kosterlitz-Thouless transition may exist for the repulsive interaction if we adjust the band parameters in the region of optimal doping.
5. Summary
We have investigated the superconductivity of electronic origin on the basis of the (single-band and three-band) two-dimensional Hubbard model. First, we employ the variational Monte Carlo method to clarify the phase diagram of the ground state of the Hubbard model. The superconducting condensation energy per site obtained by the Gutzwiller ansatz is reasonably close to experimental value
. We have examined the stability of striped and checkerboard states in the under-doped region. The relation of the incommensurability and hole density, 
, is satisfied in the
lower doped region. We have found that the 
period checkerboard spin modulation is stabilized in the two-dimensional Hubbard model with the Bi-2212 type band structure.
We have further performed investigation by using the quantum Monte Carlo method that is an exact unbiased method. We have presented an algorithm of the quantum Monte Carlo diagonalization to avoid the negative sign problem in quantum simulations of many-fermion systems. We have computed d-wave pair correlation functions. In the half-filled case 
is suppressed for the repulsive
, and when doped away from half-filling
, 
is enhanced slightly for
. It is noteworthy that the correlation function 
is indeed enhanced and is increased as the system size increases in the 2D Hubbard model. The enhancement ratio increases almost linearly 
as the system size is increased, which is an indication of the existence of superconductivity. Our criterion is that when the enhancement ratio as a function of the system size 
is proportional to a certain power of
, superconductivity will be developed. This ratio depends on 
and is reduced as 
is decreased. The dependence on the band filling shows a dome structure as a function of the electron density. In the 
system, the ratio is greater than 1 in the range
. Let us also mention on superconductivity at half-filling. Our results indicates the absence of superconductivity in the half-filling case because there is no enhancement of pair correlation functions
6. Acknowledgements
We thank I. Hase, S. Koikegami, S. Koike and J. Kondo for stimulating discussions. This work was supported by a Grant-in Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan, and CREST Program of Japan Science and Technology Agency. A part of numerical calculations was performed at the facilities in the Supercomputer Center of the Institute for Solid State Physics of the University of Tokyo.