Crossover from weakly to strongly correlated regions in the two-dimensional Hubbard model -- Off-diagonal wave function Monte Carlo studies of Hubbard model II --

The ground state of the two-dimensional (2D) Hubbard model is investigated by adopting improved wave functions that take into account intersite electron correlation beyond the Gutzwiller ansatz. The ground-state energy is lowered considerably, giving the best estimate of the ground-state energy for the 2D Hubbard model. There is a crossover from weakly to strongly correlated regions as the on-site Coulomb interaction $U$ increases. The antiferromagnetic correlation induced by $U$ is reduced for hole doping when $U$ is large, being greater than the bandwidth, thus increasing the kinetic energy gain. The spin and charge fluctuations are induced in the strongly correlated region. These antiferromagnetic and kinetic charge fluctuations induce electron pairings, which would result in high-temperature superconductivity.


I. INTRODUCTION
The mechanism of high-temperature superconductivity has been studied intensively since the discovery of cuprate high-temperature superconductors [1,2]. The correlation between electrons plays an important role because parent compounds without carriers are insulators. It is primarily important to clarify the phase diagram of electronic states in the CuO 2 plane contained commonly in cuprate high-temperature superconductors. It is obvious that an interaction with a large energy scale is necessary and responsible for the realization of hightemperature superconductivity. The Coulomb interaction obviously has a large characteristic energy scale and can possibly be a candidate of the interaction that induces high-temperature superconductivity.
A system described by the Hubbard model is a typical strongly correlated system. The 2D Hubbard model has been investigated intensively for several decades [16,[31][32][33][34][35][36][37]. The Hubbard model was first introduced to describe a metal-insulator transition [16]. Since the discovery of cuprate high-temperature superconductors, many researchers have tried to explain the occurrence of superconductivity in cuprates in terms of the 2D Hubbard model. The results of quantum Monte Carlo methods, which are considered to be exact unbiased methods, do not support the existence of high-temperature supercon-ductivity in this model [19][20][21]. In our opinion, this is because the Coulomb interaction U is not large in quantum Monte Carlo calculations, where the range of accessible U is very restricted because of the nature of the method. On the basis of the variational Monte Carlo method, a finite superconducting gap with d-wave symmetry is hardly obtained for U/t < 6 in the 2D Hubbard model [38][39][40][41][42][43].
It is necessary to improve the wave function in the variational Monte Carlo method since the Gutzwiller function only accounts for the on-site correlation. The Gutzwiller function is a starting function that should be improved to take account of correlation effects. We discuss several methods of improving the wave function, and show that an exp(−λK)-P G -type function [44][45][46] can be the best function with the lowest variational energy. This wave function and its generalization were examined thoroughly on the basis of the variational Monte Carlo method in Ref. 45, which is referred to as paper I in this paper. The variational energy is lowered considerably and can be close to the exact value.
We discuss the properties of antiferromagnetism and superconductivity in a strongly correlated region by using improved wave functions. It may be believed that the antiferromagnetic (AF) correlation is enhanced as U increases and that this will hold for large U where U is far greater than the bandwidth. This is, however, not true when holes are doped, as will be shown on the basis of improved wave functions. The AF correlation is suppressed when U is extremely large and larger than the bandwidth. A superconducting correlation is developed in this region as the AF correlation is suppressed with the increase in U . The development of a superconducting correlation is understood to be induced by spin fluctuation, which is induced by the kinetic charge fluctuation that conquers antiferromagnetism. This spin fluctuation in a strongly correlated region must be distinguished from that in a weakly correlated region. The latter is the conventional spin fluctuation that has been discussed extensively in the literature [47,48].
Thus, there is a crossover between weakly correlated and strongly correlated regions as the Coulomb repulsion U is varied. A high-temperature superconductivity will be realized in the strongly correlated region.
In the next section we discuss improved wave functions as below in correlated electron systems. In Sect. 3, we examine the stability of the AF and pairing states, on the basis of our wave functions. We give a summary in the last section.

II. HAMILTONIAN AND WAVE FUNCTIONS
A. Two-dimensional Hubbard model The single-band Hubbard model is given by where t ij are transfer integrals and U is the on-site Coulomb energy. n iσ is the number operator given by n iσ = c † iσ c iσ . The transfer integral t ij is nonzero, t ij = −t for the nearest-neighbor pair ij and t ij = −t ′ for the next-nearest-neighbor pair ij . Otherwise, t ij vanishes. We denote the number of sites as N and the number of electrons as N e . The energy unit is given by t. The second term represents the on-site Coulomb repulsive interaction between electrons with opposite spins. This simple term will illustrate profound phenomena in strongly correlated electron systems.

B. Improved wave functions
The wave function should include correlation between electrons. A first step to include the electron correlation is the well-known Gutzwiller wave function, which reads as ψ G = P G ψ 0 , where P G is the Gutzwiller operator defined by with the variational parameter g in the range of 0 ≤ g ≤ 1. P G controls the double occupancy to take account of electron correlation. ψ 0 is a trial one-particle state that is taken to be the Fermi sea ψ F S , the AF state (spindensity wave state) ψ AF , or the BCS state ψ BCS . The AF one-particle state ψ AF is given by the eigenstate of the AF trial Hamiltonian: where r i ≡ (x i , y i ) are the coordinates of site i. ∆ AF indicates the AF order parameter. In general, the transfer integrals t ij in the trial Hamiltonian H trial can be treated as variational parameters so that the ground-state energy is lowered. In this paper, however, we do not use this procedure for simplicity, and thus t ij in H trial are the same as those in the original Hamiltonian. One must improve the Gutzwiller wave function to lower the ground-state energy because only the on-site correlation is considered in the Gutzwiller ansatz. It has been proposed that the wave function can be improved by taking account of the nearest-neighbor doublon-holon correlation [49][50][51][52][53]. The wave function with the doublonholon correlation is given by Here, d j is the operator for the doubly occupied site given by d j = n j↑ n j↓ , and e j is the empty site operator given as e j = (1−n j↑ )(1−n j↓ ). η is a variational parameter in the range of 0 ≤ η ≤ 1. We put η = 1 in the non-interacting case. A Jastrow factor is defined as where n i = n i↑ + n i↓ and {h ij } are variational parameters. We can take into account intersite correlations by multiplying by P J to obtain wave functions such as A more efficient way to improve the wave function is to take account of the intersite correlation by multiplying the Gutzwiller function by the kinetic operator. A wave function of this type is written as [44,45,54] where K is the kinetic term in the Hamiltonian: K = ijσ t ij c † iσ c jσ and λ is a variational parameter to be optimized to lower the energy. The transfer integrals t ij in K are not regarded as variational parameters in this paper. This wave function has been investigated by the variational Monte Carlo method [44,45,54,55] and a perturbative method [46,[56][57][58].
This wave function can be further improved by multiplying by the Gutzwiller operator and e −λK again: The wave function ψ (m) is called the level-m wave function in this paper. This type of wave functions was first examined by the variational Monte Carlo method in paper I [45]. The wave functions ψ (2) , ψ (3) , ψ (4) ,· · · contain intersite correlations of the site-off-diagonal type such as c † iσ c jσ , whereas the Gutzwiller and Jastrow operators control only site-diagonal correlations such as n iσ n jσ ′ .
For the Gutzwiller-projected BCS-type wave function, we take ψ BCS as a trial one-particle state ψ 0 : with coefficients u k and v k appearing only in the ra- where ∆ k is the kdependent gap function and ξ k is the band dispersion given as We multiply by P Ne to extract only the state with the fixed total electron number N e . For the wave functions ψ (2) , ψ (3) ,· · · , it is difficult to impose the sitediagonal-type constraint P Ne to fix the total electron number. To introduce the off-diagonal long-range order in the wave function with off-diagonal site correlation, we perform the electron-hole transformation for down-spin electrons [54]: The up-spin electrons are unaltered, and we use the notation c k = c k↑ . If we denote the vacuum for c and d particles as |0 , for which c k |0 = d k |0 = 0, the vacuum |0 reads as |0 = k d † k |0 . Then, the BCS wave function is written as In this wave function, c and d particles are mixed and the total electron number is controlled by the chemical potential µ. We adjust the chemical potential so that we have the expectation value N e for the total electron number.
In the above transformation, the Gutzwiller operator is mapped to and the Hamiltonian is transformed to where c j and d j are the Fourier transforms of c k and d k , respectively. The averaged total electron number is written as C. Ground-state energy of improved wave functions We calculate the ground state energy to check the validity of our wave functions. Recently, a variational wave function has been proposed by introducing a large number of variational parameters in the Gutzwiller-Jastrow factor and the one-particle function [59,60]. The energy evaluated using a wave function with many variational parameters is lowered considerably compared with that using the Gutzwiller function ψ G . Our wave function ψ λ also gives a good estimate of the energy, and P G ψ λ and e −γK P G ψ λ give the best ground-state energy. We show the variational energy calculated on a 4 × 4 lattice in Table I for a half-filled case and in Table II for hole doping. The many-parameter wave function in Ref. 59 is a good wave function with much lower variational energy. The variational energy of our wave functions is also lowered considerably compared with that of the Gutzwiller wave function and exhibits the best ground-state energy. We show in Table III the expectation values of the spin correlation function S(q) for q = (π, π/2) on a 4 × 4 lattice with N e = 12, t ′ /t = 0, and U/t = 10. The improved ψ (2) shows good agreement with the exact value.
It is seen that the energy is not significantly improved only by multiplying the Gutzwiller function by the doublon-holon correlation factor P d−h compared with that of the wave function with site-off-diagonal correlation. The trial wave function P d−h P G ψ 0 was used to develop the physics of the Mott transition [49] following the suggestion that the Mott transition is due to doublonholon binding [61]. We examined the Mott transition using the wave function e −λK P G ψ 0 [55] in a previous paper because the variational energy of this wave function is much lower than that of the doublon-holon wave function. The wave function ψ λ ≡ ψ (2) clearly exhibited the Mott transition as U increases.
The other way to improve wave functions is to regard the band parameters in ψ 0 as variational parameters, that is, to take account of the band renormalization. There have been works in this direction [51,53,62,63] We do not adopt this method in this paper to minimize the number of variational parameters. In Table IV, we show the ground-state energy to compare wave functions on a 10 × 10 lattice (which we mainly consider in this paper). The level-2 wave function ψ (2) gives a good estimate of the ground-state energy.

III. CROSSOVER FROM WEAKLY TO STRONGLY CORRELATED STATES
A. Ground-state energy We evaluate the ground-state energy of the 2D Hubbard model with hole doping using ψ λ , where the oneparticle state is taken as the Fermi sea, ψ 0 = ψ F S . The details of the methods of Monte Carlo calculations with ψ (m) are shown in paper I. We show the ground-state energy as a function of U in Fig. 1, where the calculations were performed on a 10 × 10 lattice with N e = 88 and t ′ = 0. Solid circles indicate the energy obtained by the Gutzwiller function and open circles denote that I: Variational energies for 4 × 4 lattice, Ne = 16, U/t = 5, and t ′ /t = 0. The boundary conditions are periodic in one direction and antiperiodic in the other direction. The nearest-neighbor (n.n.) Jastrow function in the third row means that we considered only the nearest-neighbor Jastrow correlation factor exp(− ij hij ninj ). The oneparticle state ψF S indicates the Fermi sea. The result in the fourth row is from Ref. 59   obtained using ψ λ . Squares represent the expectation values of the kinetic term (non-interacting part of the Hamiltonian) E K ≡ ijσ t ij c † iσ c jσ ; solid squares are for the Gutzwiller function and open squares are for ψ λ . The ground-state energy is lowered considerably by the optimization with the λ parameter. The increase in kinetic energy gain is appreciable, as shown in Fig. 1. In this evaluation, we performed 5 × 10 4 Monte Carlo steps with about 200 parallel processors.
The double occupancy i n i↑ n i↓ /N ≡ E U /U N and the expectation value of the Coulomb potential term E U /N ≡ U i n i↑ n i↓ /N are shown in Fig. 2 and Fig.  3, respectively. while the expectation value E G U obtained by the Gutzwiller function is larger than E λ U obtained by using ψ λ in the weakly correlated region, E λ U becomes greater than E G U for large U in the strongly correlated region. This suggests that the charge fluctuation is larger in ψ λ than in the Gutzwiller function. To see the role of the parameter λ in ψ λ , we show the energy as a function of λ in Fig. 4. The absolute value of E K increases as λ increases, which lowers the total ground-state energy. III: Spin correlation function S(π, π/2) for 4 × 4 lattice with Ne = 12 and t ′ /t = 0. The parameters are U/t = 4 and U/t = 10. We take the one-particle state ψ0 as the Fermi sea, ψ0 = ψF S . The last column shows S(π, π/2) for U/t = 10. The boundary conditions are the same as in Table  I.
Wave function U/t = 4 U/t = 10 (π, π/2) PGψF Wave function Thus, λ induces charge fluctuation to increase the kinetic energy gain. This will bring about a crossover from weakly to strongly correlated regions as U increases.

B. Crossover and antiferromagnetic correlation
In this subsection, we investigate the stability of the antiferromagnetic state as a function of U in the case of hole doping. The one-particle state is the AF state ψ 0 = ψ AF with the order parameter ∆ AF . The AF correlation is induced by the Coulomb repulsion U and increases as U increases. In contrast, it is suppressed if U becomes larger than the bandwidth in the strongly correlated region.
Let us consider the spin correlation function S(q) given as S(q) = 1 4N ij e iq·(ri−rj ) (n i↑ − n i↓ )(n j↑ − n j↓ ). (15) S(q) has a maximum at q = (π, π). The q = (π, π) component of S(q) is shown in Fig. 5 as a function of U for N e = 88 and t ′ = 0 on a 10 × 10 lattice, where the calculations are based on ψ λ . S(π, π) has a peak near U ≃ 10t, showing the reduction in spin correlation for large U . We show the antiferromagnetic (AF) order parameter ∆ AF as a function of U in Fig. 6. The AF order parameter is included in a trial wave function ψ 0 as a variational parameter. The optimized ∆ AF shows a peak structure around U ≃ 10t and vanishes in the strongly correlated region. This is a typical behavior of ∆ AF , and the AF energy gain ∆E AF also exhibits a similar behavior. The calculation was carried out by employing the wave function ψ λ on a 10 × 10 lattice. When U is small, ∆ AF increases with the increase in U and has a maximum at U co ≃ 10t, which is on the order of the bandwidth. When U is larger than U co , ∆ AF is decreased as U is increased. This indicates that AF correlation is suppressed for extremely large U and diminishes. In the region U > U co , there is competition between the AF correlation and the charge fluctuation; this means that we must have an AF energy gain or kinetic energy gain to lower the ground-state energy. ∆ AF is reduced gradually as U increases (U > U co ) since the energy gain ∆E AF is presumably proportional to the AF exchange coupling, J ∝ t 2 /U . The AF correlation should be suppressed to obtain a kinetic energy gain for large U . Thus, we have a weak AF correlation in the strongly correlated region with U ≥ U co . This indicates that there is a large AF fluctuation in this region, brought about by charge fluctuation, where the charge fluctuation is driven by the kinetic operator K in the exponential factor exp(−λK). The correct term for the charge fluctuation is the kinetic charge fluctuation.

C. Superconductivity in strongly correlated region
It has been found that there is a large spin fluctuation being driven by the kinetic charge fluctuation. We expect that a pairing interaction is induced by this kind of large spin and charge fluctuation.
We perform calculations with the superconducting order parameter in the wave function ψ λ . To do this, we use the electron-hole transformation for electrons with the down spin as shown in the previous section. The chemical potential is used to adjust the total electron number to be equal to N e . We assume the d-wave symmetry for electron pairing and introduce the order parameter ∆. A trial state is represented by a 2N × N matrix φ, and the correlated pairing state can be formulated following the method in Ref. 54. We use a real-space formulation by using the solution of the Bogoliubov-de Gennes equation [54,64]: for a trial Hamiltonian H ijσ and F ij , where (H ijσ ) and (F ij ) are N ×N matrices in the real space representation. The matrix (F ij ) includes the gap function ∆ as matrix elements between c and d particles. The initial matrix φ is given by φ iα = u α i and φ i+N,α = v α i for i = 1, · · · , N and α = 1, · · · , N with E α > 0.
The energy plotted against N e /2 is shown in Fig. 7 for U/t = 18. The figure indicates that the minimum of the energy is at ∆/t ∼ 0.06, giving E/N ≃ −0.5172 for N e = 88. In this case, the condensation energy per site is ∆E/N ≃ 0.0057t. The pairing state e −λK P G ψ BCS with finite ∆ is never stabilized for small U such as U < 6t. The optimized superconducting order parameter increases as U increases and has a maximum at some U and greater than U co . This is shown in Fig. 8, where the superconducting order parameter ∆ and the AF order parameter are shown as functions of U .
It is appropriate to call the region U > U co the strongly correlated region because the ground state at half-filling is insulating in this region [50,55]. It is difficult to find a clear sign of superconductivity in weakly correlated regions, where U ≤ 6t, using numerical methods such as a quantum Monte Carlo method [19,20,24]. This is also the case for the variational Monte Carlo method since the superconducting condensation energy vanishes or is very small for small U . Instead, in a strongly correlated region, we obtain a clear indication of superconductivity.  Fig. 1 In the small-U region, the double occupancy of the Gutzwiller function is greater than that of ψ λ , and this is reversed in the large-U region.

IV. SUMMARY
We have investigated the 2D Hubbard model by adopting improved wave functions that take into account intersite correlations beyond the Gutzwiller ansatz. Our wave function is an exp(−λK)-P G -type wave function, which is inspired by the wave function used in quantum Monte Carlo methods and can be improved systematically by multiplying by P G and e −λK . We can improve the variational energy considerably with this wave function. Our wave function can give the best estimate of the ground-state energy of the 2D Hubbard model.
The improved wave function gives results that are qualitatively different from those obtained by the simple Gutzwiller function. In particular, the picture of the stability of the antiferromagnetic state is crucially changed when we employ wave functions with intersite correlations. The antiferromagnetic state is stable when U is larger than a critical value at half-filling. Although one may think that this also holds for the hole-doped case, this is not true away from half-filling. To obtain the kinetic energy gain, we must eliminate the antiferromagnetic order when U is extremely large because the energy lowering by the antiferromagnetic order is small, being proportional to t 2 /U . Thus, the antiferromagnetic correlation will fade away as U increases to be greater than the critical U co .
The reduction in the antiferromagnetic correlation in the strongly correlated region suggests the existence of a large antiferromagnetic spin fluctuation. The charge fluctuation induced by the kinetic operator is apprecia-  Fig. 1 The expectation value of the Coulomb interaction obtained by the Gutzwiller function E G U (indicated by circles) is larger than the expectation value E λ U (squares) obtained using ψ λ for small U . For large U , we have E λ U > E G U . Ground-state energy (circles), kinetic energy (squares), and Coulomb energy EU /N (diamonds) as functions of λ on 10 × 10 lattice for U/t = 18. The number of electrons is Ne = 88 and we set t ′ = 0.0. The boundary conditions are the same as in Fig. 1 Open symbols indicate the values obtained with the variational parameter g fixed to the optimized value at which the total energy is minimized. For solid symbols, g is also optimized. As is clear in the figure, the kinetic energy gain increases as λ increases. This lowers the ground-state energy. ble and helps electrons to form pairs. We expect that these large spin and charge fluctuations induce an effective pairing interaction between electrons, which may be able to bring about high-temperature superconductivity.
There is also spin fluctuation in the weakly correlated region, where the antiferromagnetic correlation fades away as U is reduced. The next-nearest-neighbor transfer t ′ /t is expected to be important in determining the stability of the antiferromagnetic state. It has been argued that the groundstate property, especially the stability of ordered states, crucially depends on t ′ /t [53,59] and that the antiferromagnetic state will be stabilized due to t ′ /t. This will be an interesting future subject when studied on the basis of our improved wave functions.
A discussion on the relationship between the t-J model and the Hubbard model is given here. It is well known that the energy scale is given by J = 4t 2 /U in the large-U region, as shown by mapping the Hubbard model to a t-J-like model [65]. The U -dependence of the antiferromagnetic correlation can be understood from this mapping, namely, the AF correlation is suppressed in the large-U region such as U > 8t. This suggests that the physics in the strongly correlated region is similar to that of the t-J model. In the t-J model, the exchange interaction J induces superconductivity with d-wave pairing as well as AF ordering [66]. At half-filling, the ground state has an AF long-range order for J > 0. This AF ordering shows instability due to hole doping; this occurs even for small doping. In the case of the Hubbard model for large U , say U ∼ 20t, the AF state becomes unstable upon hole doping. This indicates a similarity between the t-J model and the Hubbard model in the strongly correlated region. There is, however, a difference between the two models. For example, the superconducting gap ∆ remains finite even for large U with small J, where the AF order is suppressed. This means that we cannot well understand ∆ as a function of U only by means of J. We consider that we must take account of the charge fluctuation to understand the superconducting state for large U . We must acknowledge that fluctuations are important, and in particular, the charge fluctuation, as well as spin fluctuation, plays a significant role in stabilizing the pairing state.
The crossover behavior from weakly to strongly cor-related regions or vice versa in the 2D Hubbard model may belong to a universality class. The well-known Kondo problem exhibits a crossover from the weak coupling region to the strong coupling region as the temperature T decreases [67]. The s-d model is mapped to an electron-gas model interacting with a logarithmic interaction, from which the scaling equation for the coupling constant J was derived by the renormalization group procedure [68,69]. The renormalization group equation for the s-d model agrees with that for the 2D sine-Gordon model [70][71][72]. The one-dimensional (1D) Hubbard model is mapped onto the 2D sine-Gordon model by a bosonization procedure [73,74]. Thus, these models (2D sine-Gordon model, s-d model, and 1D Hubbard model) are in the same universality class from the viewpoint of scaling theory. We expect that, concerning the crossover between weakly and strongly correlated regions, the 2D Hubbard model is also in this class given by the 2D sine-Gordon model.
In conclusion, the crossover occurs from a weakly correlated to a strongly correlated region, and hightemperature superconductivity will appear in the strongly correlated region.