# Competition between charge and spin order in the extended Hubbard model on the triangular lattice.

## Abstract

Several new classes of compounds can be modeled in first approximation by electrons on the triangular lattice that interact through on-site repulsion as well as nearest-neighbor repulsion . This extended Hubbard model on a triangular lattice has been studied mostly in the strong coupling limit for only a few types of instabilities. Using the extended two-particle self consistent approach (ETPSC), that is valid at weak to intermediate coupling, we present an unbiased study of the density and interaction dependent crossover diagram for spin and charge density wave instabilities of the normal state at arbitrary wave vector. When dominates over and electron filling is large, instabilities are chiefly in the spin sector and are controlled mostly by Fermi surface properties. Increasing eventually leads to charge instabilities. In the latter case, it is mostly the wave vector dependence of the vertex that determines the wave vector of the instability rather than Fermi surface properties. At small filling, non-trivial instabilities appear only beyond the weak coupling limit. There again, charge density wave instabilities are favored over a wide range of dopings by large at wave vectors corresponding to superlattice in real space. Commensurate fillings do not play a special role for this instability. Increasing leads to competition with ferromagnetism. At negative values of or , neglecting superconducting fluctuations, one finds that charge instabilities are favored. In general, the crossover diagram presents a rich variety of instabilities. We also show that thermal charge-density wave fluctuations in the renormalized classical regime can open a pseudogap in the single-particle spectral weight, just as spin or superconducting fluctuations.

###### pacs:

71.10.Fd, 05.30.Fk, 71.10.-w## I Introduction

One of the outstanding problems in quantum many-body physics is to understand quasi-two dimensional systems where both electron-electron interaction and geometric frustration are important (1). The triangular lattice is the prime example where the geometry frustrates near-neighbor anti-alignment of the spins that naturally tends to occur in the presence of short-range electron-electron interaction. Studying models of interacting electrons on such a lattice is thus certainly of fundamental interest, but it is also strongly motivated by the discovery of new materials. Prime examples of these materials are organic bis-(ethylenedithio) (BEDT) Cu(CN) layered compounds (2), triangular lattice antiferromagnets of the CuCrO family (3), and transition-metal oxide materials like and . The layered cobaltates have drawn much attention because of their unconventional properties. Sodium cobaltate shows an unusually strong thermopower(4) at doping that can be suppressed drastically by applying an in-plan magnetic field(5). The observation of Curie-Weiss behavior in the magnetic susceptibility while resistivity displays metallic behavior is another puzzle(7). The system also becomes superconductor when it is diluted by water(8); (9); (10). Various types of charge- and spin orders also have been found in the system for various dopings(11); (12); (13).

consist of two-dimensional layers separated by insulating layers. The layers have Co atoms at the center of oxygen octahedra forming a 2D triangular lattice. The band structure calculation performed by Singh(14), revealed details of splitting of the bands in Co atoms. With help of this calculation and also of NMR experiments(7), one can find a rough estimate of hopping and exchange constants that would enter a two-dimensional Hubbard or model for this system. However, the modeling is complicated by the fact that band structure calculations lead to hole pockets that are not observed experimentally, a question that is still debated by several groups using, for example, the Gutzwiller approximation (15), the local density approximation plus Hubbard (16) and dynamical-mean field theory (17); (18); (19). In addition, the effect of long-range Coulomb interaction from the sodium leads to modifications to the simplest Hubbard Hamiltonian for the cobaltates(20); (21).

In this paper, we do not address the question of detailed modeling of the cobaltates or of other triangular lattice systems. Instead, we note that since several types of spin and charge density waves are observed in these materials, it is quite likely that first-neighbor repulsion , and not only on-site repulsion , must be taken into account. by itself tends to favor spin-density waves. We thus just focus on the simplest extended one-band model Hubbard model on the triangular lattice and ask a few general questions: What types of phases are typical in different doping ranges, what type of interaction favors them, and should one expect pseudogap effects.

Previous theoretical and numerical works have obtained phase diagrams for the triangular lattice in the presence of competing interactions. There are, for example, variational Monte-Carlo calculations (23); (22) for the extended Hubbard model. That work focused mostly on the presence of the Charge density wave (CDW) at filling and RVB superconductivity at . Slave boson methods were used for the and models (24); (25) to study CDW, ferromagnetism and also RVB superconductivity. Series expansion methods and cluster mean field theory (26) have also investigated CDW, Néel order, ferromagnetic order, dimer order and phase separation in a model. We will comment further on some of these calculations in the context of our own results.

The results of this paper are obtained with the recently developed Extended two-particle self-consistent approach (27); (28)(ETPSC) that is valid from weak to intermediate coupling. This method has been benchmarked against Quantum Monte-Carlo (QMC) simulations (QMC) for the extended Hubbard model on a square lattice. The approach satisfies conservation laws and the Mermin-Wagner theorem stating that no continuous symmetry can be broken at finite temperature in two dimensions. More traditional methods, such as the Random phase approximation, do not satisfy this requirement. With ETPSC, quantum renormalization of interactions (Kanamori-Brüuckner screening) is taken into account. Instability towards zero-temperature long-range order is signaled at finite temperature by crossover to the renormalized-classical regime where the correlation length grows exponentially. The wave-vector of the instability is determined self-consistently by the approach and all wave vectors are in principle allowed. No a priori selection is necessary.

Within ETPSC we can also compute the self-energy and other related quantities, such as the spectral weight that is measured in photoemission experiments (31). For the Hubbard model, it has been shown with the Two-particle self-consistent approach (TPSC) that a pseudogap can appear as precursor induced by either antiferromagnetic (31) or superconducting fluctuations (31); (33). The former (32) has been observed experimentally in electron-doped high-temperature superconductors (34). Our results demonstrate that CDW fluctuations can also induce a pseudogap. This is a relevant question experimentally given that CDW induced pseudogaps are observed and sometimes even show similarities with observations in high-temperature superconductors (35).

In the following we first introduce the model and the ETPSC methodology. We next present our numerical results, discussing various physical effects in terms of the spin and charge structure factors. We display phase diagrams that help understand how microscopic parameters favor various phases. The CDW induced pseudogap and its effect on the Fermi surface are discussed before we present an overview and a conclusion.

## Ii Model and method

We write the extended Hubbard Hamiltonian in the following form,

(1) |

where () are annihilation (creation) operators for electrons of spin at site of a triangular lattice, is the density operator, and is the hopping matrix element. The quantities and are the on-site and nearest-neighbor interactions respectively and is the chemical potential.

For the Hubbard model (), TPSC is a very reliable approach up to intermediate coupling limit. The functional derivative method is particularly convenient to obtain the TPSC equations (36). This is the method that was used to generalize TPSC to the extended Hubbard model (27); (28).

The equations that need to be solved are the following. The charge and spin response functions take the form

(2) |

(3) |

where

(4) |

are the charge and spin vertex functions and is the free response function (non-interacting susceptibility) given by

(5) |

with

(6) |

the non-interacting dispersion relation and . In the above formula is the volume of the Brillouin zone (BZ), is the Fermi function and is the non-interacting chemical potential. The pair correlation functions are related to the static structure factors by

(7) |

(8) |

where are the charge and spin component of the static structure factor. The spin resolved static structure factor is defined by and is the Fourier transform of . The quantities and entering the vertices Eq.(II) are simply the pair correlation functions at the first-neighbor distance.

Self-consistency is established by connecting the static structure factors to the response functions by the fluctuation-dissipation theorem

(9) |

where are Bosonic Matsubara frequency. Substituting the expression for the susceptibilities Eqs.(2 and 3) and the corresponding vertices Eq.(II) on the right-hand side, one can use the result to obtain the pair correlation functions and entering the vertices using their relations Eqs.(7 and 8) to the structure factors. Assuming that the functional derivatives of the pair correlation functions are known, as discussed below, we need only three equations to determine the pair correlation functions entering the vertices. This is because the Pauli principle imposes that . The equation that is dropped out is that for . This procedure and its impact on the Pauli principle is discussed in detail in Ref. (28).

Functional derivatives of the pair correlation functions with respect to density and magnetization enter the spin and charge vertices. The functional derivatives are obtained from the following equations:

(10) | |||

These equations are strictly valid only when particle-hole symmetry is satisfied. On the square lattice, it has been checked by comparisons with QMC calculations that the results are satisfactory even in the absence of the full particle-hole symmetry (28). Apparently, particle-hole symmetry due to linearization of the dispersion relations near the Fermi surface suffices. We will make this assumption for the triangular lattice where strict particle-hole symmetry is not satisfied. This is justified a posteriori by our results. Those that can be checked against variational QMC, for example, are in excellent agreement.

Finally, the self-energy needed to address the pseudogap problem is obtained following Ref. (27).

(11) |

where is the Fermionic Matsubara frequency and is the Bosonic one. We can also find the spectral function . The above formula does not assume a Migdal theorem since one of the vertices is renormalized. However, it takes into account only the longitudinal fluctuations. Transverse fluctuations could be accounted for following a generalization of the steps in Ref. (37). Since the pseudogap appears only when fluctuations are large, the longitudinal case suffices to establish the qualitative results.

Finally, the interacting chemical potential is obtained from

(12) |

## Iii Response functions

The “phase diagrams” that we present in the section that follows this one, are determined from the growth of the the spin and charge response functions as temperature decreases. When the interaction is local, , the wave vector of the instability is determined entirely from the non-interacting susceptibility, in other words from nesting properties of the Fermi surface. The introduction of near-neighbor repulsion changes this since it introduces a wave vector dependence to the vertices. In order to disentangle the various effects, we present the non-interacting susceptibilities in the first subsection and the results with interactions in the second subsection. These numerical results are obtained from Eqs. (3-II).

### iii.1 Non-interacting susceptibility

The non-interacting susceptibility Eq. 5 is determined mostly by the shape of the non-interacting Fermi surface that is in turn determined by the dispersion relation given in Eq. (6). In Fig. 1 we present the Fermi surface for increasing values of the density, respectively. The first Brillouin zone is plotted as a solid line. It is important to notice that the Fermi surface corresponding to (long dash) touches the first Brillouin zone boundary and has long parallel segments that lead to near nesting. We will show in more details that this causes a strong peak in the non-interacting response function. We also draw two wave-vectors that are often found for the most important charge or spin density waves in the system. At these wave-vectors, the charge or spin response functions often have a strong peak. The real-space modulations corresponding to these the wave-vectors are depicted in Fig. 2.

The non-interacting response function Eq. 5 is drawn in Fig. 3 for different values of densities at . The largest response is for . While one might have expected that parallel segments of the Fermi surface would have lead to a peak at a single dominant wave vector, it seems that the frustration imposes a less pronounced maximum. However, the height of the maximum at that density increases rapidly with decreasing temperature. The position of all the peaks changes only slightly with temperature. The height of the peaks for the smaller values of density does not change drastically with decreasing temperature. That fact in addition to quantum renormalization of the interactions are the main reasons for the absence of any instability at low density up to intermediate coupling.

There is a deep minimum near the K point at higher values of the density. This is the main reason for absence of commensurate spin density wave (SDW). The non-interacting response function has a peak at the commensurate wave-vector (K point) at lower value of the density but, as we just mentioned, this peak does not grow enough to produce any sort of order including SDW up to intermediate coupling. That is not the case in the strong coupling limit but that is out of reach of our approach in the density regime where Mott physics is dominant.

It is quite remarkable that the free response function shows a strong peak at the origin for , a signature for ferromagnetism at nearby densities. It seems that, as we will see, frustration on the triangular lattice favors ferromagnetism at intermediate coupling, contrary to the square lattice case.

### iii.2 Interacting response functions

In the presence of both types of interactions and , the response functions are strongly modified. Consider typical values of the interaction, and . Using the same color code as in the figure for the non-interacting case, we show in Fig. 4 the spin, Eq.(3), and in Fig. 5 the charge, Eq.(2), response functions at for the same densities as in the non-interacting case Fig. 3.

In the ordinary random phase approximation, the spin response function in influenced only by the interaction and the maxima are at the same location as the in non-interacting case. In the present approach however, the nearest-neighbor interaction also influences the spin response, introducing a wave vector dependent vertex. Hence, some of the maxima of the spin response function in Fig. 4 are not at the same wave vector as in the non-interacting case. Nevertheless, the differences are much smaller than for the charge response function appearing in Fig. 5. In the latter case, the position of the maximum is near point K (wave vector ) for all densities, in other words the charge response is dominated by the wave vector dependent vertex introduced by .

One can summarize the results for the location of the maximum spin response function as follows. In the range the tendency is towards an incommensurate spin-density wave (ISDW) while the response is maximum near zero wave vector (ferromagnetism) just above this density. Generally the height of the peaks increases when is reduced, hence nearest-neighbor repulsion does not favor spin order. For , the maximum is near point K corresponding to the same lattice structure as the CDW2 depicted in Fig. 2 except that one should replace the big or small points with up and down spins. This resembles the spin structure that would arise with ferrimagnetic order.

For the charge response in Fig. 5, tendency towards CDW2 order (K point) is robust for these values of and . The tendency is strongest at low values of the density because of a weaker effect of the frustration that leads to a dip in the non-interacting response function in Fig. 3. The effect of frustration is very strong for densities very close to . But nevertheless the height of the peak for densities around grows dramatically with decreasing temperature.

We verify the dominant effects of and discussed above, this time by fixing the filling at and changing the interactions. Fig. 6 and Fig. 7 show, respectively, the effect on the charge and spin response functions at .

For the charge response function in Fig. 6, we imply for those curves where the value of is not written. A simple comparison of the charge response function with the corresponding non-interacting susceptibility for in Fig. 3 shows the importance of the term in the vertex Eq. (II). The non-interacting susceptibility has a deep minimum on the Brillouin zone boundary. The charge response function on the other hand shows two different maxima at wave-vectors and . The CDW modulation related to these wave vectors are illustrated in Fig. 2. We will see that these instabilities occur over a wide area in the plane. In fact apart from the phase separation instability (), which occurs for negative , they are the only charge instabilities that can be found at this density. The situation changes as we change the density and one can expect to find incommensurate CDW instabilities in another region of the plane.

Moving on to the spin susceptibility, Fig. 7 shows that the presence of the term suppresses the spin response function more and more as increases, concomitant with the increase in the charge response function. In principle, one cannot find a strong maximum in both the spin and charge response functions. This is true in all one band homogeneous paramagnetic systems as dictated by Eqs. (7), (8) and (Pauli sum-rule). Indeed, the Pauli sum-rule (obtained from the Pauli principle ) connects to in such a way that an increase in one forces a decrease in the other (30).

## Iv Crossover diagrams

In mean-field theory, one normally finds finite temperature phase transitions, in contradiction with the Mermin-Wagner theorem. In ETPSC, we obtain instead at a temperature below which the correlation length begins to grow exponentially, diverging only at zero temperature. is lower than the mean-field transition temperature because of the quantum Kanamori-Brückner renormalization of the vertices. In this low temperature regime, the characteristic frequency of the growing fluctuations becomes less than temperature in dimensionless units. This is the renormalized-classical regime. Either the spin or the charge correlations grow exponentially at some characteristic wave vector that suggests which long-range order will likely be stabilized at zero temperature. Since our approach is not valid deep in the renormalized classical regime, one cannot be sure that the zero-temperature phase will be precisely that suggested by the behavior at .

The value of depends on density, and . We use to estimate . For the sake of computational efficiency, we chose the constant to be and checked that the general features do not change by choosing a larger value. This occurs because the exponential growth of the correlations is rather sudden. A detailed discussion of this issue can be found in previous publications (27); (28).

We present our results for in Figs. 8 to 13 as color (grey scale) plots in various planes of parameter space. There are four plots that present the dependence of at four densities, then two plots for the dependence at fixed . We indicate by lines of various colors and types the boundaries between regions where there is either a change in the wave vector of the growing correlations, or a change in the type of correlation, spin or charge. When we indicate a paramagnetic (Fermi liquid) region (PM), we mean that correlations did not grow, in either the spin or charge sectors, at temperatures as low as .

Fig. 8 displays at as a function of and for both positive and negative values. At negative values of either or , superconducting correlations will be competing. Since superconductivity has not been taken into account here, the results in all quadrants, except the first one, should be taken as just indicative of what may happen in the spin or charge sectors. When is negative, unless is large, there is a strong tendency to phase separation (PS), i.e. the static charge response function starts, at , to grow exponentially for wave vector . At positive and , incommensurate spin density waves (ISDW) are dominant, but and must be large enough, as expected from the absence of perfect nesting. At small and , the system remains paramagnetic (PM). Charge density waves appear at positive and only if is relatively large. Recall however that the effect of is amplified by the presence of several neighbors. The charge instability in this parameter range is of the CDW1 type illustrated in Fig. 2. Charge instabilities are further amplified at positive only if is allowed to become negative. The CDW2 pattern is allowed only in extreme conditions of large positive and large negative . This is not suprising since one can see from the non-interacting Fermi surface in Fig. 1 that the corresponding wave vector is not particularly favored by nesting. The CDW2 phase is really governed by properties of the vertex , not so much by the non-interacting Fermi surface.

When density is decreased to , the non-interacting Fermi surface becomes almost circular so at positive and the tendency to order is strongly suppressed, as can be seen from Fig. 9. Compared with the previous figure, the CDW2 vertex related instability is more robust while the CDW1 and ICDW instabilities occur in smaller regions, the ICDW existing over a larger region this time than CDW1. At larger fillings, , where again the Fermi surface becomes almost circular, similar features are observed.

As the filling decreases, the Fermi surface becomes more and more circular. Restricting ourselves to positive and , and staying at weak to intermediate coupling where our theory is valid, nothing interesting occurs. The system remains paramagnetic down to . We thus also present, in Fig. 10 and Fig. 11, results at large and where our theory is not strictly controlled. We feel these results are nevertheless interesting for two reasons. First, some of our “phase boundaries” compare favorably with results obtained from other methods. Second, the renormalized classical regime occurs at such high temperature that and may begin to control the approach.

Fig. 10 and Fig. 11 thus show the crossover diagram for, respectively, and over a wide range of positive and . The CDW2 region now appears at positive and , contrary to the results in the previous figures, as long as the stabilizing interaction is large enough. The boundary that separates CDW2 from PM in Fig. 10 is very close to QMC results (23), which gives us confidence in the validity of the results. The ferromagnetic phase is dominant when both and are large. The ICDW phase does not appear at filling (Fig. 11). The CDW2 phase is influenced to some extent not only by the vertex, but also by commensurability, as can be seen from the fact that it is more important at than at . There are competing tendencies for the CDW2 phase: 1) The density is more favorable for the commensurate CDW2 as reflected in the free response function, 2) The effect of off-site interaction is less important at lower value of the density. Based on these observations, one might surmise the presence of an optimal density where CDW2 appears over a larger area of space.

To explore in more details the density dependence, we present results as a function of and at fixed . This was studied in particular by Motrunich and Lee (22); (24). They calculated the phase diagram with different methods: a) renormalized mean field theory and variational quantum Monte Carlo with a trial wave function (22) and b) slave boson mean field theory (24). They suggest that the effect of the term is taken into the account more accurately in the first method than in the second one. In the first paper, they found that the CDW2 phase can be reached at smaller at the densities and than at other densities. This is not the case in the second paper where these densities play no special role. Since their calculations are at , this suggests that the general features of the phase diagram can be understood physically at low and high density: One needs a large to stabilize the CDW at low value of the density. Since at large value of the density, the hole plays the same role as the electron at low density, we expect at large to find the CDW as well. This argument is correct if other phases do not suppress the CDW.

We present the results of calculations at finite values of in Figs. 12 and 13. In Fig. 12, the results in the plane are for . It can easily be seen that CDW2 appears when both and are large. A ICDW region separates CDW2 from the paramagnetic, or Fermi liquid phase. There is a small area on the left of the plane where the ferromagnetic phase is stable. Fig. 13 shows the results for . CDW2 is still stable in the same region of the plane but the small FM region of the previous figure has now grown and pushed away slightly the CDW2 phase. In other words, at larger , spin fluctuations are playing a more important role, as expected. This is also shown by the appearance of an ISDW regime. It is obvious from Figs. 12 and 13 that the densities and do not play any special role, at least at these values of . This is in agreement with the results of Ref. (24). However, at large , we find more ferromagnetic spin fluctuations, a possibility that was not considered in Ref. (24). We also calculated the same phase diagram at higher value of where we find that the CDW2 region is completely swept away by ferromagnetism.

## V Fermi surface and pseudogap

It has already been shown that spin and superconducting thermal fluctuations can (30); (31); (33); (32), in two dimensions, open up a pseudogap on the Fermi surface that reflects the wave vector of the fluctuations. The same study can be performed here to check the effect of charge density wave fluctuations. We obtain the self-energy by substituting our results for the susceptibilities in Eq. (II), from which one can compute the spectral function.

We begin by the effect of spin fluctuations. In Fig. 14, we plot for , , and . In the jargon, this is known as a Momentum distribution curve (MDC). The dashed Green line is the Brillouin zone. Following the largest intensity regions, one can recognize the shape of the non-interacting Fermi surface illustrated for in Fig. 1. It is clear from this figure that the effect of the on-site interaction at this temperature is just to introduce damping. Since this is a region where ISDW appear at low temperature, the Fermi surface can be destroyed by lowering the temperature or increasing .

For better understanding, we plot in Fig. 15 the spectral function as a function of frequency at different wave vectors for the same parameters as in the previous figure. These are Energy dispersion curves (EDC). One can clearly observe the quasi-particle dispersion relation, the effect of appearing as damping. The extra features at higher frequency are precursors of extra bands that would appear in the ordered state. Indeed, the EDC are plotted in a regime where the correlation length associated with incommensurate fluctuations is three to four lattice spacings (), enough to enter the renormalized classical regime. A simple generalization of an argument presented earlier (31) shows that in that regime, the self-energy can be approximated by

(13) |

where the sum runs over all the equivalent maxima of the susceptibility, six of them for the present case. Substituting this expression into the general result for the spectral weight

(14) |

one can reproduce qualitatively the behavior in Fig. 15. With larger correlation length (larger ), the spectrum clearly splits into more bands.

The spin fluctuation induced precursor effects are also illustrated in Fig. 16 that displays the density of states for various values of and for densities near . The curve is close to the non-interacting result. For stronger interaction, the extra band of states can be understood as a precursor effect by drawing the density of states for a static modulation.

The fact that renormalized classical spin fluctuations can create a pseudogap on the Fermi surface in two dimensions is well documented so we do not show a figure corresponding to this case. We present the case of a Fermi surface pseudogap that originates from charge fluctuations. Increasing enough that CDW1 fluctuations become important, one can observe a pseudogap, as seen in for , , and displayed in Fig. 17. The pseudogap opens up at the point, in other words the spectral function at this point is suppressed. By analogy with the spin-fluctuation induced pseudogap, the symmetry equivalent points are nearly connected by the symmetry equivalent wave vectors of the CDW1 charge instability.

## Vi Discussion and conclusion

We have used ETPSC to clarify the leading instabilities of the extended Hubbard model on the triangular lattice. The method is valid from weak to intermediate coupling. It satisfies the Mermin-Wagner theorem, includes quantum (Kanamori-Brüuckner) renormalization of the vertices and does not assume a Migdal theorem for the self-energy. We interpret the entry into the renormalized classical regime (exponential growth of the correlation length) as an indication of the phase that acquires long-range order at zero temperature. Since we scan all wave vectors for both spin and charge instabilities, our method is not biased towards a restricted set of instabilities as most other studies. Superconducting instabilities, however, have not been explored.

The range of possible phases as a function of on-site interaction , nearest-neighbor interaction and filling is quite rich. Notwithstanding superconducting instabilities, negative values of and favor charge instabilities either at zero wave vector (phase separation) when dominates, or at the CDW1 and CDW2 wave vectors when dominates. In the physically more relevant regime where and are both positive, spin instabilities are favored when dominates and they occur at wave vectors that are essentially determined by the Fermi surface. That is particularly clear at large fillings where the shape of the Fermi surface is non-trivial. Larger , on the other hand, favors charge instabilities that are generally determined by the vertex itself instead of by details of the Fermi surface. This is particularly clear at small filling where the Fermi surface is essentially circular. That predominance of the vertex is why the CDW2 () phase, for example, exists over a wide range of fillings and is not favored by commensurate fillings. We find that this CDW2 phase can compete with spin instabilities, especially ferromagnetism, when is large as well. Ferromagnetic fluctuations appear in a range of doping similar to that observed for the cobaltates. That competition with ferromagnetism has not been taken into account in earlier studies (22); (24). All the phases, except the Fermi liquid one, appear at large values of the interactions when filling is small, somewhat outside the regime of validity of our approach. Nevertheless, agreement with other approaches suggest that ETPSC extrapolates in a reasonable way towards strong coupling. The disappearance of all phases except the Fermi liquid one at small coupling and small fillings () on the triangular lattice is a clear manifestation of the effects of frustration, as can be seen by contrasting with the square lattice case (28) () where this does not occur.

Finally, we also showed that charge-density wave thermal fluctuations can also induce a pseudogap. A pseudogap associated with CDW is observed (35) experimentally. If its origin is the one discussed in the present paper, it should disappear as temperature rises above that where the charge correlation length becomes of the order of the thermal de Broglie wave length (31), as observed in the spin fluctuation case in electron-doped cuprates (34).

Further studies should include the competition with superconducting fluctuations as well as more realistic modeling of the specific materials to which one wishes to apply our results. For example, the case is relevant for the layered organics.

## Vii Acknowledgments

Computations were performed on the Elix2 Beowulf cluster in Sherbrooke and on the Ms cluster of the Réseau Québécois de calcul haute performance (RQCHP). The present work was supported by the Natural Sciences and Engineering Research Council NSERC (Canada), the Fonds québécois de recherche sur la nature et la technologie FQRNT (Québec), the Canadian Foundation for Innovation CFI (Canada), the Canadian Institute for Advanced Research CIFAR, and the Tier I Canada Research Chair Program (A.-M.S.T.).

### References

- N. P. Ong and R. J. Cava, Science 305, 52 (2004).
- Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito, Phys. Rev. Lett. 91, 107001 (2003)
- S. Seki, Y. Onose, and Y. Tokura, arXiv:0801.3757
- I. Terasaki, Y. Sasago, and K. Uchinokura, Phys. Rev. B56, R12685 (1997).
- Y. Wang, n. S. Rogado, R. J. Cava, and N. P. Ong, Nature (London) 423, 425 (2003).
- K. Takada, H. Sakurai, E. Takayama-Muromachi, F. Izumi, R.A. Dilanian, T. Sasaki, Nature 422, 53 (2003).
- R. Ray, A. Ghoshray, K. Ghoshray, and S. Nakamura, Phys. Rev. B59, 9454 (1999).
- R. E. Shaak, T. Klimczuk, M. L. Foo, adn R. J. Cava, Nature 424, 527 (2003).
- F. C. Chou, J. H. Cho,, P. A. Lee, E. T. Abel, K. Matan, and Y. S. Lee, Phys. Rev. Lett. 92, 157004 (2004).
- R. Jin, B. C. Sales, P. Khalifah, and D. Mandrus, Phys. Rev. Lett. 91, 217001 (2003).
- M. L. Foo, Y. Wang, S. Watauchi, H. W. Zandbergen,T. He, R. J. Cava, and N. P. Ong, Phys. Rev. Lett. 92, 247001 (2004).
- Y. Tokura, Phys. Today 56, 50 (2003).
- A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
- D. J. Singh, Phys. Rev. B61, 13397 (2000).
- Sen Zhou, Meng Gao, Hong Ding, Patrick A. Lee, and Ziqiang Wang, Phys. Rev. Lett. 94, 206401 (2005).
- A.O. Shorikov, V.I. Anisimov, and M.M. Korshunov, arXiv:0705:1408.
- A. Liebsch, and H. Ishida, arXiv:0705:3627
- C. A. Marianetti, K. Haule, and O. Parcollet, Phys. Rev. Lett. 99, 246404 (2007).
- K. Aryanpour, W. E. Pickett, and R. T. Scalettar, Phys. Rev. B 74, 085117 (2006).
- M. Roger, D. J. P. Morris, D. A. Tennant, M. J. Gutmann, J. P. Goff, J.-U. Hoffmann, R. Feyerherm, E. Dudzik, D. Prabhakaran, A. T. Boothroyd, N. Shannon, B. Lake, and P. P. Deen, Nature 445, 631 (2007).
- M.-H. Julien, C. de Vaulx, H. Mayaffre, C. Berthier, M. Horvatić, V. Simonet, J. Wooldridge, G. Balakrishnan, M.R. Lees, D.P. Chen, C.T. Lin, P. Lejay, arXiv:0801.4095.
- O. I. Motrunich and P. A. Lee, Phys. Rev. B69, 214516 (2004).
- H. Watanabe and M. Ogata, J. Phys. Soc. Jpn. 74, 2901(2005).
- O. I. Motrunich and P. A. Lee, Phys. Rev. B70, 024514 (2004).
- G. Baskaran, Phys. Rev. Lett. 91, 097003 (2003).
- W. Zheng, J. Oitmaa, C. J. Hamer, and R. R. P. Singh, Phys. Rev. B 70, 020504(R) (2004).
- B. Davoudi and A.-M.S. Tremblay Phys. Rev. B 74, 035113 (2006).
- B. Davoudi and A.-M. S Tremblay Phys. Rev. B 76, 085115 (2007).
- Y. Zhang, and J. Callaway, Phys. Rev. B39, 9397 (1989).
- Y. M. Vilk, Liang Chen, and A.-M. S. Tremblay, Phys. Rev. B49, 13267 (1994).
- Y. M. Vilk, and A.-M. S. Tremblay, J Phys. I 7, 1309 (1997).
- B. Kyung, V. Hankevych, A.-M. Daré, A.-M. S. Tremblay, Phys. Rev. Lett. 93, 147004 (2004).
- B. Kyung, S. Allen, and A.-M. S. Tremblay, Phys. Rev. B64, 075116 (2001).
- E. Motoyama, G. Yu, I. Vishik, O. Vajk, O. Vajk, P. Mang, and M. Greven, Nature, 445, 186, (2007).
- A. A. Kordyuk, S. V. Borisenko, V. B. Zabolotnyy, R. Schuster, D. S. Inosov, R. Follath, A. Varykhalov, L. Patthey, and H. Berger, arXiv:0801:2546.
- S. Allen, A.-M. S.Tremblay,Y. M.Vilk, in Theoretical Methods for Strongly Correlated Electrons, edited by editor C. B., D. Sénéchal and editor A.-M. Tremblay (2003).
- S. Moukouri, S. Allen, F. Lemay, B. Kyung, D. Poulin, Y.M. Vilk, and A.-M. S. Tremblay, Phys. Rev. B 61, 7887 (2000)