Lattice susceptibility for 2D Hubbard Model within dual fermion method
Abstract
In this paper, we present details of the dual fermion (DF) method to study the nonlocal correction to single site DMFT. The DMFT twoparticle Green’s function is calculated using continuous time quantum monte carlo (CTQMC) method. The momentum dependence of the vertex function is analyzed and its renormalization based on the BetheSalpeter equation is performed in particlehole channel. We found a magnetic instability in both the dual and the lattice fermions. The lattice fermion susceptibility is calculated at finite temperature in this method and also in another recently proposed method, namely dynamical vertex approximation (DA). The comparison between these two methods are presented in both weak and strong coupling region. Compared to the susceptibility from quantum monte carlo (QMC) simulation, both of them gave satisfied results.
pacs:
71.10.FdI Introduction
Strongly correlated electron systems, such as heavy fermion compounds, hightemperature superconductors, have gained much attention from both theoretical and experimental point of view. The competition between the kinetic energy and strong Coulomb interaction of fermions generates a lot of fascinating phenomena. Various theoretical approaches have been developed to treat the regime of intermediate coupling. The widely used perturbative methods, such as random phase approximation (RPA), fluctuation exchange (FLEX)Bickers et al. (1989); Bickers and White (1991a), and the twoparticle selfconsistent (TPSC)Allen and Tremblay (1995); Vilk et al. (1994) method are based on the expansion in the Coulomb interaction which is only valid in weakcoupling. To go beyond the perturbative approximation and to gain insight of the correlation effects of the fermion systems, new theoretical methods are needed. Dynamical mean field theory (DMFT)Metzner and Vollhardt (1989); Hartmann (1984); George et al. (1996) is a big step forward in the understanding MetalInsulator transition.
Dynamical mean field theory maps a manybody interacting system on a lattice onto a single impurity embedded in a noninteracting bath. Such a mapping becomes exact in the limit of infinite coordination number. All local temporal fluctuations are taken into account in this theory, however spatial fluctuations are treated on the mean field level. DMFT has been proven a successful theory describing the basic physics of the MottHubbard transition. But the nonlocal correlation effect can’t always be omitted. Although, straight forward extensions of DMFTHettler et al. (1998); Kotliar et al. (2001); Okamoto et al. (2003); Potthoff (2003); Maier et al. (2005) have captured the influence of shortrange correlation, these methods are still not capable of describing the collective behavior, e.g. spin wave excitations of manybody system. At the same time, most of the numerically exact impurity solvers require a substantial amount of time to achieve a desired accuracy even on a small cluster, which makes the investigation of larger lattice to be impossible.
Recently, some efforts have been made to take the spatial fluctuations into account in different waysToschi et al. (2007); Rubtsov et al. (); Kusunose (2006); Tokar and Monnier (); Slezak et al. (). All these methods construct the nonlocal contribution of DMFT from the local twoparticle vertex. The electron selfenergy is expressed as a function of the twoparticle vertex and the singleparticle propagator. The cluster extention of DMFT considers the correlation within the small cluster. Compared to these, the diagrammatic resummation technique involved in these new methods makes them only approximately include the nonlocal corrections. While, long range correlations are also considered in these methods and the computational burden is not serious.
In this paper we will apply the method of RubtsovRubtsov et al. () to consider the vertex renormalization of the DF through the BetheSalpeter equation. Lattice susceptibility is calculated from the renormalized DF vertex.
The paper is organized as follows: In Sec. II we summarize the basic idea of the DF method and give details of the calculation. The DMFT twoparticle Green’s function and the corresponding vertex calculation are implemented in CTQMC in Sec. III. The frequency dependent vertex is modified through the BetheSalpeter Equation to obtain the momentum dependence in Sec. IV. In Sec. V we present the calculation of the lattice susceptibility and compare it with QMC results and also the works from ToschiToschi et al. (2007). The conclusions are summarized in Sec. VI, where we also present possible application.
Ii The DF method
We study the general oneband Hubbard model at two dimensions
(1) 
creates (annihilates) an electron with spin and momentum . The dispersion relation is , where is the number of lattice sites. The basic idea of the DF methodRubtsov et al. () is to transform the hopping between different sites into coupling to an auxiliary field . By doing so, each lattice site can be viewed as an isolated impurity. The interacting lattice problem is reduced to solving a multiimpurity problem which couples to the auxiliary field. This can be done using the standard DMFT calculation. After integrating out the lattice fermions one can obtain an effective theory of the auxiliary field where DMFT serves as an starting point of the expansion over the coupling between each impurity site with the auxiliary field.
To explicitly demonstrate the above idea we start from the action of DMFT which can be written as
(2) 
where is the hybridization function of the impurity problem defined by which is the action of an isolated impurity at site with the local Green’s function . Using the Gaussian identity, we decouple the lattice sites into many impurities which couple only to the field
(3)  
The equivalence of Eqs. (2) and (3) form an exact relation between the Green’s funtion of the lattice electrons and the DF.
(4) 
This relation is easily derived by considering the derivative over in the two actions. Eq. (4) allows now to solve the manybody “lattice” problem based on DMFT which is different from the straight forward cluster extension. The problem is now to solve the Green’s function of the DF . It is determined by integrating Eq. (3) over and yielding a Taylor expansion series in powers of and . The Grassmann integral ensures that and appear only in pairs associated with the lattice fermion nparticle vertex obtained from the singlesite DMFT calculation. In this paper we restrict our considerations to the twoparticle vertex .
Expanding the LuttingerWard functional in , the first two contributions to the self energy function are the diagrams shown in Fig. 1. Diagram (a) vanishes for the bare DF since this diagram exactly corresponds to the DMFT self consistency. Therefore the first nonlocal contribution is given by diagram (b). The selfenergy for these two diagrams are
(5a)  
(5b) 
Here spacetime notation is used, , . Fermionic Matsubara frequency is , bosonic frequency is where is the inverse temperature. Together with the bare DF Green’s function , the new Green’s function can be derived from the Dyson equation
(6) 
The algorithm of the whole calculation is:

Set initial value of for the first DMFT loop.

Determine the singlesite DMFT Green’s function from the hybridization function . The selfconsistency condition ensures that the first diagram of the DF selfenergy is very small.

Go through the DMFT loop once again to calculate the twoparticle Green’s function and corresponding function. The method for determining the function is implemented for both strong and weakcoupling CTQMC in the next section of this paper.

Start an inner loop calculation to determine the DF Green’s function and in the end the lattice Green’s function.

The lattice Green’s function is then given by Eq. (4) from that of the DF.

Fourier transform the momentum lattice Green’s function into real space. And from the onsite component to determine a new hybridization function which is given by Eq. (9).

Go back to the Step 3. and iteratively perform the outer loop until the hybridization doesn’t change any more.
Although diagram (a) is exactly zero for the bare DF Green’s function, it gives nonzero contribution from the second loop where the DF Green’s function is updated from Eq. (6). As a result, the hybridization function should also be updated before the next DMFT loop is performed . This is simply done by setting the local full DF Green’s function to zero, together with the condition that the old hybrization function forces the bare local DF Green’s function to be zero (), we obtain a set of equations
(7a)  
(7b) 
which yields
(8) 
This equation finally gives us the relation between the new and old hybridization function.
(9) 
In the whole calculation, the DF perturbation calculation converges quickly. The most time consuming part of this method is the DMFT calculation of the two particle Green’s function. There are some useful symmetries to accelerate the calculation. As already pointed outAbrikosov et al. (1963); Nozieres (1964), it is convenient to take the symmetric form of the interaction term. The two particle Green’s function is then a fully antisymmetric function. Such fully antisymmetric form is very useful to speed up the calculation of the two particle Green’s function. One does not need to calculate all the frequency points within the cutoff in Mastsubara space, a few special points are calculated and the values for the other points are given by that of those special points through antisymmetric property. In the DF self energy calculation, we always have the convolution type of momentum summation which is very easy to be calculated by fast fourier transform (FFT).
Iii CTQMC and twoparticle vertex
From the above analysis, the key idea of the DF method is to construct the nonlocal contribution from the auxiliary field and the DMFT twoparticle Green’s function. Therefore it is quite important to accurately determine the twoparticle vertex. Here we adapt the newly developed CTQMC methodRubtsov et al. (2005); Werner et al. (2006); Werner and Millis (2006) to calculate the two particle Green’s function .
First we briefly outline the CTQMC technique. For more details, we refer the readers toRubtsov et al. (2005); Werner et al. (2006); Werner and Millis (2006). Here we discuss the twoparticle Green’s function and some numerical implemetations in more detailed. Two variants of the CTQMC methods have been proposed based on the diagrammatic expansion. Unlike the HirschFye method, these methods don’t have a Trotter error and can approach the low temperature region easily. In the weakcoupling methodRubtsov et al. (2005) the noninteracting part of the partition function is kept and expanded the interaction term into Taylor series. Wick’s theorem ensures that the corresponding expansion can be written into a determinant at each order
(10) 
with
(11) 
where is the noninteracting action and is the Weiss field, and the oneparticle Green’s function is measused as
(12) 
In the strong coupling method the effective action is expanded in the hybridization function by integrating over the noninteracting bath degrees of freedom. Such an expansion also yields a determinant.
(13) 
Here . The action is evaluated by a Monte Carlo random walk in the space of expansion order . Therefore the corresponding hybridization matrix changes in every Monte Carlo step. One particle Green’s function is measured from the expansion of hybridization function as . is the inverse matrix of the hybridization function. Apparently one needs to calculate this inverse matrix in every update step which is time consuming, fortunately it can be obtained by the fastupdate algorithmRubtsov et al. (2005).
At the same time such a relation allows direct measurement of the Matsubara Green’s function
(14) 
Compared with the imaginary time measurement, it seems additional computational time is needed for the sum over every matrix elements . K. Haule proposed to implement such measurement in every fast update procedure which makes sure that only linear amount of time is neededHaule and Kotliar (2007).
In our calculation the Green’s function is measured in the weakcoupling CTQMC at each accepted update which greatly reduces the computational time. The weakcoupling CTQMC normally yields a higher perturbation order than the strongcoupling CTQMC. It seems that the performance of the strongcoupling CTQMC is betterGull et al. (2007). Concerning the convergence speed, the weakcoupling CTQMC is almost same as the strongcoupling one under the above implementation together with a proper choice of , since in strongcoupling CTQMC more Monte Carlo steps are needed usually in order to smooth the noise of Green’s function at imaginary time around or at large Matsubara frequency points. Furthermore, the weakcoupling CTQMC is much easier implemented for large cluster DMFT calculation, in which case the strongcoupling method needs to handle a big eigenspace. In this paper we mainly use weakcoupling CTQMC as impurity solver, while all the results can be obtained in the strongcoupling CTQMC which was used as an accuracy check.
Similarly, we adapt K. Haule’s implementation to calculate the twoparticle Green’s function in frequency space. In the weak coupling CTQMC, the noninteracting action has Gaussian form which ensures the applicability of Wick’s theorem for measuring the two particle Green’s function
(15)  
The overline indicates the Monte Carlo average. In each Monte Carlo measurement, depends on two different argument and , only in the average level, is a function of single frequency. In each fastupdate procedure, the new and old have a closed relation which ensures that one can determine the updated Green’s function from the old one . For example, adding pair of kinks and supposing before updating the perturbation order is , then it is for the new Mmatrix. The new inserted pair is at row and column.
(16)  
Here, , and have the same definition as in RefRubtsov et al. (2005). In every step, one only needs to calculates the Green’s function when the update is accepted and only a few calculations are needed. A similar procedure for removing pairs, shiftting endpoint operation can be used. Such method is also applicable in the segment picture of strongcoupling CTQMC. In the weakcoupling CTQMC, such an implementation greatly improves the calculating speed in low temperature and strong interaction regime^{1}^{1}1In fact, the improvement is more obvious for larger Mmatrices. The strong coupling CTQMC and the weak coupling CTQMC require approximately the same amount of CPU time although in the weak coupling case the average perturbation order is higher than in the strong coupling case. Once one obtains the two frequency dependent Green’s function in every monte carlo step, the twoparticle Green’s function can be determined easily from Eq. (15). The twoparticle vertex is then given from the following equation:
(17) 
where
(18) 
is the bare susceptibility. For the multiparticle Green’s function, it still can be constructed from the two frequency dependent Green’s function , but more terms appear from Wicks theorem. Simply, when set one can calculate the oneparticle Green’s funtion easily.
Iv Momentum dependece of Vertex
As mentioned earlier diagram (a) in Fig. 1 only gives the local contribution. The first nonlocal correction in the DF method is from diagram (b). Momentum dependences comes into this theory through the bubblelike diagram between the two vertices which yields the momentum dependence of the DF vertex. The natural way to renormalize vertex is through the BetheSalpeter equation. Since the DMFT vertex is only a function of Matsubara frequency, the integral over internal momentum and ensures that the full vertex only depends on the center of mass momentum . The BetheSalpeter equation in the particlehole channelAbrikosov et al. (1963); Nozieres (1964) are shown in Fig. 2.
From the construction of the DF method, we know the interaction of the DF is coming from the two particle vertex of lattice fermion which is obtained through DMFT calculation. In the BetheSalpeter equation, it plays the role of the buildingblock. The corresponding BetheSalpeter equation for these two channels are
(19a)  
(19b) 
Here, the short hand notation of spin configuration is used. represents , while is denoted by where . are the full vertices in the and channel, respectively. is the full DF Green’s function obtained from section II which is kept unchanged in the calculation of the BetheSalpeter Equation. Different from the work of S. BrenerBrener et al. (), we solve the above equations directly in momentum space with the advantage that in this way we can calculate the susceptibility for any specific center of mass momentum and it’s convenient to use FFT for investigating larger lattice.
In Eq. (19) one has to sum over the internal spin indices in the channel which is not present in channel. One can decouple the channel into the charge and spin channels which can be solved seperately, and it turns out that the spin channel vertex function is exactly same as the that in channel, see e.g. P. NozieresNozieres (1964). Such relation is true for the DMFT vertex, and was also verified for the momentum dependent vertex in the DF methodBrener et al. (). In our calculation, we have solved the channel by decoupling it to the charge and spin channel, while the channel is not used.
Once the converged momentum dependent DF vertex is obtained, one can determine the corresponding DF susceptibility in the standard way by attaching four Green’s functions to the DF vertex.
(20a)  
(20b) 
The momentum sum over and can be performed independently by FFT becasue the DF vertx only depends on the center of mass momentum .
Now the zcomponent DF spin susceptibility can be determined from the spin channel component calculated above. In Fig. 3, is shown for at temperatures (left panel) and (right panel). With the lowing down of temperature the DF susceptibility grows up, especially at wave vector . The momentum and run from to .
The susceptibility is strongly peaked at the wave vector at the low temperature case and the peak value becomes higher and higher. The magnetic instability of the DF system is indicated by the enhancement of the DF susceptiblity. The effect of momentum dependence of vertex is clearly visible in this diagram. The bare vertex which is only a function of frequency becomes momentum dependent through the BetheSalpeter equation. Later on we will see that such momentum dependent vertex plays a very important role in the calculation of the lattice fermion susceptibility.
V Lattice susceptibility
The strong antiferromagnetic fluctuation in 2D system is indicated by the enhancement of the DF susceptibility at the wave vector shown in Fig 3. This is the consequence of the deep relation between the the Green’s function of the lattice and the DF, see Eq. (4). In order to observe the magnetic instability of the lattice fermion directly, we have calculated the lattice susceptibility based on the DF method. By differentiating the partition function in Eqns. (2, 3) twice over the kinetic energy, we obtain an exact relation between the susceptibility of DF and lattice fermions. After some simplificationsBrener et al. (), it is given by
(21) 
Here cannt be interpreted as a particle propagator, it is defined as:
(22) 
Again, the sum is performed over internal momentum and frequency which is performed by FFT and rough summing over a few Matsubara points. Again as in Eq. (4), this equation established a connection between the lattice susceptibility and the DF susceptibility. From this point of view, it is easy to understand that the instability of DFs generates the instability of the lattice fermions.
One can also find relations for the higher order Green’s function of the DF and the lattice fermions in the same way. This emphasizes the similar nature of the DF and lattice fermions except that DF possess only nonlocal information, since the DMFT selfconsistency ensures that the local DF Green’s function is exactly zero.
The lattice magnetic susceptibility is calculated using the following definition
(23)  
where . represents the lattice susceptibility in order to distinguish with that of the DF.
We have used two different ways to calculate the lattice susceptibility. First we have solved the above equation using the bare vertex which is obtained from the DMFT calculation. In contrast, the second calculation was performed using the full DF vertex. In both of these calculations, the full one particle DF Green’s function was used. The momentum dependent of the DF vertex is obtained through the calculation of the BetheSalpeter equation. The lattice susceptibility is expected to be improved if we use the momentum dependence DF vertex. In this way, we can understand the effect of momentum dependence in the DF vertex.
In Fig. 4 we plotted the results for the uniform susceptibility by using both the bare and full DF vertex. The lattice QMC resultMoreo (1993) is shown for comparison. The calculation is done for and several values of temperature. The momentum sum is approximated over 32 32 points here. Both of these calculations reproduce the well known CurieWeiss law behavior. Surprisingly enough, the results for the bare vertex fit the QMC results better than that for the momentum dependent vertex. We believe that this is the finite size effect of QMCMoreo (1993). A. Moreo showed that becomes smaller when increasing the cluster size . The 4 4 cluster calculation result at the same temperature located above of that from 8 8 cluster calculation. Therefore the results obtained from the full vertex is expected to be more reliable.
The importance of the momentum dependence of the DF vertex is more clearly observed in the calculation of , see Fig. 5. Again, in this diagram QMC resultsBickers and White (1991b) are shown for comparison. The same parameters are used as in Fig. 4. The result from the DF with bare vertex does not produce the same results compared with QMC solution. Evenmore interesting, with decreasing temperature the deviation becomes larger. On the other hand, the momentum dependent vertex in the DF method gives a satisfactory answer. This shows the importance of the momentum dependence in the DF vertex function. Fig. 6 shows the evolution of against for fixed transfer frequency . The path in momentum space is shown in the inset. From this diagram we can see that reaches its maximum value at wave vector .
The comparison between the DF and QMC results shows the good performance of DF method. The DF calculation started from a single site DMFT calculation and by introducing an auxiliary field, the nonlocal information is introduced and nicely reproduces the QMC results. Our calculation could be done within four hours for each value of the temperature on average. In this sense, this method is cheap and reliable compared with the more computationally intensive lattice QMC calculation.
Similar as the DF method, Dynamical Vertex Approximation (DA)Toschi et al. (2007) is also based on the two particle local vertex. It deals with the lattice fermion directly, without introducing any auxiliary field. The perturbative nature of this method ensures its validity at weakcoupling regime. Unlike in the DF method, DA takes the irreducible two particle local vertex as building blocks.
(24a)  
(24b) 
with the spin and charge vertex defined as . The bare susceptibility is defined as
(25a)  
(25b) 
And the selfenergy is calculated through the standard SchwingerDyson equation
(26) 
Here, the full vertex is obtained by summing all the channel dependent vertices and subtracting the double counted diagrams.
(27)  
The one particle propagator is given by the DMFT lattice Green’s function where the self energy is purely local , the local Green’s function is . Then the Dyson equation gives the lattice Green’s function from the selfenergy function .
Before presenting the comparison, we take a deeper look at the analysis of Eq. (24),
(28)  
The second term in the brackets on RHS removes the local term from the bare susceptibility. The whole term in the brackets then represents only the nonlocal bare susceptibility. In order to compare with the DF method, we take the inverse form of Eq. (19)
(29)  
The above two equations are same except for the last term. Since the local DF Green’s function is zero, the bare DF susceptibility is purely nonlocal which coincides with the analysis of DA BetheSalpeter equation. Therefore, it is not surprising that these two methods generate similar results. It is not easy to perform a term to term comparison between the DF method and DA although the bare susceptibilities have no local term in both of these method. The one particle Green’s functions have different meaning in these two methods.
The lattice susceptibility within the DA method is obtained by attaching four Green’s functions on the vertex obtained in Eq. (27). There are two possible choices of the lattice Green’s function, one is the DMFT lattice Green’s function , the other one is the Green’s function constructed by the nonlocal selfenergy from the Dyson equation. In Fig. 7 and 8, we presented the DA lattice susceptibility calculated from both the DMFT lattice Green’s function labeled as DA() and the full Green’s function labeled as DA(). The DF result from the calculation with the full DF vertex is replotted for comparison. In Fig. 7, the DA susceptibility calculated from the DMFT Green’s function (DA()) is basically the same as the DF susceptibility only with some small deviation. The results for which are not shown here which nicely repeat the DF and QMC results, the deviation between the DA and the DF method becomes smaller with the increasing of temperature. The DA susceptibility is calculated from the full Green’s function (DA()) shows a different behavior at low temperature regime which reached its maximum value at . As we know, the Hubbard Model at half filling with strong coupling maps to the Heisenberg model, reasches a maximum at where is the effective spin coupling constant given as . The calculation uses the parameter which is in the intermediate coupling regime. Therefore we further calculated the lattice susceptibility at which are shown in Fig. 9.
When the temperature is greater than 0.4, the DF method and DA (DA()) generate the similar results to the QMC calculation. Reducing the temperature further, the QMC susceptibility greatly drops and shows a peak around 0.4 which coincides with the behavior of the Heisenberg model. The DF femion and DA susceptibility continuously grows up with the decreasing of temperature. Although the DA with the full Green’s function (DA()) shows a peak, it locates at which is larger than the peak position of the QMC. And DA() generated a large deviation from that of QMC. In this diagram, we only show the results of the DF approach for and the DA results for . The Bethesalpeter equation of the DA have a eigenvalue approaching one when further lowering the temperature, which makes the access of lower temperature region impossible.
Fig. 8 shows the results of DA susceptibilities at wave vector . In contrast to the comparison for results, the DA susceptibility calculated from the full Green’s function DA () yields better results than that from the calculation with the DMFT Green’s function DA (). DA () results are almost on top of the DF results, the results with DMFT Green’s function DA () is large than the DF results. The deviation becomes larger at lower temperature. Summarizing, the DA calculation using the full Green’s function generated the same result as the DF method for while failed to produce correctly. In contrast, the calculation with the DMFT Green’s function in DA nicely produced the results calculated with the DF method for while generated larger devivation for at lower temperature regime compared to that from the DF method. Together with Fig. 4 and 5, we can see that the DF fermion calculation with the full DF vertex generated basically the same results for both and compared to the results of QMC.
In both the DF method and the DA, the operation of inverting large matrices is required for solving the BetheSalpeter equation. Fig. 10 shows the leading eigenvalue of Eqns. (19) and (24). As expected, the leading eigenvalue approaches one with decreasing temperature which directly indicates the magnetic instability of 2D system. The eigenvalues corresponding to the DF fermion method always lies below of that from DA indicating the better convergence of the DF method. When the leading eigenvalues are closed to one, the matrix inversion in Eqns. (19) and (24) are ill defined, which prevents the investigation at very low temperature.
Concerning the performance of the DF method, we also calculated the uniform susceptibility at away halffilling. In the strongcoupling limit, the Hubbard model is equivalent to the Heisenberg model with coupling constant . The consequence of doping is to effectively decrease the coupling , which yields the increasing behavior of with doping. The finite size QMC calulationMoreo (1993); Chen and Tremblay (1993) observed a slightly increasing with very small doping at strong interaction or in the low temperature region. Here, we did a similar calculation at and . Since the DF method and the DA do not suffer from the finite size problem. We would expect to observe results similar to those of QMCMoreo (1993); Chen and Tremblay (1993). In DA the suseceptibility is calculated from the DMFT Green’s function and the vertex obtained from Eq. (27). As shown in Fig. 11 at , the susceptibility slightly increases in the weak dopping region where is around , DF fermion results clearly showed such behavior, DA also gave a signal of it. Further doping the system, both the DA and the DF method reproduce the decrease with doping as already seen in the QMC. With the increasing of interaction, we would expect to see the enhancement of this effect, however our calculation indicates that such increasingdecreasing behaviro dissappear. Both the DA and the DF method give the same decreasing curve which contradict to QMC resultMoreo (1993). The results will most likely be further improved by including the higher order vertex or calculating the cluster DMFT plus DF/DAHafermann et al. (2007).
Vi Conclusion
In this paper, we extended both the DF method and DA to calculate the lattice susceptibility. Both of these methods gave equally good results compared with QMC calculation at . Although they are supposed to be weakcoupling methods, at these two methods generated right results at high temperature region. While both of them failed to reproduce the Heisenberg physics at low temperature. The investigation of the lattice susceptibility suffers from hard determined matrix inversion problem at low temperature regime. The DF methods always generates smaller eigenvalues compared to DA indicating the better convergence. The implementation of DF method in momentum space greatly improves the calculational speed and makes it easier to deal with larger size lattice.
Acknowledgements.
We would like to thank the condensed matter group of A. Lichtenstein at Hamburg University for their hospitality in particular for the discussions and open exchange of data with H. Hafermann. Gang Li and Hunpyo Lee would like to thank Philipp Werner for his help in implementing the strongcoupling CTQMC code.References
 Bickers et al. (1989) N. E. Bickers, D. J. Scalapino, and S. R. White, Phys. Rev. Lett. 62, 961 (1989).
 Bickers and White (1991a) N. E. Bickers and S. R. White, Phys. Rev. B 43, 8044 (1991a).
 Allen and Tremblay (1995) S. Allen and A.M. Tremblay, J. Phys. Chem. Solids (UK) 56, 1769 (1995).
 Vilk et al. (1994) Y. Vilk, L. Chen, and A.M. Tremblay, Phys. Rev. B 49, 13267 (1994).
 George et al. (1996) A. George, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
 Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
 Hartmann (1984) E. M. Hartmann, Z. Phys. B 57, 281 (1984).
 Maier et al. (2005) T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
 Hettler et al. (1998) M. H. Hettler, A. N. TahvildarZaden, M. Jarrell, T. Pruschke, and H. R. Krishnamurthy, Phys. Rev. B 58, 7475 (1998).
 Kotliar et al. (2001) G. Kotliar, S. Savrasov, G. Pallson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
 Okamoto et al. (2003) S. Okamoto, A. Millis, H. Monien, and A. Fuhrmann, Phys. Rev. B 68, 195121 (2003).
 Potthoff (2003) M. Potthoff, Eur. Phys. J. B 36, 335 (2003).
 Toschi et al. (2007) A. Toschi, A. A. Katanin, and K. Held, Phys. Rev. B 75, 045118 (2007).
 (14) A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, eprint condmat/0612196.
 Kusunose (2006) H. Kusunose, J. Phys. Soc. Jpn. 75, 054713 (2006).
 (16) V. I. Tokar and R. Monnier, eprint condmat/0702011.
 (17) C. Slezak, M. Jarrell, T. Maier, and J. Deisz, eprint condmat/0603421.
 Abrikosov et al. (1963) A. Abrikosov, L. Gorkov, and I. Dzyaloshinski, Methods of quantum field theory in statistical physics (Dover, New York, N.Y., 1963).
 Nozieres (1964) P. Nozieres, Theory of interacting fermi systems (W. A. Benjamin, INC., New York. New York, 1964).
 Rubtsov et al. (2005) A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
 Werner et al. (2006) P. Werner, A. Comanac, L. D. Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
 Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
 Haule and Kotliar (2007) K. Haule and G. Kotliar, Phys. Rev. B 75, 155113 (2007).
 Gull et al. (2007) E. Gull, P. Werner, A. Millis, and M. Troyer, Phys. Rev. B 76, 235123 (2007).
 (25) S. Brener, H. Hafermann, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, eprint condmatt/0711.3647.
 Moreo (1993) A. Moreo, Phys. Rev. B 48, 3380 (1993).
 Bickers and White (1991b) N. E. Bickers and S. R. White, Phys. Rev. B 43, 8044 (1991b).
 Chen and Tremblay (1993) L. Chen and A.M. S. Tremblay, Phys. Rev. B 49, 4338 (1993).
 Hafermann et al. (2007) H. Hafermann, S. Brener, A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, Pis’ma v ZhETF 86, 769 (2007).