- Email: [email protected]
Systematic improvement of wavefunctions in the variational Monte Carlo method for the t–J model
Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
www.elsevier.com/locate/jpcs
Systematic improvement of wavefunctions in the variational Monte Carlo method for the t±J model Masanori Kohno a,*, Masatoshi Imada b a
Frontier Science Institute, Mitsubishi Research Institute, Inc., 3-6 Otemachi 2-chome, Chiyoda-ku, Tokyo 100-8141, Japan b ISSP, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan
Abstract We present an algorithm to systematically improve wavefunctions in the variational Monte Carlo method for the t±J model. The ground-state wavefunction is approximated by a linear combination of single-particle states with Gutzwiller projection. The single-particle states are successively generated according to the numerical path-integral renormalization group (PIRG) procedure. In the initial step, the present method reduces to the usual variational Monte Carlo method. As the number of singleparticle states increases, the wavefunction converges to the ground-state wavefunction. Applying this method to the twodimensional t±J model, we have con®rmed that the wavefunction, which is composed of a small number of single-particle states, gives a good estimate of the ground-state energy. This algorithm does not suffer from the negative sign problem and can be applied to the models with strong correlations. q 2002 Published by Elsevier Science Ltd. Keywords: A. Superconductors; C. Ab initio calculations; D. Superconductivity
1. Introduction Strong electronic correlations cause a rich variety of phenomena for which, in some cases, our intuitive and naive understanding is misleading. High-Tc superconductivity is one of the phenomena whose mechanism has not been clari®ed yet. The two-dimensional t±J model is considered as a minimal model of high-Tc cuprate superconductors [1]. However, for this model, there are some unsettled issues such as whether the model really exhibits superconductivity in the realistic parameter regime. In principle, if numerical simulations for large systems could be performed, the properties of this model would be clari®ed. However, accuracy of numerical results extrapolated in the thermodynamic limit is not enough to show de®nite solutions for subtle issues such as the Mott transition and phaseseparation boundary [2±11]. For correlated electron systems, quantum Monte Carlo methods are powerful in some cases. For example, in the two-dimensional Hubbard model, various features of the Mott transition have been clari®ed by the auxiliary ®eld * Corresponding author. Tel.: 181-3-3277-0898; fax: 181-33277-3482. E-mail address: [email protected] (M. Kohno).
quantum Monte Carlo method [12±14]. However, in the strong-correlation regime (large-U regime), the negative sign problem becomes so severe that it is hard to get reliable results. Recently, an algorithm of the Hubbard model has been proposed by Imada and Kashima [15]. The algorithm, which they call the path-integral renormalization group (PIRG) method, does not suffer from the negative sign problem and allows us to investigate the larger-U regime than the quantum Monte Carlo method [16,17]. However, simple application of the PIRG method to the twodimensional t±J model does not work well near half-®lling, as we will explain later. In this paper, we propose a prescription for this issue by incorporating the variational Monte Carlo method with the PIRG method. The variational Monte Carlo (VMC) method is a numerical approach to investigate the ground-state properties of large-sized systems [18]. In this method, a variational wavefunction is assumed. The variational parameter is tuned to lower the energy expectation value. We can obtain approximate properties of the model within the assumed variational wavefunction. However, it is not sure how well the results obtained by VMC re¯ect the true ground-state properties of the model, because the assumed variational wavefunction is not necessarily close to the ground-state wavefunction. If the variational wavefunction is appropriately chosen,
0022-3697/02/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0022-369 7(02)00047-1
1564
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
the properties of the wavefunction are close to those of the ground state. In this paper, we present an algorithm to systematically improve the variational wavefunction by combining the PIRG procedure.
2. VMC and PIRG The algorithm we will present in this paper is based on the VMC method and the PIRG method. First, we shortly review these two methods. In the VMC method for Hubbard-type models, the Gutzwiller wavefunction [19±21] is used as the simplest variational wavefunction [22±25]. The Gutzwiller wavefunction is a single-particle state with Gutzwiller projection, which is de®ned as follows uGWl ;
Y i
1 2 1 2 gni" ni# uwl;
where variational parameter g controls the mean occupation number of electrons per site. If g is unity, the Gutzwiller wavefunction reduces to the simple single-particle state uwl: If g is zero, doubly occupied sites are not allowed. For the t± J model, g should be zero, because double occupation is forbidden. Hereafter, we write the Gutzwiller projection for the t±J model as Pd. In the VMC method, expectation values of physical operators are evaluated by Monte Carlo sampling. The Gutzwiller projection is explicitly taken into account through sampling of real-space con®gurations. The starting point of our algorithm is the VMC method. We extend the VMC method by generalizing a single-particle state uwl to a linear combination of single-particle states. If a suf®ciently large number of single-particle states are retained, the wavefunction reduces to that of the ground state. In practice, we apply the PIRG procedure in order to generate single-particle states systematically and to make the wavefunction converge to that of the ground state with a small number of single-particle states. In the PIRG method [15,16], a wavefunction is expressed as a linear combination of single-particle states as follows uCl ;
m X i1
decomposed by auxiliary ®eld si as follows [26] exp2DtHU
1YX UDt ni" 1 ni# ; exp 2asi ni" 2 ni# 2 2 2 i si
p where a tanh21 tanh UDt=4: To make ef®cient convergence with small number of m, the non-orthogonality of the base uwi l i 1±m is important. Through the update of the single-particle states by the above process, the wavefunction is optimized with number of single-particle states m ®xed. After the optimization of the wavefunction with m ®xed, we increase the number of single-particle states. In this way, the wavefunction becomes closer to that of the ground state. The ef®ciency of this method for the Hubbard model is demonstrated in Refs. [15±17]. 3. Combined algorithm of VMC and PIRG for the t±J model It is possible to apply the PIRG method to the t±J model straightforwardly by taking U ! 1 limit and adding J-term. However, the number of single-particle states should be very large near half-®lling in order to get accurate results, because a simple application of the PIRG method leads to generating an orthogonal and inef®cient base in the U ! 1 limit. This is easily understood, if one considers the case of half-®lling. The decomposed Hubbard U-term in the limit of U ! 1 corresponds to a very strong local magnetic ®eld that ®xes the spin orientation of each electron. At half®lling, this strong magnetic ®eld completely determines a real-space con®guration, which is orthogonal to those with any different magnetic ®elds. In fact, it is known that the Gutzwiller wavefunction gives a lower energy than a single real-space con®guration for the t±J model [24,25]. Hence, we use the Gutzwiller wavefunction instead of the simple single-particle state in the PIRG method. The wavefunction in our algorithm is expressed as follows u Cl ; P d
c i uwi l;
where uwi l is a single-particle state and ci is a coef®cient which is determined to lower the expectation value E of the Hamiltonian H E ; kCuHuCl=kCuCl: The singleparticle states are generated as follows uw i l exp2DtHuw0 l; where uw0 l is a single-particle state. The multiplication procedure follows that of the auxiliary ®eld quantum Monte Carlo method. For example, the Hubbard U-term is
m X i1
ci uwi l;
where uwi l is a single-particle state. The single-particle states are updated in the same way as in the PIRG method. The expectation value of physical operators is evaluated by the VMC method, where Gutzwiller projection is explicitly taken into account through sampling of real-space con®gurations. It should be noted that, for the PIRG method, the expectation value is obtained from Green's functions without insertion of the real-space representation. We summarize our algorithm in the following. (1) Initial state. By the usual variational Monte Carlo method, an optimized wavefunction is constructed. This wavefunction is used as the initial state uw1 l of our
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
Fig. 1. Energy per site as a function of the number of single-particle states for the two-dimensional t±J model at J=t 0:5 for a 10-site system with eight electrons.
algorithm. We may use this wavefunction as a guiding function that produces real-space con®gurations. (2) Production of single-particle states. Single-particle state uwm l is generated from uwi l by the operation of exp2DtHuwi l; where uwi l is a single-particle state obtained before i 1±m 2 1: (3) Evaluation of energy. Matrix elements kw~ i uHuw~ j l and kw~ i uw~ j l i; j 1±m are calculated by the Monte Carlo sampling through insertion of real-space representation X kw~ i uHuw~ j l kw~ i ux1 lkx1 uHux2 lkx2 uw~ j l; x1 ;x2
where uw~ j l ; Pd uwj l: The Gutzwiller projection Pd can be explicitly taken into account in the real-space representation. Using these matrix elements, we obtain a set of coef®cients cis which Pgives the lowest energy E kCuHuCl= kCuCl for uCl i c i uw~ i l: (4) Optimization of the wavefunction. Fixing the number of single-particle states m, we replace uw~ l l with a state uw~ 0l l which is obtained by the operation of exp2DtHuwl l l 1±m: If the energy becomes lower than that of uwl l; uw~ 0l l is accepted. Otherwise, uwl l is kept. Updating single-particle states uwl l l 1±m according to the above procedure several times, we obtain the optimal set of single-particle states uwl ls and coef®cient cls with the number of singleparticle states m ®xed. Repeating the processes (2)±(4), we obtain the optimized P wavefunction uCl i c i uw~ i l:
1565
Fig. 2. Energy per site as a function of the energy variance for the two-dimensional t±J model at J=t 0:5 for a 10-site system with eight electrons.
time for determinants is proportional to Ns3 and the number of auxiliary ®elds is proportional to Ns. In order to add m states and to sweep single-particle states, we need the computation time proportional to m 2. 5. Numerical results We have carried out calculation for the two-dimensional t±J model applying the above algorithm. The t±J model is de®ned as follows H Ht 1 HJ ; Ht 2t HJ J
X 1 c~ is c~js 1 H:c: ;
ki;jls
X ki;jl
Si ´Sj 2
1 ni nj ; 4
where c~is is the annihilation operator at site i with spin s with the constraint of no double occupancy. Operators S i and ni represent the spin and number operators at site i, respectively. The summation is taken over nearest neighbor sites. Hereafter, we set t 1 as the energy unit. The numerical result for the two-dimensional t±J model is shown in Fig. 1. As the number of single-particle states increases, the energy gets closer to that of the ground state. Although the size of the Hilbert space is 3150, only a hundred states are enough to show convergence. In order to extrapolate the energy smoothly, we have calculated energy variance D E, which is de®ned as
4. Computation time
DE ; kCuH 2 uCl 2 kCuHuCl 2 =kCuHuCl 2 :
The computation time is roughly proportional to Ms £ Ns4 £ m2 ; where Ms is the number of Monte Carlo samples, Ns is the number of sites and m is the number of singleparticle states. In each evaluation process of the energy, we have to compute Ms determinants. The computation
In general, the energy difference from the ground state dE ; E 2 Eg asymptotically scales as dE / DE : We plot energy difference d E as a function of energy variance D E in Fig. 2. Energy difference d E is well proportional to D E even with small m.
1566
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
6. Summary
References
We have presented a new algorithm which can be applied to strongly correlated electron systems such as the t±J model. This method is based on the PIRG method combined with the VMC method. By introducing the correlation effect through the Gutzwiller projection, we have developed the PIRG method to be applicable to the case of strong on-site Coulomb repulsion near half-®lling. This algorithm does not suffer from the negative sign problem. The computation time is a polynomial order of the system size and the number of single-particle states. Numerical results for the twodimensional t±J model obtained by this method show that a small number of single-particle states are enough to give a good estimate of the ground-state energy. The extrapolation from a small number of states shows convergent behavior to the ground state. Developing ef®cient algorithm for t±J type models is a long-sort and dif®cult issue. Our method offers a prescription for a powerful numerical method. Although the system size we have investigated is small, the present method would be applicable to large-sized systems without severe problems concerning the memory size of computers and the computation time. This method would also be applicable not only to the t±J model but also to other strongly correlated electron models such as the Hubbard model with very large U. Application to large-sized systems and ef®ciency in comparison to other methods need further study.
[1] F.C. Zhang, T.M. Rice, Phys. Rev. B 37 (1988) 3759. [2] V.J. Emery, S.A. Kivelson, H.Q. Lin, Phys. Rev. Lett. 64 (1990) 475. [3] W.O. Putikka, M.U. Luchini, T.M. Rice, Phys. Rev. Lett. 68 (1992) 538. [4] E. Daggoto, A. Moreo, F. Ortolani, D. Poslblone, J. Breva, Phys. Rev. B 45 (1992) 10741. [5] P. Prelovsek, X. Zotos, Phys. Rev. B 47 (1993) 5984. [6] D. Poilblanc, Phys. Rev. B 52 (1995) 9201. [7] M. Kohno, Phys. Rev. B 55 (1997) 1435. [8] C.S. Hellberg, E. Manousakis, Phys. Rev. Lett. 78 (1997) 4609. [9] C.S. Hellberg, E. Manousakis, Phys. Rev. B 61 (2000) 11787. [10] W.O. Putikka, M.U. Luchini, Phys. Rev. B 62 (2000) 1684. [11] C.T. Shin, Y.C. Chen, T.K. Lee, cond-mat/0104067. [12] M. Imada, Y. Hatsugai, J. Phys. Soc. Jpn 58 (1989) 3752. [13] N. Furukawa, M. Imada, J. Phys. Soc. Jpn 61 (1992) 3331. [14] M. Imada, A. Fujimori, Y. Tokura, Rev. Mod. Phys. 70 (1998) 1039. [15] M. Imada, T. Kashima, J. Phys. Soc. Jpn 69 (2000) 2723. [16] T. Kashima, M. Imada, J. Phys. Soc. Jpn 70 (2001) 2287. [17] T. Kashima, M. Imada, cond-mat/0104348. [18] D. Ceperley, G.V. Chester, K.H. Kalos, Phys. Rev. B 16 (1977) 3081. [19] C. Gutzwiller, Phys. Rev. Lett. 10 (1963) 159. [20] C. Gutzwiller, Phys. Rev. 134 (1964) A923. [21] C. Gutzwiller, Phys. Rev. 137 (1965) A1726. [22] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 1490. [23] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 3582. [24] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 3570. [25] C. Gros, R. Joynt, T.M. Rice, Phys. Rev. B 36 (1987) 381. [26] J.E. Hirsch, Phys. Rev. B 28 (1983) 4059.
Acknowledgements The authors would like to thank T. Kashima for useful discussions. This work is supported by `Research for the Future' program from the Japan Society for the Promotion of Science under the grant number JSPS-RFTF97P01103.
www.elsevier.com/locate/jpcs
Systematic improvement of wavefunctions in the variational Monte Carlo method for the t±J model Masanori Kohno a,*, Masatoshi Imada b a
Frontier Science Institute, Mitsubishi Research Institute, Inc., 3-6 Otemachi 2-chome, Chiyoda-ku, Tokyo 100-8141, Japan b ISSP, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan
Abstract We present an algorithm to systematically improve wavefunctions in the variational Monte Carlo method for the t±J model. The ground-state wavefunction is approximated by a linear combination of single-particle states with Gutzwiller projection. The single-particle states are successively generated according to the numerical path-integral renormalization group (PIRG) procedure. In the initial step, the present method reduces to the usual variational Monte Carlo method. As the number of singleparticle states increases, the wavefunction converges to the ground-state wavefunction. Applying this method to the twodimensional t±J model, we have con®rmed that the wavefunction, which is composed of a small number of single-particle states, gives a good estimate of the ground-state energy. This algorithm does not suffer from the negative sign problem and can be applied to the models with strong correlations. q 2002 Published by Elsevier Science Ltd. Keywords: A. Superconductors; C. Ab initio calculations; D. Superconductivity
1. Introduction Strong electronic correlations cause a rich variety of phenomena for which, in some cases, our intuitive and naive understanding is misleading. High-Tc superconductivity is one of the phenomena whose mechanism has not been clari®ed yet. The two-dimensional t±J model is considered as a minimal model of high-Tc cuprate superconductors [1]. However, for this model, there are some unsettled issues such as whether the model really exhibits superconductivity in the realistic parameter regime. In principle, if numerical simulations for large systems could be performed, the properties of this model would be clari®ed. However, accuracy of numerical results extrapolated in the thermodynamic limit is not enough to show de®nite solutions for subtle issues such as the Mott transition and phaseseparation boundary [2±11]. For correlated electron systems, quantum Monte Carlo methods are powerful in some cases. For example, in the two-dimensional Hubbard model, various features of the Mott transition have been clari®ed by the auxiliary ®eld * Corresponding author. Tel.: 181-3-3277-0898; fax: 181-33277-3482. E-mail address: [email protected] (M. Kohno).
quantum Monte Carlo method [12±14]. However, in the strong-correlation regime (large-U regime), the negative sign problem becomes so severe that it is hard to get reliable results. Recently, an algorithm of the Hubbard model has been proposed by Imada and Kashima [15]. The algorithm, which they call the path-integral renormalization group (PIRG) method, does not suffer from the negative sign problem and allows us to investigate the larger-U regime than the quantum Monte Carlo method [16,17]. However, simple application of the PIRG method to the twodimensional t±J model does not work well near half-®lling, as we will explain later. In this paper, we propose a prescription for this issue by incorporating the variational Monte Carlo method with the PIRG method. The variational Monte Carlo (VMC) method is a numerical approach to investigate the ground-state properties of large-sized systems [18]. In this method, a variational wavefunction is assumed. The variational parameter is tuned to lower the energy expectation value. We can obtain approximate properties of the model within the assumed variational wavefunction. However, it is not sure how well the results obtained by VMC re¯ect the true ground-state properties of the model, because the assumed variational wavefunction is not necessarily close to the ground-state wavefunction. If the variational wavefunction is appropriately chosen,
0022-3697/02/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0022-369 7(02)00047-1
1564
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
the properties of the wavefunction are close to those of the ground state. In this paper, we present an algorithm to systematically improve the variational wavefunction by combining the PIRG procedure.
2. VMC and PIRG The algorithm we will present in this paper is based on the VMC method and the PIRG method. First, we shortly review these two methods. In the VMC method for Hubbard-type models, the Gutzwiller wavefunction [19±21] is used as the simplest variational wavefunction [22±25]. The Gutzwiller wavefunction is a single-particle state with Gutzwiller projection, which is de®ned as follows uGWl ;
Y i
1 2 1 2 gni" ni# uwl;
where variational parameter g controls the mean occupation number of electrons per site. If g is unity, the Gutzwiller wavefunction reduces to the simple single-particle state uwl: If g is zero, doubly occupied sites are not allowed. For the t± J model, g should be zero, because double occupation is forbidden. Hereafter, we write the Gutzwiller projection for the t±J model as Pd. In the VMC method, expectation values of physical operators are evaluated by Monte Carlo sampling. The Gutzwiller projection is explicitly taken into account through sampling of real-space con®gurations. The starting point of our algorithm is the VMC method. We extend the VMC method by generalizing a single-particle state uwl to a linear combination of single-particle states. If a suf®ciently large number of single-particle states are retained, the wavefunction reduces to that of the ground state. In practice, we apply the PIRG procedure in order to generate single-particle states systematically and to make the wavefunction converge to that of the ground state with a small number of single-particle states. In the PIRG method [15,16], a wavefunction is expressed as a linear combination of single-particle states as follows uCl ;
m X i1
decomposed by auxiliary ®eld si as follows [26] exp2DtHU
1YX UDt ni" 1 ni# ; exp 2asi ni" 2 ni# 2 2 2 i si
p where a tanh21 tanh UDt=4: To make ef®cient convergence with small number of m, the non-orthogonality of the base uwi l i 1±m is important. Through the update of the single-particle states by the above process, the wavefunction is optimized with number of single-particle states m ®xed. After the optimization of the wavefunction with m ®xed, we increase the number of single-particle states. In this way, the wavefunction becomes closer to that of the ground state. The ef®ciency of this method for the Hubbard model is demonstrated in Refs. [15±17]. 3. Combined algorithm of VMC and PIRG for the t±J model It is possible to apply the PIRG method to the t±J model straightforwardly by taking U ! 1 limit and adding J-term. However, the number of single-particle states should be very large near half-®lling in order to get accurate results, because a simple application of the PIRG method leads to generating an orthogonal and inef®cient base in the U ! 1 limit. This is easily understood, if one considers the case of half-®lling. The decomposed Hubbard U-term in the limit of U ! 1 corresponds to a very strong local magnetic ®eld that ®xes the spin orientation of each electron. At half®lling, this strong magnetic ®eld completely determines a real-space con®guration, which is orthogonal to those with any different magnetic ®elds. In fact, it is known that the Gutzwiller wavefunction gives a lower energy than a single real-space con®guration for the t±J model [24,25]. Hence, we use the Gutzwiller wavefunction instead of the simple single-particle state in the PIRG method. The wavefunction in our algorithm is expressed as follows u Cl ; P d
c i uwi l;
where uwi l is a single-particle state and ci is a coef®cient which is determined to lower the expectation value E of the Hamiltonian H E ; kCuHuCl=kCuCl: The singleparticle states are generated as follows uw i l exp2DtHuw0 l; where uw0 l is a single-particle state. The multiplication procedure follows that of the auxiliary ®eld quantum Monte Carlo method. For example, the Hubbard U-term is
m X i1
ci uwi l;
where uwi l is a single-particle state. The single-particle states are updated in the same way as in the PIRG method. The expectation value of physical operators is evaluated by the VMC method, where Gutzwiller projection is explicitly taken into account through sampling of real-space con®gurations. It should be noted that, for the PIRG method, the expectation value is obtained from Green's functions without insertion of the real-space representation. We summarize our algorithm in the following. (1) Initial state. By the usual variational Monte Carlo method, an optimized wavefunction is constructed. This wavefunction is used as the initial state uw1 l of our
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
Fig. 1. Energy per site as a function of the number of single-particle states for the two-dimensional t±J model at J=t 0:5 for a 10-site system with eight electrons.
algorithm. We may use this wavefunction as a guiding function that produces real-space con®gurations. (2) Production of single-particle states. Single-particle state uwm l is generated from uwi l by the operation of exp2DtHuwi l; where uwi l is a single-particle state obtained before i 1±m 2 1: (3) Evaluation of energy. Matrix elements kw~ i uHuw~ j l and kw~ i uw~ j l i; j 1±m are calculated by the Monte Carlo sampling through insertion of real-space representation X kw~ i uHuw~ j l kw~ i ux1 lkx1 uHux2 lkx2 uw~ j l; x1 ;x2
where uw~ j l ; Pd uwj l: The Gutzwiller projection Pd can be explicitly taken into account in the real-space representation. Using these matrix elements, we obtain a set of coef®cients cis which Pgives the lowest energy E kCuHuCl= kCuCl for uCl i c i uw~ i l: (4) Optimization of the wavefunction. Fixing the number of single-particle states m, we replace uw~ l l with a state uw~ 0l l which is obtained by the operation of exp2DtHuwl l l 1±m: If the energy becomes lower than that of uwl l; uw~ 0l l is accepted. Otherwise, uwl l is kept. Updating single-particle states uwl l l 1±m according to the above procedure several times, we obtain the optimal set of single-particle states uwl ls and coef®cient cls with the number of singleparticle states m ®xed. Repeating the processes (2)±(4), we obtain the optimized P wavefunction uCl i c i uw~ i l:
1565
Fig. 2. Energy per site as a function of the energy variance for the two-dimensional t±J model at J=t 0:5 for a 10-site system with eight electrons.
time for determinants is proportional to Ns3 and the number of auxiliary ®elds is proportional to Ns. In order to add m states and to sweep single-particle states, we need the computation time proportional to m 2. 5. Numerical results We have carried out calculation for the two-dimensional t±J model applying the above algorithm. The t±J model is de®ned as follows H Ht 1 HJ ; Ht 2t HJ J
X 1 c~ is c~js 1 H:c: ;
ki;jls
X ki;jl
Si ´Sj 2
1 ni nj ; 4
where c~is is the annihilation operator at site i with spin s with the constraint of no double occupancy. Operators S i and ni represent the spin and number operators at site i, respectively. The summation is taken over nearest neighbor sites. Hereafter, we set t 1 as the energy unit. The numerical result for the two-dimensional t±J model is shown in Fig. 1. As the number of single-particle states increases, the energy gets closer to that of the ground state. Although the size of the Hilbert space is 3150, only a hundred states are enough to show convergence. In order to extrapolate the energy smoothly, we have calculated energy variance D E, which is de®ned as
4. Computation time
DE ; kCuH 2 uCl 2 kCuHuCl 2 =kCuHuCl 2 :
The computation time is roughly proportional to Ms £ Ns4 £ m2 ; where Ms is the number of Monte Carlo samples, Ns is the number of sites and m is the number of singleparticle states. In each evaluation process of the energy, we have to compute Ms determinants. The computation
In general, the energy difference from the ground state dE ; E 2 Eg asymptotically scales as dE / DE : We plot energy difference d E as a function of energy variance D E in Fig. 2. Energy difference d E is well proportional to D E even with small m.
1566
M. Kohno, M. Imada / Journal of Physics and Chemistry of Solids 63 (2002) 1563±1566
6. Summary
References
We have presented a new algorithm which can be applied to strongly correlated electron systems such as the t±J model. This method is based on the PIRG method combined with the VMC method. By introducing the correlation effect through the Gutzwiller projection, we have developed the PIRG method to be applicable to the case of strong on-site Coulomb repulsion near half-®lling. This algorithm does not suffer from the negative sign problem. The computation time is a polynomial order of the system size and the number of single-particle states. Numerical results for the twodimensional t±J model obtained by this method show that a small number of single-particle states are enough to give a good estimate of the ground-state energy. The extrapolation from a small number of states shows convergent behavior to the ground state. Developing ef®cient algorithm for t±J type models is a long-sort and dif®cult issue. Our method offers a prescription for a powerful numerical method. Although the system size we have investigated is small, the present method would be applicable to large-sized systems without severe problems concerning the memory size of computers and the computation time. This method would also be applicable not only to the t±J model but also to other strongly correlated electron models such as the Hubbard model with very large U. Application to large-sized systems and ef®ciency in comparison to other methods need further study.
[1] F.C. Zhang, T.M. Rice, Phys. Rev. B 37 (1988) 3759. [2] V.J. Emery, S.A. Kivelson, H.Q. Lin, Phys. Rev. Lett. 64 (1990) 475. [3] W.O. Putikka, M.U. Luchini, T.M. Rice, Phys. Rev. Lett. 68 (1992) 538. [4] E. Daggoto, A. Moreo, F. Ortolani, D. Poslblone, J. Breva, Phys. Rev. B 45 (1992) 10741. [5] P. Prelovsek, X. Zotos, Phys. Rev. B 47 (1993) 5984. [6] D. Poilblanc, Phys. Rev. B 52 (1995) 9201. [7] M. Kohno, Phys. Rev. B 55 (1997) 1435. [8] C.S. Hellberg, E. Manousakis, Phys. Rev. Lett. 78 (1997) 4609. [9] C.S. Hellberg, E. Manousakis, Phys. Rev. B 61 (2000) 11787. [10] W.O. Putikka, M.U. Luchini, Phys. Rev. B 62 (2000) 1684. [11] C.T. Shin, Y.C. Chen, T.K. Lee, cond-mat/0104067. [12] M. Imada, Y. Hatsugai, J. Phys. Soc. Jpn 58 (1989) 3752. [13] N. Furukawa, M. Imada, J. Phys. Soc. Jpn 61 (1992) 3331. [14] M. Imada, A. Fujimori, Y. Tokura, Rev. Mod. Phys. 70 (1998) 1039. [15] M. Imada, T. Kashima, J. Phys. Soc. Jpn 69 (2000) 2723. [16] T. Kashima, M. Imada, J. Phys. Soc. Jpn 70 (2001) 2287. [17] T. Kashima, M. Imada, cond-mat/0104348. [18] D. Ceperley, G.V. Chester, K.H. Kalos, Phys. Rev. B 16 (1977) 3081. [19] C. Gutzwiller, Phys. Rev. Lett. 10 (1963) 159. [20] C. Gutzwiller, Phys. Rev. 134 (1964) A923. [21] C. Gutzwiller, Phys. Rev. 137 (1965) A1726. [22] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 1490. [23] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 3582. [24] H. Yokoyama, H. Shiba, J. Phys. Soc. Jpn 56 (1987) 3570. [25] C. Gros, R. Joynt, T.M. Rice, Phys. Rev. B 36 (1987) 381. [26] J.E. Hirsch, Phys. Rev. B 28 (1983) 4059.
Acknowledgements The authors would like to thank T. Kashima for useful discussions. This work is supported by `Research for the Future' program from the Japan Society for the Promotion of Science under the grant number JSPS-RFTF97P01103.