- Email: [email protected]
On-site charge-density wave formation in one-dimensional chains investigated by Feynman path-integral quantum Monte Carlo calculations
Volume 202, number 1,2
CHEMICAL PHYSICS LETTERS
15 January 1993
On-site charge-density wave formation in one-dimensional chains investigated by Feynman path-integral quantum Monte Carlo calculations Michael C. Bbhm, Joachim Schulte InstitutftirPhysikalische Chemie, Physikafische Chemie IIL Technische Hochschule Darmsfadt, W-6100 Darmstadt, Germany
and Luis Utrera Institute de Ciencia Materiales, Consejo Superior de Investigationes Cieniijicas, 28006 Madrid, Spain
Received 9 June 1992;in final form 22 October 1992
The formation of on-site charge-density waves (CDWs) in the one-dimensional ( ID) Hubbard chain is studied by path-integral quantum Monte Carlo (QMC) simulations as a function of the band fdling, the ratio between on-site and nearest-neighbor intersite interaction as well as the net Coulomb interaction and kinetic hopping. The dynamics of the electrons of the 1D system are quantified by the charge fluctuations ( (An;)) and the probability of double occupancy P(2) averaged over all lattice sites i. Electronic conditions favouring CDWs are summarized. It is shown that one quarter, one half and three quarter tilled bands are optimum conditions for on-site CDWs.
1. Introduction The non-conventional material properties of quasi one-dimensional ( 1D) synthetic metals have adsorbed the interest of experimentalists and theorists for more than two decades now [ I ,2]. The rich variety of accessible low-temperature ground state configurations has been a challenging problem in the design of 1D systems with specific key quantities. Quasi 1D metals can undergo phase transitions either to charge/spin density wave (C/SDW), or to singlet/ triplet superconducting (S/TS) low-temperature modifications [3]. The symmetry of the electronic wavefunction in non-distorted C/SDW states is lower than the symmetry of the corresponding (1D) lattice. For recent one-determinantal band structure calculations of 1D materials which converged into on-site and intersite CDW configurations, see refs. 74,51. In finite electronic systems symmetry broken single determinantal mean-field solutions always indicate the shortcomings of an independent-particle
description. These so-called Hartree-Fock instabilities have been studied extensively in the past years [6,7]. Note that labels like one-electron model, independent-particle scheme, Hartree-Fock approximation, etc. are used to symbolize wavefunctions described by a single Slater determinant. In comparison with symmetry broken mean-field wavefunctions observed in finite molecules, the situation is less definite if one-electron band-structure calculations of 1D solids converge into CDW type solutions. On one hand they may indicate a wavefunction of physical significance, i.e. the manifestation of the intrinsic instability of the 1D metal [ 81. On the other hand it cannot be ruled out that symmetry broken band states are a consequence of the shortcomings of one-electron band descriptions. In order to check the physical significance of symmetry broken band states numerical approaches beyond the independent-particle description are necessary. In the past two principally different techniques have been employed to investigate this fundamental
0009-2614/93/$ 06.00 Q 1993 Elsevier Science Publishers B.V. All rights reserved.
31
Volume202, number 1,2
15January 1993
CHEMICALPHYSICSLETTERS
The first approach is the so-called renormalization group (RG) method [ 91, which is an analytic model. This method has been adopted successfully by one of us to calculate the Peierls transition temperature in 1D metals [ 10,111. With the advent of high-speed computers it became possible to study CDWs in 1D chains by Feynman pathintegral quantum Monte Carlo (QMC) simulations [ 12,131. Previous QMC calculations on the ID Hubbard chain prevailingly had been restricted to the half-filled band case, which is unimportant for physically interesting strongly correlated quasi 1D synthetic metals [ 10,141. These QMC simulations already lead to remarkable insight into electronic structure factors supporting the formation of the CDWs. In this Letter, we investigate the prerequisites leading to on-site CDWs in 1D metals over a larger range in the parameter space provided by the extended Hubbard model within the theoretical framework of Feynman path-integral QMC simulations. In analogy to our recent QMC approach, which was restricted to the non-CDW domain of the Hubbard chain [ 15,161, we employ the electronic charge fluctuations ( (An;)) around the corresponding mean value ( ni) and the probability of double occupancy Pi(Z) at the ith lattice site to describe the dynamic behavior of the electronic ensemble in the 1D chain. The organization of the present Letter is as follows. In section 2 the theoretical basis of the adopted QMC method is summarized. Results and an interpretation are given in section 3. The Letter ends with a resumt in section 4.
problem.
H=-tx
i,a
(cit,c;+l,,+h.c.)+UCnjrni, 1
(1) defined by the parameters t, Uand V.t stands for the kinetic hopping integral, U and V symbolize the onsite and nearest-neighbor intersite Coulomb interactions. The c& c,,, stand for the conventional electron creation, destruction operators for spin d at the ith lattice site of the monatomic 1D chain. The associated electron number operators read nra= c&c,,, and ni = 1, n,. To derive the kinetic energy part t in (1) we adopt the conventional nearestneighbor tight-binding approximation. Application of the QMC technique to the c-number Fermion field of ( 1) requires the subdivision of the Hamiltonian into two noncommuting c fields symbolized by the suboperators If1 and Hz
+0.5U(nitnjr+ni+Itn,+,,)+~~inini+*
I
(2)
>
which are both sums of two-site interactions. The fragmentation (2) is caused by the anticommuting nature of the c field and reflects the Pauli principle. Expression (2) allows us to express the partition function Z in the form of eq. (3) with intermediate state liL> chosen to be diagonal in the occupation number representation: Z=trexp( -/3If)
2. Theoretical basis
=(,,zzL, (il Iexp(-ArH)IiL)
Subsequently only the basic principles of the employed Feynman path-integral QMC algorithm are collected. For detailed information see refs. [ 12,16,17]. The checker-board worldline realization corresponds to a path-integral decomposition of the partition function Z. The discrete path-integral is then evaluated by a heat bath Monte Carlo technique. As already mentioned we use an extended Hubbard Hamiltonian
x(iLIexp(-ArIZ)IiL_,)...(iZIexp(-ArH)Iil).
38
(3) /3 stands for 1/k,T with kB denoting the Boltzmann constant, and 7 symbolizes the imaginary time interval 0 < 7-c /I, which has been divided into L subintervals of width A7=/Z/L. As we work within the
canonical ensemble Z occurs in the denominator of the expectation value (A) of any physical observable A
Volume 202, number 1,2
CHEMICAL PHYSICS LETTERS
Mew-PHII
(Aj=tr
tr(exp-/3H) =
trMw--PWI z
(4)
’
with tr denoting the trace. The summation over the intermediate states 1iL) in (3) is performed by using a heat-bath method implemented into the checkerboard worldline algorithm. Eq. (3) for Z is the result of the approximation exp( -AS)
=exp( -Az&)
+[l+@(A~ZIH~,hl)l
exp( -AS,)
(5)
defining the fragmentation into two non-overlapping sub-Hamiltonians H, and H,;[H,,Hz] is the corresponding commutator. The error in (5) is second order in the time slices Ar and can be made negligible by making A7 small enough. Key parameters in the subsequent discussion are the electronic charge fluctuations ((An!)>={nf>-(ni>2,
(6)
which are defined by the expectation value of the number operator square minus the square of the expectation value of the number operator. In eq. (7) we give the charge fluctuations in the independentparticle limit, where they are fully determined by the electron density ( ni) at the ith atomic site: ((Anf))sCF=O.5(ni)(2-(ni)).
(7)
In the eqs. (8) and (9) we relate the probabilities Pi(n) to find n = 2, I,0 electrons at the ith lattice site to the average number of electrons (ni) and number squares {nf). Eq. (10) is the normalization condition for the Pj( n) i
Tl’Pi(n)= (nf)
)
(8)
n=l
n~oPi(n)=l.O.
(10)
The discussion in section 3 frequently refers to the mean values P(2) and ( (An’)) averaged over the i sites of the Hubbard chain. We have used these pa-
15 January I993
rameters without index i to distinguish them from the individual on-site elements labeled by the site-index i. In symmetry non-broken configurations the mean-square deviations and the probability of double occupancy change cooperatively. This one-to-one relation is lifted in CDW configurations where the on-site charge separation leads to a mutual decoupling between the charge fluctuations and the probability of double occupancy. CDWs suppress the ( (An;) ), but cause an enhancement of the double occupancies Pi(2) at lattice sites i with enlarged densities (ai); see below. The formation of the on-site CDWs is simply identified by the respective densities (ni). Previous work has shown that the “sign problem” occurring in QMC simulations is smaller for expectation values like energy and density than the statistical errors [ 181. The sign problem in QMC simulations arises if the measure of the path-integrals is not positive definite. In this case the expectation value of the sign of the corresponding measure decreases exponentially as the inverse temperature increases. For some physical properties, e.g., generalized susceptibilities, neglect of the above signs can lead to questionable results. Greater numerical problems also arise in the determination of long-range correlation functions. This has been a technical boundary to restrict the present analysis to on-site CDWs. In the subsequent QMC simulations we have always performed 25000 to 100000 QMC steps and 5000 equilibration steps. For the product ATEa value ~0.1 leads to a reliable accuracy in the computational results. By employing the above conditions we reduce the statistical error in the total energy E to values smaller than 5x 10e3 eV. The numerical accuracy in the charge fluctuations is better than 0.0 1. The 1D chain has been modelled by 48 lattice sites with cyclic boundary conditions. This construct in the QMC simulations allowed us to reproduce the results of the most comprehensive QMC calculations [ 191. Technical prerequisite of the present simulation study is an optimum program design. As a function of technical operations a width exceeding a factor of 10 in the CPU time demand is possible. A detailed description of the present QMC implementations will be reported elsewhere [ 16 I. One QMC run takes 1500 CPU s on an IBM 3090/200 E machine. Below we present the results of 110 QMC sim39
Volume202, number I,2
15January 1993
CHEMICALPHYSICSLETTERS
ulations. With the inclusion of test calculations and repeated runs to check the convergence of the method we have used more than 70 h of computer time on the above IBM equipment.
0.6 ([Ar?)) 0.4
0.2
3. Results and discussion Before discussing the present QMC results the theoretical findings of an article by Lieb and Wu [ 201 published roughly 20 years ago should be mentioned. The two authors have demonstrated that even the symmetry non-broken ground state of a strict 1D material with a half-filled dispersion is always nonmetallic for any finite two-electron interaction beyond the mean-field one. This result has been observed in the simple Hubbard model (Cl, t). Previous QMC simulations have shown that half-filling of the band is an absolutely restrictive condition [ 15,161. The non-metallic behaviour is caused by a discontinuity in the chemical potential p at ( ni} = 1.0, always implying the property of an insulator. In our previous QMC calculations the nonmetallic nature of (ni) = 1.0 1D chains has been verified by changes in the slope of the ( (An f ) > and Pi( 2) curves at ( &) = 1.0. In the following model calculations we have adopted a constant on-site repulsion U of 2.1 eV, a value which is characteristic for many quasi 1D metals formed by organometallic stacking units [21]. The t numbers employed amount to 0.25 and 1.0 eV. The first parameter is again typical for the sizeably correlated quasi 1D synthetic metals. The nearestneighbor intersite interaction has been modified as V= U, 0.75 U, 0.5 U, 0.25 U, 0, respectively. If the effective interaction V,, in the 1D Hubbard chain is symbolically written as Vefl= U- 2 V we have an attractive two electron interaction for V= U and v=o.75u [9,10]. In fig. 1 we have displayed the ( (An:) ) curves for the series with weaker electronic correlations (t= I .O eV) as a function of the average on-site density (n) , which coincides with the band occupancy in the non-dimerized Hubbard chain, The QMC simulations predict weaker on-site CDWs in the halffilled band sector (n) = 1.0 for the V= U and V=0.75Useries. All other QMC simulations lead to an even symmetry non-broken charge distribution 40
0’ 0
I
0.5
I
I
1
15
(n)
2
Fig. 1. Averagechargefluctuations ( (An‘) ) as a function of the average on-site density (n), From top ( 1) to bottom (5) we have employed the following V/ Uratios: 1.O,0.75, 0.5,0.25 and 0; this sequence refers to the (n) points with the maximum observed fluctuations. t amounts to 1.Oand U to 2.1 eV.
((n} = ( ni) at any center i). The on-site CDW is accompanied by a sizeable suppression in the charge fluctuations ( (An:)). In the density interval from 0.5~(n)
Volume 202, number 42
CHEMICAL PHYSICS LETTERS
1 PIZ)
0.8
0.6 0.k 0.2 0
Fig. 2. Average probability of double occupancy P(2) for the curves 1 to 5 of fig. 1.The top diagram corresponds to t= 1.0 eV, the bottom one to t=0.25 eV.
charge fluctuations. This proportionality is lifted for the respective mean values P( 2) and ( (An’) ) derived for on-site CDWs. The mutual atomic charge separation causes an overall suppression of the electronic charge fluctuations ( (An’) ). Consequently enhanced double occupancy is accompanied by suppressed averaged ( (An’) ) numbers. This behavior becomes evident when comparing fig. 2 with fig. 1. It could be expected a priori that the Pi( 2) curve for V=OSU would be continuous over the whole studied (n,) interval. The non-metallic ground state at (n)=l.O for the V=OSU and V=O model is accompanied by an additional localization of the charge carriers. The associated microscopic outcome is a reduction in the P, (2) numbers. In the lower half of fig. 2 we have summarized the P(2) curves for the strongly correlated series with t~0.25 eV. Within the (n) region between 0.5 and 1.5 the QMC simulations predict a CDW ground state for V= CJand V= 0.7 5U. Again this is accompanied by a P( 2) enhancement. The V=O.5 U curve is of particular importance. Only for (n} =0.5 and 1.5 the QMC simulation converges into an on-site CDW leading to a change in the slope of the corresponding P( 2) curve at the latter (n) points. Moving away from (n) = 0.5 and 1.5 towards the half-
I5 January 1993
filled band sector (ni} = I .Oa on-site symmetry nonbroken state is restored in the Hubbard chain. The {n) = 1.0 point finally indicates a non-metallic ground state caused by the aforementioned p discontinuity. The Y=O.5U curve in the stronger correlated series is the only example which exhibits three crossover points in the ground state of the 1D chain as a function of the average electron density. This clearly indicates that one quarter, half and threequarters filled 1D Hubbard chains occupy an exceptional position in their ability to form on-site CDWs. Microscopically, the CDW formation is here strongly supported by pinning effects. For a recent explanation of this phenomenon derived in the renormalization group formalism see refs. [ 10,111. It has been emphasized previously by one of us that the charge fluctuations ( (An:)) accessible for chains with ( ni) =0.5 or 1.5 at any center i correspond to a maximum also in the perfectly correlated limit describing mixed valency [ lo]. This maximum in the ( (An’)) numbers is reversed into the opposite in the case of on-site CDWs coupled to an averaged density (n} =0.5 or 1.5, where the fluctuations are completely suppressed. The corresponding behavior is symbolized in fig. 3 where we have displayed the on-site densities (ni> in two typical CDW configurations. The top representation refers to V= U, tc0.25 eV and (n) = 1.5 and shows a perfectly or-
2
(n,)- 2
In,) I.5 (n,)=l
1
0.5
1
H
Fig. 3. Schematic representation of on-site CDWs. The diagrams give the respective on-side density ( ni) at the ith center. We have used the following parameters of the Hubbard Hamiltonian Crc2.1 eV V=U, t~0.25 eV, (n)=l.S (top), (n)=l.O (bottom).
41
Volume202, number I ,2
CHEMICALPHYSICSLETTERS
dered on-site CDW with strictly alternating centers i occupied by two or one electron. The bottom diagram refers to the on-site populations derived for a half-tilled band. The vand t parameters coincide with the above example. In contrast to the ( ni) = 1.5 lattice the CDW is less ordered here. The strong alternation between ( ni) = 0 and (n,} = 2.0 is perturbed by sites with ( ni) of about 1.0. Fig. 3 is characteristic for the present QMC simulations. We always find maximum CDW ordering at ( n) = 0.5, 1.5 and perturbations in the (ni} = 1.0 case. In fig. 4 we have collected the total energy E for the two t series studied. The upper diagram corresponds to tz0.25 eV and the bottom one to t= 1.0 eV. Except for the V= 0.5 U dispersions all E curves exhibit a change in their slope at ( ni) = 1.O in the t= 1.0 series. The electronic origin of this behavior has been explained in detail. The results of the present QMC simulations can be summarized as follows. We would like to mention that most of the subsequent topics are new. In some cases they corroborate and generalize previous information already available in the literature. (i) The formation of on-site CDWs is only pos-
-2’ 0
I 0.5
I 1
I 15
i 2
(lx Fig. 4. Total energy E of the one-dimensional Hubbard chain as a function of the average on-site density (n). From top to bottom we have adopted the following V/V ratios: 1.0, 0.75, 0.5, 0.25 and 0, respectively. Upper representation: f= 1.0 eV, lower one t = 0.25 eV. Cralways refers to 2. I eV. In the inlaid diagrams the total energy E around { n ) = 1 .O is shown in an enlarged scale.
42
15January 1993
sible in a density interval between 0.5 and 1.5. For lower and higher net populations no on-site CDW stabilization has been observed. (ii) As an underlying prerequisite the following inequalities are necessary: V>2t, V20.5U. The inequality between V and 2t has been derived by additional model calculations not summarized here in detail. The ratio between Vand t is a strict boundary for all densities ( ni) # 0; it is attenuated in the halftilled band case. Previous constructions of phase diagrams of 1D systems by using the RG approach have been possible only in the phase spanned by U and V. (iii) In extension to the results derived by Lieb and Wu in the simply Hubbard scheme (U, t) we have demonstrated that even a half-tilled interacting system can be compatible with a metallic ground state. The necessary conditions are electronic correlations weaker than a characteristic threshold value and U= 2 V. (iv) As a function of the U/V ratio and the average on-site density we have derived new stability domains/phase transitions in the extended Hubbard chain. For U= 2 V and strong correlations we have demonstrated that the (n) ~0.5 and 1.5 points exhibit a singular behavior as they stabilize on-site CDWs. A third phase transition is detected at ( ni) = 1.Owhere the system is an insulator even in connection with a symmetry-non-broken wavefunction (discontinuity of the chemical potential). (v) For the first time it has been shown generally that on-side CDW formation in 1D chains occurs via doubling of the unit cell dimension irrespective of the average density. We have thus detected a new type of instability in the 1D chain at wave vectors that are not simply related to the Fermi vector. (vi) Finally we emphasize that the exceptional properties of 1D chains with one-quarter and three quarters tilled conduction bands in their ability to stabilize CDWs has been demonstrated in many experimental and previous theoretical work [ l3,10,11,14].
4. Final resume By using a Feynman path-integral QMC approach we have investigated those material properties of the strict ID Hubbard-chain which support the forma-
Volume202,number 1,2
CHEMICAL PHYSICS LETTERS
tion of on-site CDWs. Important key conditions have already been summarized in section 3 and will not be reviewed again. In section 1 we have mentioned that one reason to perform the present series of computer simulations is the hope to judge the validity of CDW type solutions derived by one-determinantal mean-field crystal orbital (CO) methods for halffilled 1D chains. Of course there is no direct one-toone correspondence between conventional CO models and the extended Hubbard scheme used in the present QMC approach. Nevertheless it is possible to interpolate qualitatively between the two stages of theoretical sophistication. Semi-empirical CO calculations have shown that the CDW formation is already enhanced with decreasing hopping elements in the mean-field representation [ 4 1. This agrees with the present “exact” QMC data; see the different behaviour of the weaker (t= 1.0 eV) and stronger (t= 0.25 eV) correlated model series. The exceptional position of half-filled bands in the stronger correlated regime has been emphasized. Here V slightly larger than 0.5U is a sufficient condition for the formation of CDWs. Although the CDW is here less ordered than in the one-quarter and three quarters filled band occupancy they are more highly ordered than for any other fractional band filling. Adoption of on-site and intersite Coulomb repulsion integrals typically found in quasi 1D synthetic metals strongly supports the assumption that one-determinantal CO solutions of the CDW type are of physical relevance in this class of compounds [ 4,5].
Acknowledgement
This work has been supported by the Deutsche Forschungsgemeinschaft through a Heisenberg grant (MCB), the Fonds der Chemischen Industrie, a doctoral fellowship of the TH Darmstadt (JS) and the Departemento de Postgrado y Especializacion (LU).
15 January 1993
One of us (MCB) is grateful to Dr. H. Kijppel and Dr. M.Yu. Lavrenziev for many enlightening discussions. The drawings kindly have been prepared by Mr. G. Wolf.
References [l] L.N. Bulaevskii, Advan. Phys. 37 (1988) 433. [ 21 International Conference on Science and Technology of Synthetic Metals, Tiibingen 1990, Synth. Met. 41-43 (1991). [ 31 D. Jerome and H.J. Schultz, Advan. Phys. 3 I ( 1982) 399. [4] M.C. Biihm, J. Chem. Phys. 81 (1984) 855; Z. Naturforsch. 39a (1984) 807; Z. Physik B 56 (1984) 99. [ 51 M.Yu. Lavrenziev, H. Kiippe1andM.C. Biihm, Chem. Phys., in press. [6]P.-0. Liiwdin, Rev. Mod. Phys. 35 (1963) 496; Advan. Chem. Phys. 14 (1969) 283, and references therein. [7] J. Katriel, Intern. J. Quantum. Chem. 6 (1972) 541; L. Pietronero and F. Zirilli, Nuovo Cimento 24B ( 1974) 57. [8] G.A. Toombs, Pbys. Rept. 40 (1978) 181. [9] J. Solyom, Advan. Phys. 28 (I 979) 201. [lO]M.C.BBhm,J.Chem.Pbys.94 (1991) 5631. [ 111M.C. BGhmand A. Staib, Chem. Phys. 155 ( 1991) 27. [ 121J.E. Hirsch and D.J. Scalapino, Phys. Rev. B 29 ( 1984) 5554. [ 131J.E. Hirsch, in: Low-dimensional conductors and superconductors, eds. D, Jerome and L.G. Caron (Plenum Press, New York, 1987) p. 7 I. [ 141MC. Biihm, J. Chem. Phys. 95 (1991) 9300. [ 151L. Utrera, J. Schulte and M.C. BShm, Chem. Phys. Letters 191 (1992) 299. 1161M.C. Btihm, J. Schulte and L. Utrera, Mol. Phys., in press. 1171J.E. Hirsch, D.J. Scalapina, R.L. Sugar and R. Blankenbecler, Phys. Rev. B 26 (1982) 5033. [ 181E.Y. Lob Jr., J.E.Gubenatis, R.T. scalettar, S.R.White, D.J. Scalapino and R.L. Sugar, Phys. Rev. B 41 (1990) 9301I [ 191W.R. Somsky, D.K. Campbell, H.-Q. Lin, X. Wang and J.E. Gubematis, Synth. Metals41-43 (1991) 3531. [20] E.H. Lieb and F.Y. Wu, Phys. Rev. Letters 20 (1968) 1445. 1211 MC. Bijhm, Lecture notes in chemistry, Vol. 45. Onedimensional organometallic materials (Springer, Berlin, 1987).
43
CHEMICAL PHYSICS LETTERS
15 January 1993
On-site charge-density wave formation in one-dimensional chains investigated by Feynman path-integral quantum Monte Carlo calculations Michael C. Bbhm, Joachim Schulte InstitutftirPhysikalische Chemie, Physikafische Chemie IIL Technische Hochschule Darmsfadt, W-6100 Darmstadt, Germany
and Luis Utrera Institute de Ciencia Materiales, Consejo Superior de Investigationes Cieniijicas, 28006 Madrid, Spain
Received 9 June 1992;in final form 22 October 1992
The formation of on-site charge-density waves (CDWs) in the one-dimensional ( ID) Hubbard chain is studied by path-integral quantum Monte Carlo (QMC) simulations as a function of the band fdling, the ratio between on-site and nearest-neighbor intersite interaction as well as the net Coulomb interaction and kinetic hopping. The dynamics of the electrons of the 1D system are quantified by the charge fluctuations ( (An;)) and the probability of double occupancy P(2) averaged over all lattice sites i. Electronic conditions favouring CDWs are summarized. It is shown that one quarter, one half and three quarter tilled bands are optimum conditions for on-site CDWs.
1. Introduction The non-conventional material properties of quasi one-dimensional ( 1D) synthetic metals have adsorbed the interest of experimentalists and theorists for more than two decades now [ I ,2]. The rich variety of accessible low-temperature ground state configurations has been a challenging problem in the design of 1D systems with specific key quantities. Quasi 1D metals can undergo phase transitions either to charge/spin density wave (C/SDW), or to singlet/ triplet superconducting (S/TS) low-temperature modifications [3]. The symmetry of the electronic wavefunction in non-distorted C/SDW states is lower than the symmetry of the corresponding (1D) lattice. For recent one-determinantal band structure calculations of 1D materials which converged into on-site and intersite CDW configurations, see refs. 74,51. In finite electronic systems symmetry broken single determinantal mean-field solutions always indicate the shortcomings of an independent-particle
description. These so-called Hartree-Fock instabilities have been studied extensively in the past years [6,7]. Note that labels like one-electron model, independent-particle scheme, Hartree-Fock approximation, etc. are used to symbolize wavefunctions described by a single Slater determinant. In comparison with symmetry broken mean-field wavefunctions observed in finite molecules, the situation is less definite if one-electron band-structure calculations of 1D solids converge into CDW type solutions. On one hand they may indicate a wavefunction of physical significance, i.e. the manifestation of the intrinsic instability of the 1D metal [ 81. On the other hand it cannot be ruled out that symmetry broken band states are a consequence of the shortcomings of one-electron band descriptions. In order to check the physical significance of symmetry broken band states numerical approaches beyond the independent-particle description are necessary. In the past two principally different techniques have been employed to investigate this fundamental
0009-2614/93/$ 06.00 Q 1993 Elsevier Science Publishers B.V. All rights reserved.
31
Volume202, number 1,2
15January 1993
CHEMICALPHYSICSLETTERS
The first approach is the so-called renormalization group (RG) method [ 91, which is an analytic model. This method has been adopted successfully by one of us to calculate the Peierls transition temperature in 1D metals [ 10,111. With the advent of high-speed computers it became possible to study CDWs in 1D chains by Feynman pathintegral quantum Monte Carlo (QMC) simulations [ 12,131. Previous QMC calculations on the ID Hubbard chain prevailingly had been restricted to the half-filled band case, which is unimportant for physically interesting strongly correlated quasi 1D synthetic metals [ 10,141. These QMC simulations already lead to remarkable insight into electronic structure factors supporting the formation of the CDWs. In this Letter, we investigate the prerequisites leading to on-site CDWs in 1D metals over a larger range in the parameter space provided by the extended Hubbard model within the theoretical framework of Feynman path-integral QMC simulations. In analogy to our recent QMC approach, which was restricted to the non-CDW domain of the Hubbard chain [ 15,161, we employ the electronic charge fluctuations ( (An;)) around the corresponding mean value ( ni) and the probability of double occupancy Pi(Z) at the ith lattice site to describe the dynamic behavior of the electronic ensemble in the 1D chain. The organization of the present Letter is as follows. In section 2 the theoretical basis of the adopted QMC method is summarized. Results and an interpretation are given in section 3. The Letter ends with a resumt in section 4.
problem.
H=-tx
i,a
(cit,c;+l,,+h.c.)+UCnjrni, 1
(1) defined by the parameters t, Uand V.t stands for the kinetic hopping integral, U and V symbolize the onsite and nearest-neighbor intersite Coulomb interactions. The c& c,,, stand for the conventional electron creation, destruction operators for spin d at the ith lattice site of the monatomic 1D chain. The associated electron number operators read nra= c&c,,, and ni = 1, n,. To derive the kinetic energy part t in (1) we adopt the conventional nearestneighbor tight-binding approximation. Application of the QMC technique to the c-number Fermion field of ( 1) requires the subdivision of the Hamiltonian into two noncommuting c fields symbolized by the suboperators If1 and Hz
+0.5U(nitnjr+ni+Itn,+,,)+~~inini+*
I
(2)
>
which are both sums of two-site interactions. The fragmentation (2) is caused by the anticommuting nature of the c field and reflects the Pauli principle. Expression (2) allows us to express the partition function Z in the form of eq. (3) with intermediate state liL> chosen to be diagonal in the occupation number representation: Z=trexp( -/3If)
2. Theoretical basis
=(,,zzL, (il Iexp(-ArH)IiL)
Subsequently only the basic principles of the employed Feynman path-integral QMC algorithm are collected. For detailed information see refs. [ 12,16,17]. The checker-board worldline realization corresponds to a path-integral decomposition of the partition function Z. The discrete path-integral is then evaluated by a heat bath Monte Carlo technique. As already mentioned we use an extended Hubbard Hamiltonian
x(iLIexp(-ArIZ)IiL_,)...(iZIexp(-ArH)Iil).
38
(3) /3 stands for 1/k,T with kB denoting the Boltzmann constant, and 7 symbolizes the imaginary time interval 0 < 7-c /I, which has been divided into L subintervals of width A7=/Z/L. As we work within the
canonical ensemble Z occurs in the denominator of the expectation value (A) of any physical observable A
Volume 202, number 1,2
CHEMICAL PHYSICS LETTERS
Mew-PHII
(Aj=tr
tr(exp-/3H) =
trMw--PWI z
(4)
’
with tr denoting the trace. The summation over the intermediate states 1iL) in (3) is performed by using a heat-bath method implemented into the checkerboard worldline algorithm. Eq. (3) for Z is the result of the approximation exp( -AS)
=exp( -Az&)
+[l+@(A~ZIH~,hl)l
exp( -AS,)
(5)
defining the fragmentation into two non-overlapping sub-Hamiltonians H, and H,;[H,,Hz] is the corresponding commutator. The error in (5) is second order in the time slices Ar and can be made negligible by making A7 small enough. Key parameters in the subsequent discussion are the electronic charge fluctuations ((An!)>={nf>-(ni>2,
(6)
which are defined by the expectation value of the number operator square minus the square of the expectation value of the number operator. In eq. (7) we give the charge fluctuations in the independentparticle limit, where they are fully determined by the electron density ( ni) at the ith atomic site: ((Anf))sCF=O.5(ni)(2-(ni)).
(7)
In the eqs. (8) and (9) we relate the probabilities Pi(n) to find n = 2, I,0 electrons at the ith lattice site to the average number of electrons (ni) and number squares {nf). Eq. (10) is the normalization condition for the Pj( n) i
Tl’Pi(n)= (nf)
)
(8)
n=l
n~oPi(n)=l.O.
(10)
The discussion in section 3 frequently refers to the mean values P(2) and ( (An’)) averaged over the i sites of the Hubbard chain. We have used these pa-
15 January I993
rameters without index i to distinguish them from the individual on-site elements labeled by the site-index i. In symmetry non-broken configurations the mean-square deviations and the probability of double occupancy change cooperatively. This one-to-one relation is lifted in CDW configurations where the on-site charge separation leads to a mutual decoupling between the charge fluctuations and the probability of double occupancy. CDWs suppress the ( (An;) ), but cause an enhancement of the double occupancies Pi(2) at lattice sites i with enlarged densities (ai); see below. The formation of the on-site CDWs is simply identified by the respective densities (ni). Previous work has shown that the “sign problem” occurring in QMC simulations is smaller for expectation values like energy and density than the statistical errors [ 181. The sign problem in QMC simulations arises if the measure of the path-integrals is not positive definite. In this case the expectation value of the sign of the corresponding measure decreases exponentially as the inverse temperature increases. For some physical properties, e.g., generalized susceptibilities, neglect of the above signs can lead to questionable results. Greater numerical problems also arise in the determination of long-range correlation functions. This has been a technical boundary to restrict the present analysis to on-site CDWs. In the subsequent QMC simulations we have always performed 25000 to 100000 QMC steps and 5000 equilibration steps. For the product ATEa value ~0.1 leads to a reliable accuracy in the computational results. By employing the above conditions we reduce the statistical error in the total energy E to values smaller than 5x 10e3 eV. The numerical accuracy in the charge fluctuations is better than 0.0 1. The 1D chain has been modelled by 48 lattice sites with cyclic boundary conditions. This construct in the QMC simulations allowed us to reproduce the results of the most comprehensive QMC calculations [ 191. Technical prerequisite of the present simulation study is an optimum program design. As a function of technical operations a width exceeding a factor of 10 in the CPU time demand is possible. A detailed description of the present QMC implementations will be reported elsewhere [ 16 I. One QMC run takes 1500 CPU s on an IBM 3090/200 E machine. Below we present the results of 110 QMC sim39
Volume202, number I,2
15January 1993
CHEMICALPHYSICSLETTERS
ulations. With the inclusion of test calculations and repeated runs to check the convergence of the method we have used more than 70 h of computer time on the above IBM equipment.
0.6 ([Ar?)) 0.4
0.2
3. Results and discussion Before discussing the present QMC results the theoretical findings of an article by Lieb and Wu [ 201 published roughly 20 years ago should be mentioned. The two authors have demonstrated that even the symmetry non-broken ground state of a strict 1D material with a half-filled dispersion is always nonmetallic for any finite two-electron interaction beyond the mean-field one. This result has been observed in the simple Hubbard model (Cl, t). Previous QMC simulations have shown that half-filling of the band is an absolutely restrictive condition [ 15,161. The non-metallic behaviour is caused by a discontinuity in the chemical potential p at ( ni} = 1.0, always implying the property of an insulator. In our previous QMC calculations the nonmetallic nature of (ni) = 1.0 1D chains has been verified by changes in the slope of the ( (An f ) > and Pi( 2) curves at ( &) = 1.0. In the following model calculations we have adopted a constant on-site repulsion U of 2.1 eV, a value which is characteristic for many quasi 1D metals formed by organometallic stacking units [21]. The t numbers employed amount to 0.25 and 1.0 eV. The first parameter is again typical for the sizeably correlated quasi 1D synthetic metals. The nearestneighbor intersite interaction has been modified as V= U, 0.75 U, 0.5 U, 0.25 U, 0, respectively. If the effective interaction V,, in the 1D Hubbard chain is symbolically written as Vefl= U- 2 V we have an attractive two electron interaction for V= U and v=o.75u [9,10]. In fig. 1 we have displayed the ( (An:) ) curves for the series with weaker electronic correlations (t= I .O eV) as a function of the average on-site density (n) , which coincides with the band occupancy in the non-dimerized Hubbard chain, The QMC simulations predict weaker on-site CDWs in the halffilled band sector (n) = 1.0 for the V= U and V=0.75Useries. All other QMC simulations lead to an even symmetry non-broken charge distribution 40
0’ 0
I
0.5
I
I
1
15
(n)
2
Fig. 1. Averagechargefluctuations ( (An‘) ) as a function of the average on-site density (n), From top ( 1) to bottom (5) we have employed the following V/ Uratios: 1.O,0.75, 0.5,0.25 and 0; this sequence refers to the (n) points with the maximum observed fluctuations. t amounts to 1.Oand U to 2.1 eV.
((n} = ( ni) at any center i). The on-site CDW is accompanied by a sizeable suppression in the charge fluctuations ( (An:)). In the density interval from 0.5~(n)
Volume 202, number 42
CHEMICAL PHYSICS LETTERS
1 PIZ)
0.8
0.6 0.k 0.2 0
Fig. 2. Average probability of double occupancy P(2) for the curves 1 to 5 of fig. 1.The top diagram corresponds to t= 1.0 eV, the bottom one to t=0.25 eV.
charge fluctuations. This proportionality is lifted for the respective mean values P( 2) and ( (An’) ) derived for on-site CDWs. The mutual atomic charge separation causes an overall suppression of the electronic charge fluctuations ( (An’) ). Consequently enhanced double occupancy is accompanied by suppressed averaged ( (An’) ) numbers. This behavior becomes evident when comparing fig. 2 with fig. 1. It could be expected a priori that the Pi( 2) curve for V=OSU would be continuous over the whole studied (n,) interval. The non-metallic ground state at (n)=l.O for the V=OSU and V=O model is accompanied by an additional localization of the charge carriers. The associated microscopic outcome is a reduction in the P, (2) numbers. In the lower half of fig. 2 we have summarized the P(2) curves for the strongly correlated series with t~0.25 eV. Within the (n) region between 0.5 and 1.5 the QMC simulations predict a CDW ground state for V= CJand V= 0.7 5U. Again this is accompanied by a P( 2) enhancement. The V=O.5 U curve is of particular importance. Only for (n} =0.5 and 1.5 the QMC simulation converges into an on-site CDW leading to a change in the slope of the corresponding P( 2) curve at the latter (n) points. Moving away from (n) = 0.5 and 1.5 towards the half-
I5 January 1993
filled band sector (ni} = I .Oa on-site symmetry nonbroken state is restored in the Hubbard chain. The {n) = 1.0 point finally indicates a non-metallic ground state caused by the aforementioned p discontinuity. The Y=O.5U curve in the stronger correlated series is the only example which exhibits three crossover points in the ground state of the 1D chain as a function of the average electron density. This clearly indicates that one quarter, half and threequarters filled 1D Hubbard chains occupy an exceptional position in their ability to form on-site CDWs. Microscopically, the CDW formation is here strongly supported by pinning effects. For a recent explanation of this phenomenon derived in the renormalization group formalism see refs. [ 10,111. It has been emphasized previously by one of us that the charge fluctuations ( (An:)) accessible for chains with ( ni) =0.5 or 1.5 at any center i correspond to a maximum also in the perfectly correlated limit describing mixed valency [ lo]. This maximum in the ( (An’)) numbers is reversed into the opposite in the case of on-site CDWs coupled to an averaged density (n} =0.5 or 1.5, where the fluctuations are completely suppressed. The corresponding behavior is symbolized in fig. 3 where we have displayed the on-site densities (ni> in two typical CDW configurations. The top representation refers to V= U, tc0.25 eV and (n) = 1.5 and shows a perfectly or-
2
(n,)- 2
In,) I.5 (n,)=l
1
0.5
1
H
Fig. 3. Schematic representation of on-site CDWs. The diagrams give the respective on-side density ( ni) at the ith center. We have used the following parameters of the Hubbard Hamiltonian Crc2.1 eV V=U, t~0.25 eV, (n)=l.S (top), (n)=l.O (bottom).
41
Volume202, number I ,2
CHEMICALPHYSICSLETTERS
dered on-site CDW with strictly alternating centers i occupied by two or one electron. The bottom diagram refers to the on-site populations derived for a half-tilled band. The vand t parameters coincide with the above example. In contrast to the ( ni) = 1.5 lattice the CDW is less ordered here. The strong alternation between ( ni) = 0 and (n,} = 2.0 is perturbed by sites with ( ni) of about 1.0. Fig. 3 is characteristic for the present QMC simulations. We always find maximum CDW ordering at ( n) = 0.5, 1.5 and perturbations in the (ni} = 1.0 case. In fig. 4 we have collected the total energy E for the two t series studied. The upper diagram corresponds to tz0.25 eV and the bottom one to t= 1.0 eV. Except for the V= 0.5 U dispersions all E curves exhibit a change in their slope at ( ni) = 1.O in the t= 1.0 series. The electronic origin of this behavior has been explained in detail. The results of the present QMC simulations can be summarized as follows. We would like to mention that most of the subsequent topics are new. In some cases they corroborate and generalize previous information already available in the literature. (i) The formation of on-site CDWs is only pos-
-2’ 0
I 0.5
I 1
I 15
i 2
(lx Fig. 4. Total energy E of the one-dimensional Hubbard chain as a function of the average on-site density (n). From top to bottom we have adopted the following V/V ratios: 1.0, 0.75, 0.5, 0.25 and 0, respectively. Upper representation: f= 1.0 eV, lower one t = 0.25 eV. Cralways refers to 2. I eV. In the inlaid diagrams the total energy E around { n ) = 1 .O is shown in an enlarged scale.
42
15January 1993
sible in a density interval between 0.5 and 1.5. For lower and higher net populations no on-site CDW stabilization has been observed. (ii) As an underlying prerequisite the following inequalities are necessary: V>2t, V20.5U. The inequality between V and 2t has been derived by additional model calculations not summarized here in detail. The ratio between Vand t is a strict boundary for all densities ( ni) # 0; it is attenuated in the halftilled band case. Previous constructions of phase diagrams of 1D systems by using the RG approach have been possible only in the phase spanned by U and V. (iii) In extension to the results derived by Lieb and Wu in the simply Hubbard scheme (U, t) we have demonstrated that even a half-tilled interacting system can be compatible with a metallic ground state. The necessary conditions are electronic correlations weaker than a characteristic threshold value and U= 2 V. (iv) As a function of the U/V ratio and the average on-site density we have derived new stability domains/phase transitions in the extended Hubbard chain. For U= 2 V and strong correlations we have demonstrated that the (n) ~0.5 and 1.5 points exhibit a singular behavior as they stabilize on-site CDWs. A third phase transition is detected at ( ni) = 1.Owhere the system is an insulator even in connection with a symmetry-non-broken wavefunction (discontinuity of the chemical potential). (v) For the first time it has been shown generally that on-side CDW formation in 1D chains occurs via doubling of the unit cell dimension irrespective of the average density. We have thus detected a new type of instability in the 1D chain at wave vectors that are not simply related to the Fermi vector. (vi) Finally we emphasize that the exceptional properties of 1D chains with one-quarter and three quarters tilled conduction bands in their ability to stabilize CDWs has been demonstrated in many experimental and previous theoretical work [ l3,10,11,14].
4. Final resume By using a Feynman path-integral QMC approach we have investigated those material properties of the strict ID Hubbard-chain which support the forma-
Volume202,number 1,2
CHEMICAL PHYSICS LETTERS
tion of on-site CDWs. Important key conditions have already been summarized in section 3 and will not be reviewed again. In section 1 we have mentioned that one reason to perform the present series of computer simulations is the hope to judge the validity of CDW type solutions derived by one-determinantal mean-field crystal orbital (CO) methods for halffilled 1D chains. Of course there is no direct one-toone correspondence between conventional CO models and the extended Hubbard scheme used in the present QMC approach. Nevertheless it is possible to interpolate qualitatively between the two stages of theoretical sophistication. Semi-empirical CO calculations have shown that the CDW formation is already enhanced with decreasing hopping elements in the mean-field representation [ 4 1. This agrees with the present “exact” QMC data; see the different behaviour of the weaker (t= 1.0 eV) and stronger (t= 0.25 eV) correlated model series. The exceptional position of half-filled bands in the stronger correlated regime has been emphasized. Here V slightly larger than 0.5U is a sufficient condition for the formation of CDWs. Although the CDW is here less ordered than in the one-quarter and three quarters filled band occupancy they are more highly ordered than for any other fractional band filling. Adoption of on-site and intersite Coulomb repulsion integrals typically found in quasi 1D synthetic metals strongly supports the assumption that one-determinantal CO solutions of the CDW type are of physical relevance in this class of compounds [ 4,5].
Acknowledgement
This work has been supported by the Deutsche Forschungsgemeinschaft through a Heisenberg grant (MCB), the Fonds der Chemischen Industrie, a doctoral fellowship of the TH Darmstadt (JS) and the Departemento de Postgrado y Especializacion (LU).
15 January 1993
One of us (MCB) is grateful to Dr. H. Kijppel and Dr. M.Yu. Lavrenziev for many enlightening discussions. The drawings kindly have been prepared by Mr. G. Wolf.
References [l] L.N. Bulaevskii, Advan. Phys. 37 (1988) 433. [ 21 International Conference on Science and Technology of Synthetic Metals, Tiibingen 1990, Synth. Met. 41-43 (1991). [ 31 D. Jerome and H.J. Schultz, Advan. Phys. 3 I ( 1982) 399. [4] M.C. Biihm, J. Chem. Phys. 81 (1984) 855; Z. Naturforsch. 39a (1984) 807; Z. Physik B 56 (1984) 99. [ 51 M.Yu. Lavrenziev, H. Kiippe1andM.C. Biihm, Chem. Phys., in press. [6]P.-0. Liiwdin, Rev. Mod. Phys. 35 (1963) 496; Advan. Chem. Phys. 14 (1969) 283, and references therein. [7] J. Katriel, Intern. J. Quantum. Chem. 6 (1972) 541; L. Pietronero and F. Zirilli, Nuovo Cimento 24B ( 1974) 57. [8] G.A. Toombs, Pbys. Rept. 40 (1978) 181. [9] J. Solyom, Advan. Phys. 28 (I 979) 201. [lO]M.C.BBhm,J.Chem.Pbys.94 (1991) 5631. [ 111M.C. BGhmand A. Staib, Chem. Phys. 155 ( 1991) 27. [ 121J.E. Hirsch and D.J. Scalapino, Phys. Rev. B 29 ( 1984) 5554. [ 131J.E. Hirsch, in: Low-dimensional conductors and superconductors, eds. D, Jerome and L.G. Caron (Plenum Press, New York, 1987) p. 7 I. [ 141MC. Biihm, J. Chem. Phys. 95 (1991) 9300. [ 151L. Utrera, J. Schulte and M.C. BShm, Chem. Phys. Letters 191 (1992) 299. 1161M.C. Btihm, J. Schulte and L. Utrera, Mol. Phys., in press. 1171J.E. Hirsch, D.J. Scalapina, R.L. Sugar and R. Blankenbecler, Phys. Rev. B 26 (1982) 5033. [ 181E.Y. Lob Jr., J.E.Gubenatis, R.T. scalettar, S.R.White, D.J. Scalapino and R.L. Sugar, Phys. Rev. B 41 (1990) 9301I [ 191W.R. Somsky, D.K. Campbell, H.-Q. Lin, X. Wang and J.E. Gubematis, Synth. Metals41-43 (1991) 3531. [20] E.H. Lieb and F.Y. Wu, Phys. Rev. Letters 20 (1968) 1445. 1211 MC. Bijhm, Lecture notes in chemistry, Vol. 45. Onedimensional organometallic materials (Springer, Berlin, 1987).
43