Superconductivity from repulsion: Variational results for the Hubbard model in the limit of weak interaction
Abstract
The twodimensional Hubbard model is studied for small values of the interaction strength ( of the order of the hopping amplitude ), using a variational ansatz well suited for this regime. The wave function, a refined Gutzwiller ansatz, has a BCS meanfield state with wave symmetry as its reference state. Superconducting order is found for densities (but not for ). This resolves a discrepancy between results obtained with the functional renormalization group, which do predict superconducting order for small values of , and numerical simulations, which did not find superconductivity for . Both the gap parameter and the order parameter have a domelike shape as a function of with a maximum for . Expectation values for the energy, the particle number and the superconducting order parameter are calculated using a linkedcluster expansion up to second order in . In this way large systems (millions of sites) can be readily treated and well converged results are obtained. The gap parameter exhibits an intriguing powerlaw increase as a function of . Finitesize effects are very important for lattices with , in agreement with Anderson’s criterion for the minimum size of a superconductor, but even for irregularities persist because the level spacing at the Fermi energy (the HOMOLUMO gap) fluctuates as a function of for the tightbinding band structure.
I Introduction
The Hubbard model was introduced more than fifty years ago by Anderson Anderson_59 (), Gutzwiller Gutzwiller_63 (), Hubbard Hubbard_63 () and Kanamori Kanamori_63 () for discussing magnetic ordering in narrowband materials. Later on the model was used for describing the Mott metalinsulator transition Brinkman_70 () and it even served as a microscopic basis for Landau’s Fermi Liquid theory of He Vollhardt_84 (). A dramatic upsurge of interest set in when Anderson postulated that the Hubbard model on a square lattice embodies the essentials of cuprate superconductors, by reproducing both the insulating antiferromagnetic phase of the parent compounds and the superconducting phase observed upon doping Anderson_87 (). Still today, it is not clear to what extent Anderson’s postulate is corroborated by experiments.
That superconductivity can arise from purely repulsive interactions was already shown in 1965 by Kohn and Luttinger Kohn_65 (). Using perturbation theory for a continuum model with weak shortrange and purely repulsive interaction, they found superconductivity for an unconventional symmetry of the order parameter. During the last thirty years a lot of effort has been spent to find out whether pairing by repulsion is also realized in the twodimensional Hubbard model, for positive values of the onsite coupling parameter .
For large values of there is a general consensus that superconductivity, preferentially with wave symmetry, exists for a density close to one particle per site, i.e., for a nearly halffilled band Lee_06 (); Scalapino_12 (). Variational methods have been particularly helpful for estimating the energy gap and the superconducting order parameter, and for studying the competition between antiferromagnetism and superconductivity Edegger_07 (); Ogata_08 (). Cluster extensions of Dynamical MeanField Theory have also found evidence for (wave) pairing for intermediate to large values of Maier_05 (), but different schemes yield rather different results Rohringer_18 ().
For small values of , the method of the functional renormalization group has been particularly helpful for tracking the instabilities of the twodimensional Hubbard model Metzner_12 (). Initially, superconducting instabilities were detected through divergences of the effective pairing interaction in the normal phase Zanchi_00 (); Halboth_00b (); Honerkamp_01c (); Binz_03 (). More recently, techniques were developed to continue the procedure into the superconducting phase, either by using “partial bosonization” Friederich_11 () or by combining the scaling in the normal phase with meanfield theory in the ordered phase Eberlein_14 (). While calculations based on the functional renormalization group consistently predict superconductivity (wave pairing close to half filling), they cannot provide quantitative results for the energy gap or for the order parameter.
Variational calculations do make quantitative predictions, but, besides from being to some extent biased, they usually rely on Monte Carlo sampling, therefore they are limited to modest system sizes (typically sites for a square lattice) and suffer from statistical uncertainties. This should not be a problem if the energy gap due to superconductivity is large enough, i.e., for moderate to large values of . However, for small values of , where variational Monte Carlo calculations are unable to find any evidence for superconductivity, one has to worry about the reliability of the method.
In Gutzwiller’s celebrated variational ansatz Gutzwiller_63 () an operator acts on an uncorrelated state to reduce double occupancy. A simple modification refines this ansatz by adding an operator involving the kinetic energy Baeriswyl_87 (); Otsuka_92 (). The refined wave function is particularly well suited for treating the small limit Otsuka_92 (); Dzierzawa_95 (); Tahara_08 (). Superconductivity of the twodimensional Hubbard model has been explored by applying variational Monte Carlo to this ansatz Eichenberger_07 (); Yanagisawa_16 ().
In the present paper a slightly modified version of the refined Gutzwiller wave function is worked out perturbatively. Our method is restricted to a relatively small region of coupling strengths, , where is the hopping amplitude between nearestneighbor sites. Nevertheless, the results are clearcut, especially so because very large system sizes can be treated easily (millions of sites). Superconductivity is found away from half filling (but not at half filling), with a largely increased order parameter as compared to results obtained with the standard Gutzwiller ansatz. The dependence of the gap parameter is consistent with a power law. This dependence differs from what has been found in RPAtype theories for spinfluctuation exchange. We also find that superconductivity is not necessarily produced by a lowering of potential energy, but depending on filling it may also be due to a lowering of the “kinetic energy” (the expectation value of the hopping term). Therefore nonBCS features are not necessarily a privilege of the large region of the Hubbard model, i.e., of the doped Mott insulator, but they may also show up already for small values of , i.e., in the “itinerant part” of the phase diagram.
The paper is organized as follows. The variational ansatz is introduced in Section II, where also the linkedcluster expansion is explained. Section III presents the variational groundstate energy for the normal state (vanishing gap parameter), in comparison to a straightforward perturbative expansion in powers of . In Section IV the formalism is applied to the superconducting state (finite gap parameter). Some details about the minimization procedure are also given. Section V deals with the condensation energy and the delicate balance between kinetic and potentialenergy lowering. Results for the gap parameter are discussed in detail in Section VI, including its “domelike” dependence on particle density as well as an unconventional dependence on . Section VII introduces the superconducting order parameter and shows that it has a similar domelike shape as the gap parameter. The dependence on system size is examined in Section VIII. At half filling superconductivity disappears in the thermodynamic limit, while away from half filling a criterion is found for the characteristic size above which the gap parameter has converged (to a finite value). Irregularities below this characteristic size are attributed to fluctuations in the energy gap between HOMO and LUMO levels. The paper is summarized in Section IX, where also a few problems are listed which could be studied in the future. Some technicalities are discussed in Appendices A and D. In Appendix B an effective model is treated which has a superconducting meanfield ground state with wave symmetry. Appendix C presents results for the standard Gutzwiller ansatz.
Ii Variational approach
ii.1 Hamiltonian
We consider the (repulsive) Hubbard Hamiltonian , with nearestneighbor hopping
(1) 
and onsite repulsion
(2) 
where the operator () creates (annihilates) an electron at site with spin projection or , and . We restrict ourselves to a square lattice with sites and a lattice constant 1. The hopping amplitude is taken as the unit of energy, i.e., we put . In Fourier space we have
(3) 
where the Kronecker symbol is equal to 1 if is a reciprocal lattice vector and 0 otherwise. The tightbinding spectrum is given by
(4) 
The identity
(5) 
will be extensively used later on.
ii.2 Variational ansatz
The (standard) Gutzwiller ansatz reads
(6) 
where measures the number of doubly occupied sites, is a variational parameter and is the ground state of (the filled Fermi sea). To deal with ordering phenomena, such as antiferromagnetism or superconductivity, one uses, instead of , the ground state of a symmetrybreaking meanfield Hamiltonian as the reference state.
The ansatz (6) has been widely adopted in the limit , where double occupancy is completely suppressed and the Hubbard Hamiltonian can be replaced by its large limit, which corresponds to the  model Edegger_07 (); Ogata_08 (). The case of finite with a superconducting reference state has been treated both numerically, using Monte Carlo sampling Giamarchi_91 (), and partly analytically, by a perturbative expansion Kaczmarczyk_13 ().
Unfortunately, the ansatz (6) has its weaknesses. Both the energy and the momentum distribution are at odds with exact results in one dimension, where the Gutzwiller wave function can be analyzed exactly Metzner_89 (). A remedy was proposed already in the early eighties in terms of an additional prefactor which strengthens the correlations between doubly occupied and empty sites (“doublonholon binding”) Kaplan_82 (). This additional term turned out to be important for intermediate to large values of , but not for small . Variational Monte Carlo calculations with this modified wave function yield a superconducting ground state in the intermediate to strong coupling regime and for not too large doping, with clear deviations from BCS behavior Yokoyama_04 (); Yokoyama_13 (). A Jastrow factor producing longrange chargecharge correlations Yokoyama_90 () has also been proposed. It can lead to longrange order in the absence of a symmetrybreaking field Capello_05 (); Kaneko_16 (). Another way of improving the ansatz (6) is to modify the reference state . For instance, instead of using two parameters for a superconducting meanfield state (the gap parameter and the “chemical potential”), one can independently vary the BCS amplitudes and in this way introduce a huge number of variational parameters (of the order of ) Tahara_08 (). It has also been proposed to incorporate a “backflow term” to improve the accuracy and to account for the kinetic exchange for large values of Tocchio_08 ().
Our ansatz reads
(7) 
where is the ground state of the meanfield Hamiltonian (with energy eigenvalue ) and is an additional variational parameter. Eq. (7) differs slightly from the ansatz used in previous (variational Monte Carlo) studies Eichenberger_07 (); Yanagisawa_16 (), where the “kinetic energy” was used in the exponent. The choice of instead of is very convenient for a perturbative evaluation of expectation values, as will become clear below.
ii.3 Linked cluster expansion
When Gutzwiller introduced his wave function Gutzwiller_63 (), he adopted the linked cluster expansion for calculating expectation values. Here we show that the expansion can also be used for our ansatz (7).
The variational parameters plus those defining the meanfield state are determined by minimizing the energy
(8) 
for a given average number of particles
(9) 
where . Eq. (7) can be written as
(10) 
where we have introduced the notation
(11) 
for any operator . Correspondingly, we have
(12) 
and, because of the linkedcluster theorem,
(13) 
where means that only those diagrams have to be taken into account where either both exponentials are connected to the operator or one of the two is connected to and the two exponential operators are connected to each other. We carry the expansion out to second order in for the hopping term and to first order for the interaction . This is justified because the optimized correlation parameter is linear in for and hence the secondorder contribution to the interaction part would be proportional to , i.e., negligible at this order of the expansion.^{1}^{1}1This is strictly valid only in the absence of symmetry breaking, but the argument should remain valid as long as . We find
(14) 
where we have introduced the notation
(15) 
for averages with respect to the meanfield ground state. The average particle number is calculated in the same way, and we obtain
(16) 
If is an eigenstate of all the connected terms vanish and . This is the case for or for an antiferromagnetic reference state, but not for a BCS state, for which Eq. (16) together with the constraint of a fixed yields a nontrivial relation between the variational parameters. For the contributions and also vanish.
For a fixed meanfield state and a fixed parameter , the energy (14) is easily minimized with respect to . Obviously has to be small enough, otherwise the limitation to second order is no longer a good approximation. How small? A simple argument can be given by looking at the problem of two particles on two sites, for which the Gutzwiller ansatz is exact. One readily finds that the secondorder expansion corresponds to the replacement
(17) 
For this approximation is very good, better than , and it is still better than for . This simple estimate agrees with an explicit comparison between a full evaluation of the Gutzwiller ansatz (variational Monte Carlo) and the secondorder expansion for the onedimensional PeierlsHubbard model, showing good agreement for .Jeckelmann_94 () In the present paper we limit ourselves to the region where does not exceed values of the order of 0.6. As will be shown later, this criterion implies that should be lower than about 1.2 for the full ansatz (7), while the Gutzwiller wave function (6) admits values up to 3.
Iii Normal state
We discuss first the symmetric case, , . The various terms of the expansion (14) are readily calculated by Wick decomposition. The zerothorder term is given by
(18) 
where
(19) 
is the Fermi function (which does not depend on for a nondegenerate state, where all levels are either fully occupied or unoccupied), and
(20) 
is the particle density. It is convenient to introduce also the Fermi function of holes,
(21) 
The firstorder contributions are illustrated by diagrams in Fig. 1. The lines represent
(22) 
or
(23) 
where is the singleparticle spectrum measured with respect to the Fermi energy . Thus the parameter acts like a soft energy cutoff, which renders correlation effects strongest close to the Fermi surface. We know already that the contribution , which represents the term , has to vanish. This follows also from the fact that diagram involves the factor . The same is true for the contribution , which therefore also vanishes. To evaluate the contribution , we transform the threefold momentum sum into a single sum over lattice sites using Eq. (5). We obtain
(24) 
where ( or )
(25) 
The secondorder contributions are illustrated in Fig. 2. Diagrams and yield vanishing contributions, while diagram can be written as a sum over lattice sites, by defining the quantities
(26) 
We obtain
(27) 
The minimization of Eq. (14) with respect to yields the correlation energy
(28) 
which is readily evaluated numerically because the individual terms, Eqs. (24) and (27), are simple sums. The result has still to be minimized with respect to the parameter . It turns out that depends weakly on (for small ), but more sensitively on the particle density . Some results for the particular case of are given in Table 1 for a square lattice. Both and are largest at half filling.
1.0  0.59633  0.19805  0.012169 
0.9  0.57690  0.19286  0.011774 
0.8  0.54594  0.18435  0.010787 
0.7  0.51750  0.17618  0.009406 
0.6  0.49427  0.16917  0.007774 
0.5  0.47666  0.16349  0.006014 
It is instructive to compare these results with the exact secondorder term deduced by perturbation theory Metzner_89 (). The latter can be written as
(29) 
The results shown in Table 2 confirm that our variational ansatz (Table 1) reproduces the exact secondorder term to a high precision ( at half filling, for ) and performs much better than the Gutzwiller ansatz ( at half filling and 86% for ). This large improvement of the correlation energy by the parameter also holds for larger values of Otsuka_92 (); Dzierzawa_95 (). We notice that the values for the variational ansatz including are much larger than those for the Gutzwiller ansatz (where ). This happens because with increasing the contribution of the region away from the Fermi surface is reduced and thus the cost in band energy due to the reduction of double occupancy is lowered, and assumes higher values than for .
The relatively large values of the correlation parameter for small values of reflect the fact that small bare interactions do not necessarily imply weak correlations.
1.0  0.14717  0.010220  0.012562 
0.9  0.14639  0.009961  0.012072 
0.8  0.14457  0.009236  0.010995 
0.7  0.14216  0.008139  0.009553 
0.6  0.13944  0.006771  0.007880 
0.5  0.13652  0.005246  0.006092 
Iv Superconducting state
iv.1 Reference state
For superconductivity the reference state is the ground state of the meanfield Hamiltonian
(30) 
where and the gap parameter must have an appropriate symmetry, such as wave or wave. In this paper we restrict ourselves to wave symmetry, i.e.,
(31) 
The parameter is chosen in such a way that the average particle number is equal to a fixed value. For the correlated ground state (7) can be identified with the (true) chemical potential only for and .
is diagonalized by the Bogoliubov transformation
(32) 
if is chosen as
(33) 
where
(34) 
is the excitation spectrum. The meanfield Hamiltonian now reads
(35) 
Its ground state is the vacuum for quasiparticles, . It is then easy to see that
(36) 
Three different functions appear in the Wick decomposition, the momentum distribution functions
(37) 
and the “Gor’kov function”
(38) 
iv.2 Secondorder expansion
We are now prepared for carrying out explicitly the expansion (14) for a superconducting reference state. The contribution of zeroth order is given by
(39) 
where
(40) 
For an order parameter with  or wave symmetry the “average Gor’kov function” vanishes if both the lattice and the boundary conditions have fourfold rotational symmetry. In order to cope with slight deviations from this symmetry for finite system sizes (due to periodicantiperiodic boundary conditions) we retain in the analytical expressions. For the large system sizes considered here the breaking of the fourfold rotational symmetry has very little effect. In recent variational Monte Carlo studies Ido_18 (), with periodicantiperiodic boundary conditions and up to 24, striped phases have been found to be slightly more stable than homogeneous superconductivity. Because breaking of the fourfold rotational symmetry is expected to favor stripes and to weaken wave superconductivity, one may ask whether these results survive for larger system sizes or for symmetric boundary conditions.
The contribution of first order in has three terms
(41) 
where comes from the hopping part of the Hamiltonian and from the interaction. They correspond to the three diagrams of Fig. 1, where a line can represent any of the three functions , and are given explicitly by
(42) 
where
(43) 
with
(44) 
The triple summation in is replaced by a simple summation over lattice sites using Eq. (5) together with the Fourier transform
(45) 
and correspondingly for the other functions. We get
(46) 
We now turn to the secondorder contribution in Eq. (14). The diagrams are grouped according to the three general structures of Fig. 2 but, because of the various possibilities for lines (representing , or ) and Hartree bubbles (density per spin or average Gor’kov function) there are many different specific diagrams, namely 29 of type , 18 of type and 16 of type . Nevertheless, the result can be presented in a relatively compact form, as shown in Appendix A. The numerical evaluation of the various terms requires only simple summations, either in  or in space.
The secondorder expansion of the particle number, Eq. (16) is effectuated in the same way. To deduce the corresponding formulae one simply has to replace by 1 and by 0 in the expression for the energy.
iv.3 Numerical procedure
In the numerical calculations we have considered finite quadratic arrays of size with up to 4000. Thus the number of points in the Brillouin zone is . Periodicantiperiodic boundary conditions have been used, in order to reduce level degeneracies. For a given density and a given system size the particle number is chosen in such a way that there be no ambiguity in the reference state at . For and this is the case for because with this choice the “highest occupied molecular orbital” (HOMO) is completely full and the “lowest unoccupied molecular orbital” (LUMO) is completely empty and there is no degeneracy in the reference state. For a lattice of sites a particle number of 8000 would not lead to a fullshell situation, the closest number satisfying this criterion is . For a finite gap parameter the constraint of a fixed particle number has to be satisfied very accurately, because the energy gained by pairing, the “condensation energy”, is much smaller than the correlation energy. The results presented below have been obtained with a precision of at least for the density .
Four parameters have to be determined by minimizing the energy for a fixed density, namely , , and . To reduce the complexity of the problem, we use the fact that the parameters and interfere weakly. Therefore we determine the optimal value of initially, i.e., for . We have checked that a full variational treatment of all parameters would only slightly increase the stability of the superconducting state. This is illustrated in Table 3, where the full optimization is shown to enhance the gap parameter by about . Correspondingly, the condensation energy increases slightly. In what follows, we will restrict ourselves to the “initial optimization”. Some examples of optimized parameters are given in Table 4 for . The parameter is practically independent in this range, while the correlation parameter is proportional to and the gap parameter varies much more strongly. A more detailed discussion of the dependence of the gap parameter will be given in Section VI.
1.0  0.198053  0.0004247  0.198882  0.0004252 
0.8  0.184349  0.0026698  0.186988  0.0027491 
0.6  0.169172  0.0020630  0.170660  0.0020878 
0.6  0.184349  0.327026  0.0004204 
0.8  0.184349  0.436102  0.0010891 
1.0  0.184349  0.546204  0.0026698 
1.2  0.184349  0.660916  0.0067672 
V Energetics
v.1 Energy gain
After the initial minimization with respect to for , the energy is calculated for a fixed gap parameter and minimized analytically with respect to , while is determined by the constraint of a fixed density. This yields the parameters , and the energy difference
(47) 
Fig. 3 shows this quantity for , and four different system sizes. Negative values of imply that the system is unstable with respect to superconductivity. Both for and for , exhibits two minima, but the two curves differ appreciably from each other and from curves for larger system sizes. By contrast, the results for and are quite similar, with a single minimum at . These finitesize effects will be discussed in more detail in Section VIII.
For comparison, is also presented in Fig. 3 for the Gutzwiller ansatz (). A minimum is again found, but its position is an order of magnitude below that obtained for finite . Correspondingly, the energy gain is nearly two orders of magnitude smaller. Therefore the variational parameter enhances both the correlation energy (as shown in Section II) and the energy gain due to superconductivity. This is an important observation because one could imagine a poor correlation energy to be compensated by an artificially large condensation energy.
v.2 Condensation energy
The minimization of the energy yields the optimized values of the parameters and . The optimized gap parameter ^{2}^{2}2From now on the symbol will be used for the optimized gap parameter. will be detailed in Section VI. Fig. 4 shows the condensation energy, which we define as
(48) 
in agreement with BCS theory Bardeen_57 (); Schrieffer_64 (). first increases when moving away from half filling (i.e., upon doping the parent halffilled system), passes through a maximum for and then decreases. The values differ little when passing from to , except at half filling where decreases with increasing .
In BCS theory the condensation energy is related to the density of states and the gap parameter by the simple formula
(49) 
Our results (solid line in Fig. 4) differ markedly from the BCS prediction (dashed line in Fig. 4). To find out whether the discrepancy can be attributed to the symmetry of the gap function (wave due to a local attraction in BCS theory, wave in the present case), we have calculated the condensation energy for an extended Hubbard model with nearestneighbor attraction using the meanfield ground state of section IV.1. Details are given in Appendix B. The results agree quite well with Eq. (49), which seems to hold approximately for any BCStype theory. But why is the condensation energy so much larger in the present case than what would be predicted by BCS theory? We attribute the difference to the correlation energy, which involves not only the region very close to the Fermi energy, but also band states further away. In fact, the correlation parameter  and therefore the correlation energy  increases with , as can be verified by comparing Tables 1 and 4 (for ). The disparity would be even more pronounced if the fully optimized value of would be used for the superconducting state.
v.3 Conventional or unconventional?
Unconventional superconductivity is not a sharply defined concept. Sometimes it is associated with an unconventional symmetry of the order parameter, and sometimes the emphasis is on properties deviating markedly from BCS predictions or on the nonphonon glue mediating the effective attraction Sigrist_91 (); Stewart_17 (). Superconductivity in the repulsive Hubbard model is unconventional in several respects, in its order parameter (wave close to half filling, wave further away), in the mechanism (no phonons by assumption, maybe exchange of spin fluctuations or no glue at all) and also in deviations from BCS predictions. Here we discuss the question whether pairing is due to a decrease in potential energy, as in BCS theory, or rather due to an unconventional decrease in kinetic energy.
For the reduced BCS Hamiltonian it is easy to convince oneself that the kinetic energy cannot be lowered by pairing. For this model the normal state corresponds to the filled Fermi sea, which has the lowest possible kinetic energy for a given number of particles. Any interaction effect must then lead to an increase of kinetic energy, and superconductivity can only occur if this increase is overcompensated by a decrease in potential energy.
The issue whether the condensation energy is generated by a gain in potential energy, as in BCS theory, or by a gain in kinetic energy has been addressed in the frameworks of spinfluctuation exchange Yanase_05 (), cluster dynamical meanfield theory Maier_04 (); Gull_12 () and variational methods Yokoyama_04 (); Yokoyama_13 (); Kaczmarczyk_13 (); Tocchio_16 (). Quite generally, an unconventional gain in kinetic energy is found for large values of and/or weak doping, while a conventional gain in potential energy is obtained for heavy doping and/or not too large , but the detailed predictions differ somewhat. For instance, different variational wave functions may give different answers for the same values of and for the same density Eichenberger_07 ().
It is straightforward to calculate individually the changes in potential and kinetic energies due to superconductivity within the present approach. The results, shown in Fig. 5 for , are quite surprising, because the kinetic energy is lowered for heavy doping (), while for weak doping () there is a gain in potential energy, contrary to what is typically found in numerical calculations for large values of . There is no contradiction with the arguments given for the reduced BCS Hamiltonian, because the normal state () is correlated and has a kinetic energy exceeding its minimum value.
The corresponding results for the Gutzwiller ansatz () are slightly different, as shown in Appendix C.
Vi Gap parameter
The optimized gap parameter is given in Fig. 6 for as a function of density for various system sizes. It shows a similar domelike shape as the condensation energy (depicted in Fig. 4), although less pronounced. Again the results vary little with system size from to except at half filling where they decrease steadily. Finitesize scaling (discussed in Section VIII) indicates that at half filling the gap parameter vanishes in the thermodynamic limit.
The gap parameter for the Gutzwiller ansatz () is much smaller for than the data for shown in Fig. 6, but the densitydependence is rather similar, as shown in Fig. 12 of Appendix C for and .
An interesting question is how varies with . Adopting the RPA expression for the effective interaction induced by the exchange of spin fluctuations Scalapino_95 (); Kondo_01 () in the small limit, one obtains an attraction proportional to and a BCS behavior for the gap parameter,
(50) 
The same dependence on has been found for the critical temperature (which is expected to be proportional to ) using a renormalization group approach Raghu_10 ().
To see whether the dependence of Eq. (50) also comes out from our variational method, we have calculated the gap parameter for various values of . These are limited to a relatively small region because for too large values of the secondorder expansion breaks down and for too small values of the tiny condensation energy would require a higher precision than what we used in the calculations. Fig. 7 shows results for and . A linear relationship between and fits the data very well, giving
(51) 
with and . By contrast, the functional behavior (50) does not fit these data. Results for are very similar, but because of the larger values of the correlation parameter (as compared to those for ) the range of values has to be reduced. It is interesting to note that is of the order of the interaction strength where a Mott transition would occur in the absence of antiferromagnetic ordering Fratino_17 ().
Vii Order parameter
The physical interpretation of the gap parameter is far from obvious for any variational treatment of a superconductor with strong correlations. Thus it cannot simply be associated with a pseudogap, although this may yield an appealing picture of weakly doped cuprates Anderson_04 (). In fact, the value of depends quite strongly on the choice of the wave function Eichenberger_07 (). By contrast, the order parameter can be sharply defined as an expectation value, and it depends much less on details Baeriswyl_09 (). In “canonical” calculations (fixed particle number) the order parameter is deduced from the longdistance behavior of the pairpair correlation function Vondelft_01 (). In the present “grandcanonical” approach it is defined as the pair amplitude on nearestneighbor sites ,
(52) 
where
(53) 
creates a Cooper pair. The expansion in powers of proceeds in exactly the same way as for the hopping term and we find
(54) 
The zerothorder term is just given by the Gor’kov function, as in BCS theory. The firstorder term reads
(55) 
and is represented by diagram of Fig. 1. The secondorder contributions correspond to the diagrams of Fig. 2 and are given explicitly in Appendix D.
The numerical evaluation of the order parameter is again straightforward, as only simple  and sums have to be calculated. The result for and is shown in Fig. 8. has a maximum at and is to a large extent given by the zerothorder contribution (the Gor’kov function). For additional results for larger system sizes agree with the asymptotic behavior found for the gap parameter, which will be discussed below.
Viii Finitesize effects
To highlight the size dependence of the gap parameter we have plotted vs. for three different densities in Fig. 9. Both for and for , shows a rather irregular behavior for , but then it approaches a finite limiting value for the largest system sizes. By contrast, the data for are smooth for and well described by the size dependence , as evidenced in Fig. 10. This confirms that in the thermodynamic limit the superconducting gap vanishes at half filling. A similar behavior has been observed for the Gutzwiller ansatz (not shown).
Already in 1959 Anderson wondered what happens to a superconducting material if its size shrinks more and more Anderson_59b (). He argued that superconductivity ceases as soon as the characteristic level spacing becomes larger than the energy gap of the bulk system. In the nineties, Anderson’s question was investigated thoroughly, both in spectroscopic experiments on ultrasmall aluminium particles and theoretically using the exact solution of the reduced BCS Hamiltonian Vondelft_01 (). Anderson’s estimate was confirmed, at the same time the critical size was found to signal a crossover rather than a true transition.
We cannot directly apply Anderson’s criterion because the gap function (31) depends on with nodes for . Nonetheless, the gap parameter can be used as an average scale, because it is the root mean square of the gap function,
(56) 
A further difficulty is that, due to the anisotropy of , the level spacing, or more precisely the HOMOLUMO gap (the separation between the lowest unoccupied and the highest occupied singleparticle levels for ), changes irregularly with system size. This is demonstrated in Fig. 11 for a density . By contrast, the HOMOLUMO gap is smooth for and can be given analytically,
(57) 
Apparently, a smooth HOMOLUMO gap implies a smooth superconducting gap, while a chaotic dependence of on system size produces irregularities in .
The horizontal line in Fig. 11 corresponds to the value of the gap parameter for , and its crossing with the HOMOLUMO line gives a critical length of , according to Anderson’s criterion. Clearly, the fluctuations in are larger than for , but even for they remain substantial. Therefore it is not surprising that in the results for (Fig. 9) the irregularities become small only for .
From this analysis we obtain a simple criterion for the size above which the results for the gap parameter are reliable (although not necessarily converged to the thermodynamic limit, cf. Fig. 10. The superconducting gap parameter should be a few times larger than the level spacing . This criterion is expected to be valid for any numerical treatment of the twodimensional Hubbard model on a finite lattice, be it quantum Monte Carlo, variational treatments or dynamical meanfield theory. For of the order of the bandwidth () the typical size of the superconducting gap is 0.1. This is also the typical size of the HOMOLUMO gap for a lattice. It is then understandable that calculations carried out with lattices of this size cannot produce conclusive results for , where the gap parameter is much smaller.
Ix Summary and outlook
In this paper the repulsive Hubbard model has been scrutinized for wave superconductivity, using a refined Gutzwiller ansatz for the ground state. The operator , which simply reduces double occupancy of a reference state (a BCS meanfield state in the present case), was supplemented by a term , which involves the BCS meanfield Hamiltonian . This ansatz is well suited for treating the small region. On the one hand, it admits a linkedcluster expansion in powers of the correlation parameter , as in the case of the standard Gutzwiller ansatz. On the other hand, the energy is greatly improved by the additional term, as evidenced by a comparison with the expansion of the exact groundstate energy. Moreover, this approach has the advantage that large system sizes (millions of sites) can be treated to a very high precision. The drawback is that the method is limited to a relatively small interval of coupling strengths, , where (for lattices which are not larger than a few million sites) and (to keep smaller than about 0.6).
The following main results were obtained.

Superconductivity with wave symmetry exists away from half filling () with a maximum stability for .

Both the gap parameter and the order parameter are an order of magnitude larger than in the case of the Gutzwiller ansatz.

Correspondingly, the condensation energy also has a maximum for and is larger by nearly two orders of magnitude than in the Gutzwiller ansatz.

The pairing is due to a lowering of potential energy for the “underdoped” case and due to the lowering of kinetic energy for the “overdoped” case, in contrast to what is generically reported for large values of .

Our findings resolve a discrepancy between variational Monte Carlo studies, which did not find signatures of superconductivity below some critical value of (of the order of 4) and perturbative treatments such as the functional renormalization group or the fluctuationexchange approximation, which do find a superconducting instability for small values of . The main reason for the failure of variational Monte Carlo calculations is that the system sizes that can be treated are not large enough for the small gap parameters found for .

In the special case of a halffilled band () and vanish as the system size tends to infinity.

For other fillings a limiting finite value for and is obtained if exceeds a certain size. The (Andersontype) criterion is that the bulk gap ( for ) should be a few times larger than the gap between HOMO and LUMO levels.

For smaller lengths fluctuations in lead to irregularities in the results.

The dependence of the gap (for ) is consistent with a power law, , where and is of the order of the critical value for the Mott transition at half filling.
Four variational parameters have been inserted into the trial wave function, the gap parameter , the “chemical potential” , the correlation parameter , and , the inverse of a soft energy cutoff. The parameter can also be interpreted as a characteristic imaginary time. The exact secondorder contribution, Eq. (29), is indeed an integral over the imaginary time . The corresponding variational term, Eq. (24), looks very similar, but the integral is replaced by the integrand at “time” . Since this characteristic time strengthens superconductivity, one may wonder whether its role is to introduce retardation in some effective interaction. It would be worthwhile to study this question thoroughly.
The power law contradicts simple theories based on the exchange of spin fluctuations. It would be great if it would be possible to prove or disprove this unusual behavior, for instance using the functional renormalization group in the ordered phase.
Other questions could also be addressed, such as wave superconductivity (which is expected to dominate for small densities ), the competition between superconductivity and antiferromagnetism (which is expected to be weak in view of the vanishing superconducting order parameter at half filling, where antiferromagnetic ordering is strongest) or the effect of nextnearestneighbor hopping (which could bring back superconductivity at half filling).
Acknowledgment
I have profited from a close collaboration with Mikheil Menteshashvili at an early stage of this work and also from useful discussions with Florian Gebhard and Jörg Bünemann.
Appendix A Secondorder terms for the energy
The secondorder contributions of Eq. (14) can be grouped according to the diagrams of Fig. 2, and we write
(58) 
Proceeding as in Section III, we introduce the quantities
(59) 
where
(60) 
All three diagrams of Fig. 2 contribute if the gap parameter is finite and we write, correspondingly,
(61) 
We find
(62) 
where and are given by Eqs. (43) and (44), respectively, and , and so on. , , and are simple sums, but also in and there are no true double sums because the  and dependent terms can be handled independently.
For periodic boundary conditions, where is even under a reflection by the diagonal and is odd, the above expressions are simplified. We have chosen periodicantiperiodic boundary conditions, for which this symmetry does not hold as long as remains finite. Therefore we have used the full expressions (62) in the computations.
Appendix B BCS pairing for nearestneighbor attraction
The extended Hubbard model with repulsive onsite and attractive nearestneighbor interactions, defined by the Hamiltonian
(63) 
(, , ), is unstable with respect to wave pairing at the BCS meanfield level. Using the ansatz of Section LABEL:31, one readily finds the expectation value
(64) 
where , and are defined, respectively, by Eqs. (37), (38) and (40).
n  V  

0.8  0.5  0.010131  0.8159  
0.8  0.8  0.089255  0.8096  
0.8  1.0  0.176589 