Correlated-Electron Systems and High Temperature Superconductivity

We present recent theoretical results on superconductivity in correlated-electron systems, especially in the two-dimensional Hubbard model and the three-band d-p model. The mechanism of superconductivity in high-temperature superconductors has been extensively studied on the basis of various electronic models and also electron-phonon models. In this study we investigate the properties of superconductivity in correlated-electron systems by using numerical methods such as the variational Monte Carlo method and the quantum Monte Carlo method. The Hubbard model is one of basic models for strongly correlated electron systems, and is regarded as the model of cuprate high temperature superconductors. The d-p model is more realistic model for cuprates. The superconducting condensation energy obtained by adopting the Gutzwiller ansatz is in reasonable agreement with the condensation energy estimated for YBa_2Cu_3O_7. We show the phase diagram of the ground state using this method. We have further investigated the stability of striped and checkerboard states in the under-doped region. The relationship of the hole density x and incommensurability \delta, \delta\sim x, is satisfied in the lower doping region, as indicated by the variational Monte Carlo calculations. We have performed a variational Monte Carlo simulation on a two-dimensional t-t'-t"-U Hubbard model with a Bi-2212 type band structure and found that the 4\times 4 period checkerboard spin modulation, that is characterized by multi Q vectors, is stabilized. We also present a new algorithm of the diagonalization quantum Monte Carlo method that is a method for the evaluation of expectation values without the negative sign difficulty. We show that the pair correlation function is not enhanced at half-filling, and is indeed enhanced with hole doping.


I. 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][2][3][4][5], heavy fermions [6][7][8][9] and organic conductors [10]. The phase diagram for the typical high-T c cuprates is shown in Fig.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 Fig.2 there is the anomalous metallic region where the susceptibility and 1/T 1 show a peak above T c 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-T c cuprates have the d-wave symmetry in the hole-doped materials. [11,12] Several evidences of d-wave pairing symmetry were provided for the electron-doped cuprates Nd 2−x Ce x CuO 4 . [13][14][15] Thus it is expected that the superconductivity of electronic origin is a candidate for the high-T c superconductivity. We can also expect that the origin of d-wave superconductivity lies in the onsite Coulomb interaction of the two-dimensional Hubbard model.
The antiferromagnetism should also be examined in correlated electron systems.
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 U , that is given by [62][63][64][65][66] where U is the strength of the on-site Coulomb interaction and the exponent x 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. x 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 6 ≤ U ≤ 12, 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 YBa 2 Cu 3 O 7 . 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 δ, δ ∼ x, is satisfied in the lower doping region. This is consistent with the results by neutron scattering measurements. To examine the stability of a 4 × 4 checkerboard state, we have performed a variational Monte Carlo simulation on a two-dimensional t−t ′ −t ′′ −U Hubbard model with a Bi-2212 type band structure. We found that the 4 × 4 period checkerboard checkerboard spin modulation that is characterized by multi Q 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][41][42]. We present the results obtained by a method, quantum Monte Carlo diagonalization, without the negative sign problem.

The Hubbard Hamiltonian is
where c † iσ and c iσ denote the creation and annihilation operators of electrons, respectively, and n iσ = c † iσ c iσ is the number operator. The second term represents the onsite Coulomb interaction which acts when the two electrons occupy the same site. The numbers of lattice sites and electrons are denoted as N and N e , respectively. The electron density is n e = N e /N . In the non-interacting limit U = 0, 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 t ij = 0, each site is occupied by the up-or down-spin electron, or is empty. The non-zero t ij induces the movement of electrons that leads to a metallic state id N e = N . The ground state is probably insulating at half-filling N e = N if U is sufficiently large.
If t ij = t are non-zero only for the nearest-neighbor pairs, the Hubbard Hamiltonian is transformed to the following effective Hamiltonian for large U/t: [67] where a iσ = c iσ (1 − n i,−σ ) and j + µ and j + µ ′ 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 J = 4t 2 /U . 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][71][72] and conformal field theory. [73][74][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 onedimensional Hubbard model. A possibility of superconductivity is a hot topic as well as the magnetism and metal-insulator transition for the two-and threedimensional Hubbard model.
The three-band Hubbard model that contains d and p orbitals has also been investigated intensively with respect to high temperature superconductors. [24,64,[77][78][79][80][81][82][83][84][85][86][87][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 CuO 2 planes which are contained usually in the crystal structures of high-T c superconductors. The network of CuO 2 layer is shown in Fig.3. The parameters of the three-band Hubbard model are given by the Coulomb repulsion U d , energy levels of p electrons ǫ p and d electron ǫ d , and transfer between p orbitals given by t pp . Typical parameter values for the three-band (d-p) Hubbard model are shown in Table I.
The Hamiltonian of the three-band Hubbard model is written as [24,47,48,80] x andŷ represent unit vectors along x and y directions, respectively. p † i±x/2,σ and p i±x/2,σ denote the operators for the p electrons at site R i ±x/2. Similarly p † i±ŷ/2,σ and p i±ŷ/2,σ are defined. U d denotes the strength of Coulomb interaction between d electrons. For simplicity we neglect the Coulomb interaction among p electrons in this paper. Other notations are standard and energies are measured in t dp units. The number of cells is denoted as N for the three-band Hubbard model. In the non-interacting case (U d = 0) the Hamiltonian in the k-space is written as: where d kσ (d † kσ ), p xkσ (p † xkσ ) and p ykσ (p † ykσ ) are operators for d-, p x -and p y -electron of the momentum k and spin σ, respectively.
In the limit t dp where t ef f ≃ t 2 dp /(ǫ p − ǫ d ). J K = 4t ef f gives the antiferromagnetic coupling between the neighboring d and p electrons. In real materials (ǫ p − ǫ d )/t dp is not so large. Thus it seems that the mapping to the t-J model is not necessarily justified.

III. VARIATIONAL MONTE CARLO STUDIES
In this Section we present studies on the twodimensional Hubbard model by using the variational where ij denotes summation over all the nearestneighbor bonds and jℓ means summation over the next nearest-neighbor pairs. t is our energy unit. The dispersion is given by Our trial wave function is the Gutzwiller-projected wave functions defined as where P G is the Gutzwiller projection operator given by g is a variational parameter in the range from 0 to unity and j labels a site in the real space. P Ne is a projection operator which extracts only the sites with a fixed total electron number N e . Coefficients u k and v k in ψ BCS appear in the ratio defined by where ξ k = ǫ k − µ and ∆ k is a k-dependent gap function. µ is a variational parameter working like the chemical potential. c kσ is the Fourier transform of c jσ . The wave functions ψ n and ψ s are expressed by the Slater determinants for which the expectations values are evaluated using a Monte Carlo procedure. [43,44,92] ψ s is written as where a(j, ℓ) = (1/N ) Then ψ s is written using the Slater determinants as where A(j 1 · · · i Ne/2 , ℓ 1 · · · ℓ Ne/2 ) is the Slater determinant defined by A(j 1 · · · ℓ Ne/2 ) = a(j 1 , ℓ 1 ) a(j 1 , ℓ 2 ) · · · a(j 1 , ℓ Ne/2 ) a(j 2 , ℓ 1 ) a(j 2 , ℓ 2 ) · · · a(j 2 , ℓ Ne/2 ) · · · · · · · · · · · · a(j Ne/2 , ℓ 1 ) a(j Ne/2 , ℓ 2 ) · · · a(j Ne/2 , ℓ Ne/2 ) .
In the process of Monte Carlo procedure the values of cofactors of the matrix in eq.(18) are stored and corrected at each time when the electron distribution is modified. We optimized the ground state energy with respect to g, ∆ k and µ for ψ s for ψ s . For ψ n the variational parameter is only g. We can employ the correlated measurements method [93] in the process of searching optimal parameter values minimizing E g . 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][95][96] Note that the Gutzwiller projection operator is written as where α = log(1/g). Then using the discrete Hubbard-Stratonovich transformation, the Gutzwiller operator is the bilinear form: where cosh(2a) = e α/2 . The Hubbard-Stratonovich auxiliary field s i takes the values of ±1. The norm ψ n |ψ n is written as where V σ (s, α) is a diagonal N ×N matrix corresponding to the potential V σ (s, α) is written as where diag(a, · · · ) denotes a diagonal matrix with elements given by the arguments a, · · · . The elements of (φ σ 0 ) ij (i = 1, · · · , N ; j = 1, · · · , N e /2) are given by linear combinations of plane waves. For example, 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 III.J.

B. Superconducting Condensation Energy
We study the cases of the d-, extended s-(s * -) and s-wave gap functions in the following: In t d dependence of the SC condensation energy ∆E/2N for the two-chain Hubbard model in the bulk limit. [54] 0. [46] E g /N is plotted as a function of ∆ for three types of gap functions shown above. We impose the periodic and the antiperiodic boundary conditions for x-and ydirection, respectively. This set of boundary conditions is chosen so that ∆ k does not vanish for any k-points occupied by electrons. E g was obtained as the average of the results of several Monte Carlo calculations each with 5 × 10 7 steps. E g /N has minimum at a finite value of ∆ ≃ 0.08 in the case of the d-wave gap function.
The energy gain ∆E g in the superconducting state is called the SC condensation energy in this paper. ∆E g /N is plotted as a function of 1/N in Fig.5 in order to examine the size dependence of the SC energy gain. [97] Lattice sizes treated are from 8 × 8 to 22 × 22. The electron  density n e is in the range of 0.80 ≥ n e ≤ 0.86. Other parameters are −0.20 ≤ t ′ ≤ 0.0 and U/t = 8 in t units. Bulk limit N → ∞ of SC condensation energy E cond was obtained by plotting as a function of 1/N . The linear fitting line indicates very clearly that the bulk limit remains finite when −0.25 ≤ t ′ ≤ −0.10 and n e ≥ 0.84.
When n e = 0.86, t ′ = −0.20 and U = 8, the bulk-limit E cond is E cond = 0.00117/site ≃ 0.60 meV/site, where t = 0.51eV is used. [ . 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 E cond = 0.026t ≃ 13 meV/site at n e = 0.84 for J = 4t 2 /U = 0.5 and t ′ = 0. [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 t/U 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 U in Fig.6. The condensation energy E cond = ∆E/N is increased as U/t is increased as far as U/t ≤ 12. In the strong coupling region U > 8t, we obtain the large condensation energy.

C. 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) 2 CuO 4 (LSCO) is estimated at 0.029meV/(Cu site) or 0.00008 in units of t which is much smaller than that for YBCO. [103] The band parameter values of LSCO were estimated as t ′ = −0.12 and t ′′ = 0.08. [104] This set corresponds roughly to E cond ≃ 0.0010. 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 E cond is much reduced. Therefore, such a coexistence allows us to qualitatively understand the SC E cond in LSCO.
On the other hand, Tl2201 (T c = 93K) and Hg1201 (T c = 98K) band calculations by Singh and Pickett [107] give very much deformed Fermi surfaces that can be fitted by large |t ′ | such as t ′ ∼ −0.4. For Tl2201, an Angular Magnetoresistance Oscillations (AMRO) work [108] gives information of the Fermi surface, which allows to get t ′ ∼ −0.2 and t ′′ ∼ 0.165. 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 t ′ ∼ −0.2 and t ′′ ∼ 0.175. For such a deformed Fermi surface, E cond in the bulk limit is reduced considerably. [111,112] Therefore, the SC E cond calculated by VMC indicates that the Fermi surface of LSCO-type is more favorable for high T c . The lower T c in LSCO may be attributed to the coexistence with antiferromagnetism of stripe type.

D. 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][114][115][116] where c † jℓσ (c jℓσ ) is the creation (annihilation) operator of an electron with spin σ at the ℓth site along the jth chain (j = 1, 2). t is the intrachain nearest-neighbor transfer and t d is the interchain nearest-neighbor transfer energy. The energy is measured in t units. The energy minimum was given when the components of the gap function ∆ k take finite values plotted in Fig.7 for the lattice of 20 × 2 sites with 34 electrons imposing the periodic boundary condition. [54] Each component of ∆ k was optimized for U 0 = 8 and t d = 1.8. 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 ∆ k of the antibonding band (k y = π) is much larger than that of the bonding band (k y = 0). In order to reduce the computation cpu time, ∆ k of each band was forced to take a fixed value specific to each band, i.e. ∆ 1 for the bonding band and ∆ 2 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. ∆ 1 and ∆ 2 take opposite sign, which is similar to that of the d x 2 −y 2 gap function.
The energy gain ∆F 2c remains finite in the bulk limit when 1.2 < t d < 1.6. The SC condensation energy per site in the bulk limit is plotted as a function of t d in Fig.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 ∆E 2c is 0.0008 which is of the same order of magnitude as the maximum condensation energy obtained for the 2D Hubbard model. [46]

E. Condensation Energy in the d-p Model
The SC energy gain for the d-p model, namely, threeband Hubbard model in eq.(4) has also been evaluated using the variational Monte Carlo method. For the threeband model the wave functions are written as where α kσ is the linear combination of d kσ , p xkσ and p ykσ 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 g, effective level differenceǫ p −ǫ d , chemical potential µ and superconducting order parameter ∆ are the variational parameters.
The similar results to the single-band Hubbard model were obtained as shown in Fig.9 for t pp = 0.0, U d = 8 and ǫ d − ǫ p = 2 in t dp units where the calculations were performed in the hole picture. [24] The SC condensation energy for the three-band model is E cond ≃ 0.0005t dp ≃ 0.75 meV per site in the optimally doped region. We set t dp = 1.5 eV as estimated in Table I. There is a tendency that E cond increases as ǫ d − ǫ p increases which is plotted in Fig.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 III.K.

F. 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 Summation over k and k ′ in eq.(33) is performed over the filled k-points, as in the calculation of the normal state energy. Q is the SDW wave vector given by (π, π) and ∆ AF is the SDW potential amplitude. As shown in Fig.11, the energy gain per site in the SDW state rises very sharply from n e ∼ 0.84 for t ′ = 0. [46] At n e ∼ 0.84 it is slightly larger than that in the SC state, and at n e = 0.80 there is no more stable SDW state. Thus the boundary between the SDW and the SC states is given at n e ∼ 0.84. The results of the bulk limit calculations indicate that the energy gain in the SC state at n e = 0.84 takes the extremely small value and the value at n e = 0.80 vanishes as N → 0. Hence the pure d-wave SC state possibly exists near the boundary at n e ∼ 0.84, 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 t pp = 0.0 and ǫ p − ǫ d = 2 in the plane of U and the hole density in Fig.12 where AF denotes the antiferromagnetic region. [47] The value ǫ p −ǫ d = 2 is taken from the estimations by cluster calculations. [89][90][91] The phase boundary Since the three-band Hubbard model contains several parameters, we must examine the parameter dependence of the energy of SDW state. The energy gain ∆E SDW in the SDW state is shown in Fig.13 as a function of doping ratio for several values of ∆ dp ≡ ǫ p − ǫ d . ∆E SDW increases as ∆ dp increases as expected. In Fig.14   U d -dependencies of ∆E SDW are presented. The SDW phase extends up to 30 percent doping when U d is large. It follows from the calculations that the SDW region will be reduced if ǫ p − ǫ d and U d decrease.
From the calculations for the SDW wave functions, we should set ǫ p − ǫ d and U d small so that the SDW phase does not occupy a huge region near half-filling. In tivity in this region suggests that the electric conduction is approximately one dimensional. [124] The angleresolved photo-emission spectroscopy measurements also suggested a formation of two sets of one-dimensional 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 0.05 < x < 0.
while for the open circles ∆ i,i+x = ∆cos|(Qx(xi +x/2))| and possibility of superconductivity that occurs in the stripe state is a subject of great interest. [127][128][129][130] The incommensurate magnetic scattering spots around (π, π) were observed in the SC phase in the range of 0.05 < x < 0.13 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 neutron-scattering experiments Phase diagram of the d-p model based on the Gutzwiller wave function. [106] have also revealed that a diagonal spin modulation occurs across the insulating spin-glass phase in La 2−x Sr x CuO 4 for 0.02 ≤ x ≤ 0.05, 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 Fig.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][132][133][134][135][136][137][138][139][140][141] The stripe SDW state has been studied theoretically by using the meanfiled theory. [132][133][134][135][136] They found that the stripe state appears when an incommensurate nesting becomes 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 where φ MF coexist is a mean-field wave function. The effective mean-field Hamiltonian for the coexisting state is [105] represented by where the diagonal terms describe the incommensurate spin-density wave state: where µ is the chemical potential. The vertical stripe state is represented by the charge density n i and the spin density m i that are spatially modulated as where Y ℓ denotes the position of vertical stripes. The amplitude α is fixed by i n i = N e . The off-diagonal terms in eq. (38) are defined in terms of the d-wave SC gap as whereê = ±x, ±ŷ (unit vectors). We consider two types of the spatially inhomogeneous superconductivity: antiphase and in-phase defined as and respectively. Here, q = (0, 2πδ) and δ is a incommensurability given by the stripe's periodicity in the y direction with regard to the spin. The anti-phase (in-phase) means that the sign if the superconducting gap is (is not) changed between nearest domain walls. The wave function ψ 0 coexist is made from the solution of the Bogoliubov-de Gennes equation represented by The Bogoliubov quasiparticle operators are written in the form Then the coexistence wave function is written as [105,142] ψ 0 for constants C and C ′ . The calculations are performed for the wave function ψ coexist = P G ψ 0 coexist . The variational parameters are µ, m, g, ξ c and ξ s . The system parameters were chosen as t ′ = −0.20 and U = 8 suitable for cuprate superconductors. It has been shown that the "anti-phase" configuration is more stable than the "inphase" one. [105] Here, the system parameters are t ′ = −0.2 and U = 8. We use the periodic boundary condition in the xdirection and anti-periodic one in the y-direction. In Fig.19, we show the total energy of the coexistent state, E coexist , as a function of the SC gap ∆ for the cases of anti-phase and in-phase. The SC condensation energy ∆E coexist is estimated as 0.0008t per site at the hole density 0.125 on the 12 × 8 lattice with periodic boundary condition in x-direction and antiperiodic one in y-direction. ∆E coexist in the coexistence state is defined as the decrease of energy due to finite ∆. If we use t ∼ 0.5eV, this is evaluated as ∼ 0.4meV. The SC condensation energy per site is shown as a function of hole density in Fig.20. One finds that ∆E coexist 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 t ′ in the trial function as a variational parametert ′ . In Fig.21, optimized values oft ′ /t for the striped state are shown as a function of U/t. Thet ′ /t increases as U/t increases. We also mention that the optimized t ′′ /t almost vanishes. The renormalized Fermi surface oft ′ /t = −0.30, −0.40 and −0.50 are plotted in Fig.22. The system is a 16 × 16 lattice with t ′ /t = −0.2 and the electron density 0.875. As U/t is incresed, the Fermi surface is more deformed. We show the the gradient of the momentum distribution function, |∇n k |, calculated in the optimized stripe state in Fig.23. The brighter areas coincide with the renormalized Fermi surface witht ′ /t = −0.31 andt ′′ /t = 0.0 for U/t = 8.
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 Fig.24 presents the region of stable AF phase in the plane of t pp and ∆ dp = ǫ p − ǫ d . For large ∆ dp = ǫ p − ǫ d , 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 x = 1/16 is shown in Fig.25 where the 4-lattice stripe state has higher energy than that for 8-lattice stripe for all the values of t pp .
The Bogoliubov-de Gennes equation is extended to the case of three orbitals d, p x and p y . (H ijσ ) and (F ij ) are now 3N ×3N matrices. The energy of the state with double order parameters ∆ and m is shown in Fig.26 on the 16×4 lattice at the doping rate 1/8. The SC condensation energy per site is evaluated as ∼ 0.00016t dp for U d = 8, t pp = 0.4 and ǫ p − ǫ d = 2. If we use t dp ∼ 1.5eV, [89][90][91] we obtain ∆E coexist ∼ 0.24meV 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 energy for x = 0.2, 0.125 0.08333 and 0.0625 in Fig.27. We set the parameters as ǫ p − ǫ d = 2 and t pp = 0.4 in t dp units, which is reasonable from the viewpoint of the density of states and is remarkably in accordance with cluster estimations [89][90][91], and also in the region of eight-lattice periodicity at x = 1/8. We have carried out the Monte Carlo calculations up to 16 × 16 sites (768 atoms in total). In the overdoped region in the range of 0.18 < x < 0.28, we have the uniform d-wave pairing state as the ground state. The periodicity of spatial variation increases as the doping rate x decreases proportional to 1/x. In the figure we have the 12-lattice periodicity at x = 0.08333 and 16-lattice periodicity at x = 0.0625. For x = 0.2, 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 5.0×10 7 ) is required to get convergent expectation values for each set of parameters.

H. Diagonal Stripe States in the Light-Doping Region
Here we examine whether the relationship δ ∼ x 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, 0.03 < x < 0.07, revealed that the position of incommensurate magnetic peaks changed from (1/2, 1/2 ± δ) to (1/2 ± δ ′ , 1/2 ± δ ′ ) as x 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 δ ∼ x that holds in the vertical stripe state in the low doping region.
In Fig.29, we show the incommensurability of the most stable stripe state as a function of x. 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 x critical between the diagonal and vertical stripe states lies at or above 0.0625 in the case of U = 8 and t ′ = −0.2. The following factors may give rise to slight changes of the phase boundary x critical : the diagonal stripe state may be stabilized in the lowtemperature-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 t ′ generated by the Cu-O buckling in the LTO phase.

I. Checkerboard States
A checkerboard-like density modulation with a roughly 4a × 4a period (a is a lattice constant) has also been observed by scanning tunneling microscopy (STM) experiments in Bi 2 Sr 2 CaCu 2 O 8+δ , [146] Bi 2 Sr 2−x La x CuO 6+δ , [147] and Ca 2−x Na x CuO 2 Cl 2 (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 ρ i and spin density m i are spatially modulated as where ρ ℓ and m ℓ 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 Q s 1 = (π, π ± 2πδ) and Q s 2 = (π ± 2πδ, π) describe a spin vertical checkerboard state, where two diagonal domain walls are orthogonal. Diagonal wave vectors Q s 1 = (π ± 2πδ, π ± 2πδ) and Q s 2 = (π ± 2πδ, π ∓ 2πδ) lead to a spin diagonal checkerboard state with a 1/δ-period. The hole density forms the charge vertical checkerboard pattern with vertical wave vectors Q c 1 = (0, ±4πδ) and Q c 2 = (2π ± 4πδ, 2π) in which the hole density is maximal on the crossing point of magnetic domain walls in the spin diagonal checkerboard state. If δ = 1/8 is assumed, the charge modulation pattern is consistent with the 4a × 4a 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 4 × 4 period checkerboard spin modulation. [151] In Fig.30(a), we show the condensation energies of some heterogeneous states, (E normal − E hetero )/N site , with fixing the transfer energies t ′ = −0.32 and t ′′ = 0.22 suitable for Bi-2212. The system is a 16 × 16 lattice with the electron-filling ρ = N e /N site = 0.875. The energy of the normal state E normal is calculated by adopting m ℓ = ρ ℓ = 0. In our calculations, the condensation energies of both bond-centered 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 spinvertical checkerboard (VC) with ρ ℓ = 0 is most stable, as shown in Fig.30(a). The measured expectation value of the spin density is shown in Fig.30(b).

J. 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] ψ Jastrow = jℓ( =0) σσ ′ [1 − (1 − g(ℓ))n jσ n j+ℓσ ′ ]ψ n , (54) for the variational parameter g(ℓ). We have investigated the 2D Hubbard model by using the Jastrow-Gutzwiller 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 ψ (m) defined in the following way: [95] and so on, where K denotes the kinetic part of the Hamiltonian and V denotes the on-site Coulomb interaction. λ, λ ′ , λ ′′ ,α, α ′ and α ′′ are variational parameters. It is considered that ψ (m) approaches the true ground state wave function as m grows larger since the ground state wave function is written as for large β = ǫ 1 + · · · + ǫ m and small ǫ i (i = 1, · · · , m). If we can extrapolate the expectation values from the data for ψ (1) , ψ (2) , · · · , we can evaluate correct expectation values. The computations are performed by using the discrete Hubbard-Stratonovich transformation as described in Section III.A. In the evaluation of the expectation values we generate the Monte Carlo samples by the importance sampling [95] with the weight function |w| = |w ↑ w ↓ | where Since the Monte Carlo samplings are generated with the weight |w|, the expectation values are calculated with the sign of w 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 Fig.31 the energy is shown as a function of 1/m where the SDW and normal states are chosen as the initial state ψ 0 . The extrapolated values for different initial states coincide with each other within Monte Carlo statistical errors. The energy expectation values as a function of U for 8 × 8 square lattice are presented in Fig.32 for ψ n = ψ (1) , ψ AF , ψ (3) . 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 U = 8. 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 [155]: The up-spin electrons are unaltered.
We show the energy versus 1/m in Fig.33 for ψ (m) and ψ (m) s . From an extrapolation to the limit m → ∞, 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  [96] lower than that in the normal state for each m. The energy gain in the SC state remains the same order after the multiplication of correlation factors.

K. Tc and ǫp − ǫ d
Relationships between T c and structural features in cuprate high-temperature superconductors are very interesting. Torrance and Metzger found the first such relationship between T c and the Madelung potential difference ∆V M . [156] Here ∆V M is the potential difference between Cu and O sites in the CuO 2 plane. T c was found to increase with decreasing ∆V M . There is an interesting tendency of increasing T c 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][63][64][65][66]. There have been many similar works by making some kind of approximation such as Random phase approximation (RPA) [157][158][159], fluctuation-exchange approximation (FLEX) [160][161][162][163], effective spin-fluctuations [4,164,165], and perturbation theory in terms of U [166][167][168]. An application was made for Sr 2 RuO 4 where we need to consider the multi-band structure including α and β orbitals [169], and also to the three-dimensional d-p model [170]. In our formulation the gap function is written as The exponent x indicates the strength of superconductivity. The results are in Figs. 34 and 35. [171] As shown in the figure, for positive ǫ p − ǫ d , with increase of ǫ p − ǫ d the exponent x increases monotonously. This means the increase of superconducting gap and so of T c , and is consistent with the wide-range tendency of the variational Monte Carlo calculation. [24,172] This tendency can be understood in terms of where u 0 k is the weight ofd electrons. This clearly indicates that increase of ǫ p − ǫ d leads to the increase of U ef f and subsequently of x. In the case of ǫ p − ǫ d < 0, we take account of finite Coulomb repulsion U p on oxygen sites. The effective interaction coming from U p is similarly given by the susceptibility with the weight of p electrons. The results of x with U p /U d = 6/8 indicates that all four types of even parity (b 1g , b 2g , a 1g and a 2g ) SC strength values increase, so that T c is raised, as the absolute value |ǫ p − ǫ d | increases in this region. This result shows that U p also plays an important role as well.
Let us give a discussion on this result. Increase of |ǫ p − ǫ d | in the region of ǫ p − ǫ d < 0 means decrease of is the second electron affinity of oxygen atom and I Cu 3 is the third electron ionicity of copper atom and e 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 CuO 2 plane, cuprate superconductors are known to increase T c . [173] The accompanying raise of ǫ d should tend to increase T c .
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 ǫ p or the raise of ǫ d relatively. This indicates that T c will be increased by the Coulomb interaction between p electrons.

A. 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 where c † jσ (c jσ ) is the creation (annihilation) operator of an electron with spin σ at the j-th site and n jσ = c † jσ c jσ . t ij is the transfer energy between the sites i and j. t ij = t for the nearest-neighbor bonds. For all other cases t ij = 0. U is the on-site Coulomb energy. The number of sites is N and the linear dimension of the system is denoted as L. The energy unit is given by t and the number of electrons is denoted as N e .
In a Quantum Monte Carlo simulation, the ground state wave function is where ψ 0 is the initial one-particle state represented by a Slater determinant. For large τ , e −τ H will project out the ground state from ψ 0 . We write the Hamiltonian as H = K +V where K and V are the kinetic and interaction terms of the Hamiltonian in Eq.(64), respectively. The wave function in Eq. (65) is written as for τ = ∆τ · M . Using the Hubbard-Stratonovich transformation [27,94], we have for (tanha) 2 = tanh(∆τ U/4) or cosh(2a) = e ∆τ U/2 . The wave function is expressed as a summation of the oneparticle Slater determinants over all the configurations of the auxiliary fields s j = ±1. The exponential operator is expressed as where we have defined for The ground-state wave function is where φ m is a Slater determinant corresponding to a configuration m = {s i (ℓ)} (i = 1, · · · , N ; ℓ = 1, · · · , M ) of the auxiliary fields: The coefficients c m are constant real numbers: c 1 = c 2 = · · · . The initial state ψ 0 is a one-particle state. If electrons occupy the wave numbers k 1 , k 2 , · · · , k Nσ for each spin σ, ψ 0 is given by the product ψ ↑ 0 ψ ↓ 0 where ψ σ 0 is the matrix represented as [31]     e ik1·r1 e ik2·r1 · · · · · · e ikN σ ·r1 e ik1·r2 e ik2·r2 · · · · · · · · · · · · · · e ik1·rN e ik2·rN · · · · · ·     .

(74)
N σ is the number of electrons for spin σ. In actual calculations we can use a real representation where the matrix elements are cos(k i · r j ) or sin(k i · r j ). In the real-space representation, the matrix of V σ ({s i }) is a diagonal matrix given as φ σ m is an N × N σ matrix given by the product of the matrices e −∆τ Kσ , e Vσ and ψ σ 0 . The inner product is thereby calculated as a determinant [38], The expectation value of the quantity Q is evaluated as If Q is a bilinear operator Q σ for spin σ, we have The expectation value with respect to the Slater determinants φ σ m Q σ φ σ n is evaluated using the single-particle Green's function [31,38], In the above expression, 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 |P mn |; the resulting relationship is, where sign(a) = a/|a| and 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, where n MC is the number of samples. The sign problem is an issue if the summation of sign(P mn ) vanishes within statistical errors. In this case it is indeed impossible to obtain definite expectation values.

Basic Method and Optimization
Quantum Monte Carlo diagonalization (QMD) is a method for the evaluation of Q σ without the negative sign problem. [41] A bosonic version of this method was developed before in Ref. [174]. The configuration space of the probability P mn in Eq.(84) is generally very strongly peaked. The sign problem lies in the distribution of P mn in the configuration space. It is important to note that the distribution of the basis functions φ m (m = 1, 2, · · · ) is uniform since c m are constant numbers: c 1 = c 2 = · · · . In the subspace {φ m }, selected from all configurations of auxiliary fields, the right-hand side of Eq.(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, where c m (m = 1, 2, · · · ) are variational parameters. In order to minimize the energy the equation ∂E/∂c n = 0 (n = 1, 2, · · · ) is solved for, If we set the eigen-equation is for u = (c 1 , c 2 , · · · ) t . Since φ m (m = 1, 2, · · · ) are not necessarily orthogonal, A is not a diagonal matrix. We diagonalize the Hamiltonian A −1 H, 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 M is large enough, the wave function in Eq.(72) will approach the exact groundstate wave function, ψ exact , as the number of basis functions, N states , is increased. If the number of basis functions is large enough, the wave function will approach, ψ exact , as M 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]: where v denotes the variance defined as and Q exact is the expected exact value of the quantity Q. The simplest procedure for optimizing the groundstate wave function is to increase the number of basis states {φ m } by random sampling. First, we set τ and M , for example, τ = 0.1, 0.2, · · · , and M = 20, 30, · · · . We denote the number of basis functions as N states . We start with N states = 100 ∼ 300 and then increase up to 10000. This procedure can be outlined as follows: A1. Generate the auxiliary fields s i (i = 1, · · · , N ) in B σ ℓ ({s i })) randomly for ℓ = 1, · · · , M for φ m (m = 1, · · · , N states ), and generate N states basis wave function {φ m }.
A2. Evaluate the matrices H mn = φ m Hφ n and A mn = φ m φ n , and diagonalize the matrix A −1 H to obtain ψ = m c m φ m . Then calculate the expectation values and the energy variance.

A3.
Repeat the procedure from A1 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 {s i (ℓ)} (i = 1, · · · , N ; ℓ = 1, · · · , M ) in φ n that has the large weight |c n | 2 to generate a new basis function φ ′ n . The new basis function φ ′ n 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].

Ground State Energy and Correlation Functions
The energy as a function of the variance is presented in Figs.36, 37 and 38 for 4 × 4, 6 × 2 and 6 × 6, respectively. To obtain these results the genetic algorithm was employed to produce the basis functions except the open symbols in Fig.37. The 4 × 4 where N e = 10 in Fig.36 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 4 × 4, 6 × 2 and 6 × 6 systems in Table II.
The Fig. 39 is the momentum distribution function n(k), for 14×14 sites where the results for the Gutzwiller VMC and the QMD are indicated. The Gutzwiller function gives the results that n(k) increases as k approaches k F 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.

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 U > 0. The ground state is insulating at half-filling and metallic for less than half-filling. The Fig.  40 is the spin and charge correlation functions, S(k) and C(k), as a function of the wave number, for the 1D Hubbard model where N = 80. The 2k F singularity can be clearly identified where the dotted line is for U = 0. The  [38]. The column VMC is the results obtained for the optimized variational wave function ψ (2) λ except for 6 × 2 for which ψ (1) λ is employed. The QMC results are from Ref. [35]. Exact results are obtained using diagonalization [177].
Size The spin correlation function S(k) for the Hubbard ladder is presented in Fig.41, where U = 4 and t d = 1. S(k) is defined as where R iℓ denotes the site (i, ℓ) (ℓ=1,2). We use the convention that k = (k, k y ) where k y = 0 and π indicate the lower band and upper band, respectively. There are four singularities at 2k F 1 , 2k F 2 , k F 1 −k F 2 , and k F 1 +k F 2 for the Hubbard ladder, where k F 1 and k F 2 are the Fermi wave numbers of the lower and upper band, respectively.
It has been expected that the charge gap opens up as U turns on at half-filling for the Hubbard ladder model. In Fig.42 the charge gap at half-filling is shown as a function of U . The charge gap is defined as where E(N e ) is the ground state energy for the N e electrons. The charge gap in Fig.42 was estimated using the extrapolation to the infinite system from the data for the 20 × 2, 30 × 2, and 40 × 2 systems. The data suggest the exponentially small charge gap for small U or the existence of the critical value U c in the range of 0 ≤ U c < 1.5, below which the charge gap vanishes.

Magnetization in 2D Hubbard Model
The ground state of the 2D Hubbard model at halffilling is antiferromagnetic for U > 0 because of the nesting due to the commensurate vector Q = (π, π). The Gutzwiller function predicts that the magnetization increases rapidly as U increases and approaches m = 1 for large U . In Fig.43 the QMD results are presented for m as a function of U . The previous results obtained using

C. Pair Correlation Function
The pair correlation function D αβ is defined by where ∆ α (i), α = x, y, denote the annihilation operators of the singlet electron pairs for the nearest-neighbor sites:  [28]. The squares are the Gutzwiller-VMC results [43] and gray solid circles show the 3rd λ-function (ψ λ ) VMC results carried out on the 8 × 8 lattice [95]. The diamond symbol is the value from the twodimensional Heisenberg model where m = 0.615 [179,180].
Hereα is a unit vector in the α(= x, y)-direction. We consider the correlation function of d-wave pairing: where i and i + ℓ 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 D yy and D yx on the lattice  over the superconductive pairing correlation and thus the pairing correlation function is suppressed as the Coulomb repulsion U is increased. The Fig.48 exhibits the d-wave pairing correlation function P d on 8 × 8 lattice as a function of the distance. The P d 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 U = 4. We then turn to the case of less than half-filling. We show the results on 8 × 8 with electron number N e = 54. We show P d as a function of the distance in Fig.49 (N e = 54). In the scale of this figure, P d for U > 0 is almost the same as that of the non-interacting case, and is enhanced slightly for large  U . Our results indicate that the pairing correlation is not suppressed and is indeed enhanced by the Coulomb interaction U , and its enhancement is very small.
In the low density case, the enhancement is also suppressed being equal to 1. In Fig.52, the enhancement ratio is shown as a function of the electron density n e for U = 4. A dome structure emerges even in small systems. The square in Fig.52 indicates the result for the half-filled  The circles indicate the results for ne ∼ 0.8, where we use the periodic boundary condition in both the x and y directions, and the chemical potential is set at the center of the level spacing between adjacent energy levels. The lowest dotted line is for U = 0 (ne ∼ 0.75), which is fitted by a logarithmic curve, that is, χ s pair ∼ log(L). We show χ s pair for ne ∼ 0.9 and U = −4 by squares, where the boundary condition is antiperiodic in one direction and periodic in the other direction. case with t ′ = −0.2 on 8×8 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 t ′ < 0 enhances P d compared to the case with t ′ = 0 on 8 × 8 lattice. P d is, however, not enhanced over the non-interacting case at half-filling. This also holds for 10 × 10 lattice where the enhancement ratio ∼ 1. This indicates the absence of superconductivity at half-filling.

D. Spin Susceptibility
We proposde a method to compute the magnetic susceptibility at absolute zero (T = 0). [178] We add the source term H 1 to the Hamiltonian as follows j↑ c j↓ e iq·Rj + h.c. = g(S + −q + h.c.), (101) where g is a small real number of the order 10 −3 or 10 −4 . We calculate − S − q /g in the ground state, which is, as shown by the linear response theory, the magnetic susceptibility in the limit of small g, where G ret is the retarded Green function and χ −+ (q, ω) is the dynamical susceptibility, Indeed, the above formula gives the correct spin susceptibility χ −+ (q, ω = 0) on the finite lattice for the noninteracting case, which is given by k (f (ξ k+q ) − f (ξ k ))/(ξ k − ξ k+q ) with the Fermi distribution function f (ξ). We calculate − S − q /g by using the quantum Monte Carlo method to obtain χ −+ (q, ω = 0).
We examine the results obtained for the susceptibilities. Figure 53 shows the spin susceptibility χ(Q) = χ −+ (Q, ω = 0) for Q = (π, π) on a 6 × 2 lattice as a function of 1/m or the energy variance v. 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 χ −+ stag by adding the source term H 1 = g j (c † j↑ c j↓ (−1) jx+jy + h.c.) to the Hamiltonian, where j = (j x , j y ). Here we set the periodic and antiperiodic boundary conditions in the x and y directions, respectively, to avoid a numerical difficulty caused by the degeneracy between states k and k + Q where Q = (π, π). It has been shown that a long-range spin correlation exists in the ground state of the half-filled Hubbard model with t ′ = 0 for U > 0 [40,41,179,180]. In the case U = 0, χ −+ stag exhibits a double logarithmic behavior (log(L)) 2 . χ −+ stag is shown as a function of L in Fig.54 for U = 2, 3, and 4. The obtained values are well fitted by L 4 and χ −+ stag diverges in the limit of a large system size L: This result is consistent with the existence of the longrange spin correlation for U > 0 [179,180]. The degree of divergence of χ −+ stag is beyond the criterion of the Kosterlitz-Thouless transition, and thus the long-range order represented by χ −+ stag belongs to a different category. The L 4 behavior of χ stag is consistent with the predictions of perturbation theory in the 2D non-linear sigma model at low temperatures [181].

E. Pair Susceptibility
In this section we consider a method to evaluate the pair susceptibility χ pair at T = 0. In order to compute the pair susceptibility, we use an electron-hole transformation for the down spin, c i↓ = d † i , whereas the up-spin electrons are unaltered, c i↑ = c i . For the on-site s-wave pairing, the source term is given by the following expression where a ±x = 1 and a ±y = −1. The s-wave and d-wave pair susceptibility are respectively: (a µ c † i+µ d i + h.c.). (107) Using the Fourier transformation, the source term for the pair potential is written as follows H a 1 = g k z k (c † −k↓ c † k↑ + h.c.) for a = s or d with the kdependence factor z k . If we define ∆ k = c k↑ c −k↓ , then for a small value of g, we have the following where for b k = c k↑ c −k↓ and b k (t) = e iHt b k e −iHt . On the basis of analytic continuation, using the thermal Green function, ∆ k is written as In the noninteracting system, this formula exhibits logarithmic divergence on the finite lattice L × L: χ pair = A |z k | 2 log(cL) with constants A and c, which can be confirmed by numerical estimations on finite systems. In the Kosterlitz-Thouless theory, the susceptibility is scaled as follows [182,183]: χ ∼ ξ 2−η , where ξ is the coherence length. ξ is of order L on a lattice L × L if long-range coherence exists. The exponent η is expected to be 0 at absolute zero. Thus χ scales as χ ∼ L 2 in the ground state if the Kosterlitz-Thouless transition occurs at some temperature.
First, we investigate χ s pair for the attractive Coulomb interaction U < 0. 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 Fig.55 show that the size dependence for t ′ = 0 and n e ∼ 0.8 is which is consistent with previous studies, and shows the existence of a Kosterlitz-Thouless transition for the attractive interaction. At near half-filling, χ s pair is more enhanced than that at n e ∼ 0.8. Second, let us investigate the d-wave pair susceptibility χ d pair 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 twodimensional systems. We compute χ d pair at an electron density n e ∼ 0.87, a value near that of optimally doped high-temperature cuprates. We set t ′ = −0.2. Figure 56 shows χ d pair as a function of L 2 for U = 2, 3, 4, and 5 with t ′ = −0.2 and n e ∼ 0.87. This shows that if U 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.

V. SUMMARY
We have investigated the superconductivity of electronic origin on the basis of the (single-band and threeband) 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 0.17∼ 0.26meV/site. We have examined the stability of striped and checkerboard states in the underdoped region. The relation of the incommensurability and hole density, δ ∼ x, is satisfied in the lower doped region. We have found that the 4×4 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 P d is suppressed for the repulsive U > 0, and when doped away from half-filling N e < N , P d is enhanced slightly for U > 0. It is noteworthy that the correlation function P d is indeed enhanced and is increased as the system size increases in the 2D Hubbard model. The enhancement ratio increases almost linearly ∝ L 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 L is proportional to a certain power of L, superconductivity will be developed. This ratio depends on U and is reduced as U is decreased. The dependence on the band filling shows a dome structure as a function of the electron density. In the 10 × 10 system, the ratio is greater than 1 in the range 0.3 < n e < 0.9. Let us also mention on superconductivity at half-filling. Our result indicates the absence of superconductivity in the half-filling case because there is no enhancement of pair correlation functions.
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.