Superconductivity from repulsion: Variational results for the 2D Hubbard model in the limit of weak interaction

Superconductivity from repulsion: Variational results for the Hubbard model in the limit of weak interaction

Dionys Baeriswyl Department of Physics, University of Fribourg,
Chemin du Musée 3, CH-1700 Fribourg, Switzerland
International Institute of Physics, 59078-400 Natal-RN, Brazil

The two-dimensional 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 mean-field 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 dome-like 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 linked-cluster 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 power-law increase as a function of . Finite-size 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 HOMO-LUMO gap) fluctuates as a function of for the tight-binding 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 narrow-band materials. Later on the model was used for describing the Mott metal-insulator 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 short-range 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 two-dimensional Hubbard model, for positive values of the on-site 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 half-filled 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 Mean-Field 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 two-dimensional 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 mean-field 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 two-dimensional 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 nearest-neighbor sites. Nevertheless, the results are clear-cut, 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 RPA-type theories for spin-fluctuation 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 non-BCS 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 linked-cluster expansion is explained. Section III presents the variational ground-state 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 potential-energy lowering. Results for the gap parameter are discussed in detail in Section VI, including its “dome-like” 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 dome-like 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 mean-field 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 nearest-neighbor hopping


and on-site repulsion


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


where the Kronecker symbol is equal to 1 if is a reciprocal lattice vector and 0 otherwise. The tight-binding spectrum is given by


The identity


will be extensively used later on.

ii.2 Variational ansatz

The (standard) Gutzwiller ansatz reads


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 symmetry-breaking mean-field 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 (“doublon-holon 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 long-range charge-charge correlations Yokoyama_90 () has also been proposed. It can lead to long-range order in the absence of a symmetry-breaking 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 mean-field 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


where is the ground state of the mean-field 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 mean-field state are determined by minimizing the energy


for a given average number of particles


where . Eq. (7) can be written as


where we have introduced the notation


for any operator . Correspondingly, we have


and, because of the linked-cluster theorem,


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 second-order contribution to the interaction part would be proportional to , i.e., negligible at this order of the expansion.111This is strictly valid only in the absence of symmetry breaking, but the argument should remain valid as long as . We find


where we have introduced the notation


for averages with respect to the mean-field ground state. The average particle number is calculated in the same way, and we obtain


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 non-trivial relation between the variational parameters. For the contributions and also vanish.

For a fixed mean-field 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 second-order expansion corresponds to the replacement


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 second-order expansion for the one-dimensional Peierls-Hubbard 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 zeroth-order term is given by




is the Fermi function (which does not depend on for a non-degenerate state, where all levels are either fully occupied or unoccupied), and


is the particle density. It is convenient to introduce also the Fermi function of holes,

Figure 1: First-order diagrams for the expansion (14). Dots represent two-particle vertices, squares single-particle vertices.

The first-order contributions are illustrated by diagrams in Fig. 1. The lines represent




where is the single-particle spectrum measured with respect to the Fermi energy . Thus the parameter acts like a soft energy cut-off, 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


where ( or )

Figure 2: Second-order diagrams with symbols as in Fig. 1.

The second-order 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


We obtain


The minimization of Eq. (14) with respect to yields the correlation energy


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
Table 1: Variational parameters , and correlation energy per site for , and various particle densities .

It is instructive to compare these results with the exact second-order term deduced by perturbation theory Metzner_89 (). The latter can be written as


The results shown in Table 2 confirm that our variational ansatz (Table 1) reproduces the exact second-order 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
Table 2: Variational parameter and correlation energy for the Gutzwiller ansatz, in comparison with the correlation energy obtained by second-order perturbation theory (, ).

Iv Superconducting state

iv.1 Reference state

For superconductivity the reference state is the ground state of the mean-field Hamiltonian


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.,


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


if is chosen as




is the excitation spectrum. The mean-field Hamiltonian now reads


Its ground state is the vacuum for quasi-particles, . It is then easy to see that


Three different functions appear in the Wick decomposition, the momentum distribution functions


and the “Gor’kov function”


iv.2 Second-order 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




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 periodic-antiperiodic 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 periodic-antiperiodic 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


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






The triple summation in is replaced by a simple summation over lattice sites using Eq. (5) together with the Fourier transform


and correspondingly for the other functions. We get


We now turn to the second-order 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 second-order 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 . Periodic-antiperiodic 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 full-shell 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
Table 3: Gap parameters in the cases of full () and “initial” optimization , for , and three different densities .
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
Table 4: Variational parameters for , and several values of .

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


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 finite-size effects will be discussed in more detail in Section VIII.

Figure 3: Energy gain due to superconductivity as a function of the gap parameter for and . Symbols correspond to different system sizes, (circles), (squares), (down-pointing triangles) and (up-pointing triangles). The dashed line represents results obtained with the Gutzwiller ansatz (for .)

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 222From 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


in agreement with BCS theory Bardeen_57 (); Schrieffer_64 (). first increases when moving away from half filling (i.e., upon doping the parent half-filled 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 .

Figure 4: Condensation energy for and different system sizes. Down-pointing triangles: , squares: , circles and solid line: , up-pointing triangle: . The dashed line represents the BCS prediction, Eq. (49).

In BCS theory the condensation energy is related to the density of states and the gap parameter by the simple formula


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 nearest-neighbor attraction using the mean-field 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 BCS-type 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 non-phonon 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 spin-fluctuation exchange Yanase_05 (), cluster dynamical mean-field 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.

Figure 5: Changes in kinetic (up-pointing triangles), potential (down-pointing triangles) and total energies (circles) for and .

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 dome-like 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. Finite-size scaling (discussed in Section VIII) indicates that at half filling the gap parameter vanishes in the thermodynamic limit.

Figure 6: Gap parameter as a function of electron density for and different system sizes. Down-pointing triangles: , squares: , circles and dashed line: , up-pointing triangle: .

The gap parameter for the Gutzwiller ansatz () is much smaller for than the data for shown in Fig. 6, but the density-dependence 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,


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 second-order 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


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 ().

Figure 7: Gap parameter as a function of for and . The solid line is a linear fit through the data points.

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 long-distance behavior of the pair-pair correlation function Vondelft_01 (). In the present “grand-canonical” approach it is defined as the pair amplitude on nearest-neighbor sites ,




creates a Cooper pair. The expansion in powers of proceeds in exactly the same way as for the hopping term and we find


The zeroth-order term is just given by the Gor’kov function, as in BCS theory. The first-order term reads


and is represented by diagram of Fig. 1. The second-order contributions correspond to the diagrams of Fig. 2 and are given explicitly in Appendix D.

Figure 8: Order parameter for and as a function of the density . Circles and the dashed line represent the second-order expansion, squares the zeroth-order contribution.

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 zeroth-order 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 Finite-size 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).

Figure 9: Variation of the gap parameter with system size (). Circles and solid line: , squares: , triangles: .
Figure 10: Size dependence of the gap parameter for at half filling (). The dashed line corresponds to the HOMO-LUMO gap , given by Eq. (57). The solid line is a linear fit through those data points for which exceeds .

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,


A further difficulty is that, due to the anisotropy of , the level spacing, or more precisely the HOMO-LUMO gap (the separation between the lowest unoccupied and the highest occupied single-particle levels for ), changes irregularly with system size. This is demonstrated in Fig. 11 for a density . By contrast, the HOMO-LUMO gap is smooth for and can be given analytically,


Apparently, a smooth HOMO-LUMO gap implies a smooth superconducting gap, while a chaotic dependence of on system size produces irregularities in .

Figure 11: Finite-size behavior of the HOMO-LUMO gap for . The solid line is a coarse-grained version of the data (circles), while the dashed line stands for the bulk value of the gap parameter for .

The horizontal line in Fig. 11 corresponds to the value of the gap parameter for , and its crossing with the HOMO-LUMO 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 two-dimensional Hubbard model on a finite lattice, be it quantum Monte Carlo, variational treatments or dynamical mean-field 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 HOMO-LUMO 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 mean-field state in the present case), was supplemented by a term , which involves the BCS mean-field Hamiltonian . This ansatz is well suited for treating the small region. On the one hand, it admits a linked-cluster 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 ground-state 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.

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

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

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

  4. 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 .

  5. 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 fluctuation-exchange 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 .

  6. In the special case of a half-filled band () and vanish as the system size tends to infinity.

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

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

  9. 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 cut-off. The parameter can also be interpreted as a characteristic imaginary time. The exact second-order 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 next-nearest-neighbor hopping (which could bring back superconductivity at half filling).


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 Second-order terms for the energy

The second-order contributions of Eq. (14) can be grouped according to the diagrams of Fig. 2, and we write


Proceeding as in Section III, we introduce the quantities




All three diagrams of Fig. 2 contribute if the gap parameter is finite and we write, correspondingly,


We find


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 periodic-antiperiodic 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 nearest-neighbor attraction

The extended Hubbard model with repulsive on-site and attractive nearest-neighbor interactions, defined by the Hamiltonian


(, , ), is unstable with respect to -wave pairing at the BCS mean-field level. Using the ansatz of Section LABEL:31, one readily finds the expectation value


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