ABSTRACT

We present recent theoretical results on superconductivity in correlated-electron systems, especially in the twodimensional 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₂Cu₃O₇. 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. Holes doped in a half-filled square lattice lead to an incommensurate spin and charge density wave. The relationship of the hole density and incommensurability, , is satisfied in the lower doping region, as indicated by the variational Monte Carlo calculations for the two-dimensional Hubbard model. 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. We have performed a variational Monte Carlo simulation on a two-dimensional --- Hubbard model with a Bi-2212 type band structure and found that the period checkerboard spin modulation, that is characterized by multi Q vectors, is indeed stabilized. We have further performed an investigation by using a quantum Monte Carlo method, which is a numerical method that can be used to simulate the behavior of correlated electron systems. We present a new algorithm of the quantum Monte Carlo diagonalization that is a method for the evaluation of expectation value without the negative sign problem. We compute pair correlation functions and show that pair correlation is indeed enhanced with hole doping.

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-T_c 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 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 -wave symmetry in the hole-doped materials [11,12]. Several evidences of -wave pairing

symmetry were provided for the electron-doped cuprates Nd_2–xCe_xCuO₄ [13-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 -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 YBa₂Cu₃O₇. 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 CuO₂ planes which are contained usually in the crystal structures of high- superconductors. The network of CuO₂ 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

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

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 YBa₂Cu₃O₄ (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 t²/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)₂CuO₄ (LSCO) is estimated at 0.029 meV/(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 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].

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,

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 La_2–x–yNd_ySr_xCuO₄ [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 La_2–xSr_xCuO₄ (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 La_2–xSr_xCuO₄ 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

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,

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

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

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 Bi₂Sr₂CaCu₂O_8+d, [146] Bi₂Sr_2–xLa_xCuO_6+d, [147] and Ca_2–xNa_xCuO₂Cl₂ (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 Sr₂RuO₄ 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

(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 CuO₂ 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

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.

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

(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.

REFERENCES