Finite-temperature superconducting correlations of the Hubbard model
aa r X i v : . [ c ond - m a t . s t r- e l ] J u l Finite-temperature superconducting correlations of the Hubbard model
Ehsan Khatami,
1, 2
Richard T. Scalettar, and Rajiv R. P. Singh Department of Physics, University of California, Davis, CA 95616, USA Department of Physics and Astronomy, San Jose State University, San Jose, CA 95192, USA
We utilize numerical linked-cluster expansions (NLCEs) and the determinantal quantum MonteCarlo algorithm to study pairing correlations in the square lattice Hubbard model. To benchmarkthe NLCE, we first locate the finite-temperature phase transition of the attractive model to asuperconducting state away from half filling. We then explore the superconducting properties ofthe repulsive model for the d -wave and extended s -wave pairing symmetries. The pairing structurefactor shows a strong tendency to d -wave pairing and peaks at an interaction strength comparableto the bandwidth. The extended s -wave structure factor and correlation length are larger at highertemperatures but clearly saturate as temperature is lowered, whereas the d -wave counterparts,which start off lower at high temperatures, continue to rise near half filling. This rise is even moredramatic in the d -wave susceptibility. The convergence of NLCEs breaks down as the susceptibilitiesand correlation lengths become large, so we are unable to determine the onset of long-range order.However, our results extend the conclusion, previously restricted to only magnetic and chargecorrelations, that NLCEs offer unique window into pairing in the Hubbard model at strong coupling. PACS numbers: 71.10.Fd, 74.72.-h, 67.85.-d, 05.10.-a
I. INTRODUCTION
Despite several decades of intensive theoreticalresearch, the question of whether a non-local attractioncan dominate in a fermionic Hubbard model with localrepulsive interaction has remained largely unansweredfor parameters relevant to cuprate high-temperaturesuperconductors.
Controlled theoretical approachesconfirm this possibility, however, only when the strengthof the local repulsion is much smaller than the hoppingamplitude of fermions on a square lattice. Numerical methods provide important data forstrongly-correlated quantum Hamiltonians, and, inparticular, for phenomena like superconductivity,magnetism, and Mott metal-insulator transitions.Although many developments have made theseapproaches increasingly powerful over the last decade,significant limitations remain, especially for fermions.The density matrix renormalization group, andrelated techniques, function best in one dimension.Diagrammatic quantum Monte Carlo techniques are restricted to weak-coupling regimes. Determinantquantum Monte Carlo (DQMC), and clusterextensions of the dynamic mean-field theory arelimited to real space or momentum space clusters of tensto hundred of sites. Moreover, the “sign problem” remains an unsolved problem which limits accessibletemperatures unless special symmetries prevail.These limitations emphasize the need for continuedalgorithm development. Recently developed numericallinked-cluster expansions (NLCEs) are especiallypromising as an approach to access strong couplingregimes, which are inaccessible to QMC methods, asa consequence both of the sign problem and of largeand even diverging statistical fluctuations. For instance,analysis of magnetic correlations and Mott phases intrapped atoms on optical lattices, where strong coupling is present at the cloud edge, would not havebeen possible without NLCEs.A natural next step is the application of NLCEs tosuperconductivity. In this Rapid Communication, weshow this method can be developed and successfully usedto study the pairing correlations in the square latticeHubbard model, H = − t X h ij i σ c † iσ c jσ + U X i n i ↑ n i ↓ − µ X iσ n iσ , (1)where c iσ ( c † iσ ) annihilates (creates) a fermion with spin σ on site i , n iσ = c † iσ c iσ is the number operator, U is theonsite Coulomb interaction, and t is the near neighborhopping integral. We set k B = 1, and t = 1 as the unitof energy throughout the paper.We complement our NLCE results with those obtainedfrom (numerically unbiased) DQMC simulations on alarge lattice. We find excellent agreement between thetwo in parameter regions accessible to both, and showthat the lowest temperatures achievable in the NLCE aresimilar to, or often lower than, those of the DQMC. Forthe attractive model ( U < U equalto the bandwidth, the s -wave pairing structure factorof the attractive model away from half filling showsdivergent behavior at low temperatures, and points toa finite transition temperature that is consistent withfindings of previous large-scale DQMC studies. Forthe repulsive model, we consider several values of U and doping and study pairing in the nonlocal channelsof extended s -wave ( s ∗ -wave) and d -wave. While thestructure factor for the former symmetry tends tosaturate at increasingly high temperatures as the dopingis increased, for the latter symmetry, no such tendencyis observed. We examine results at 10% doping moreclosely and find that the low-temperature structure factoris maximum around U = 8. On the other hand, thepair-field susceptibility, while larger for smaller valuesof U in the intermediate temperature region, shows asharp upturn at the lowest accessible temperatures forthe largest interactions considered. II. NUMERICAL METHODS
In NLCEs, an extensive property of the lattice model,when normalized to the number of sites, is expressed inthe thermodynamic limit in terms of contributions fromfinite clusters of various sizes and topologies that can beembedded in the lattice. Thus, NLCEs use the same basisas the high-temperature expansions (HTEs). However,the calculation of the extensive quantities at the levelof individual clusters is left to an exact numericalmethod, such as exact diagonalization, as opposed to aperturbative expansion in terms of inverse temperature inthe HTEs. A typical expansion involves clusters up to acertain size that are chosen according to a self-consistentcriterion (see below). Despite the lack of an explicitsmall parameter, having a finite number of clusters in theseries inevitably leads to the loss of convergence below acertain temperature, where the correlations in the systemextend beyond a length of the order of the largest sizesconsidered. However, the exact treatment of clustersleads to convergence temperatures that are lower thanthose of HTE with a comparable number of orders.Similar to the analytic Pad´e approximations usedextensively in HTEs, here we take advantage oftwo numerical resummation techniques to improve theconvergence of our series at low temperatures. Weuse the Euler algorithm to resum the last 4-6 termsof the series or the Wynn algorithm with 3 and 4cycles of improvement (details of these techniques canbe found in Ref. 19). We then take the average of fourvalues, the last two orders after the Euler and the lasttwo orders after the Wynn transformations, as our bestestimate. To quantify our confidence in the accuracy ofthe resummed results, we define a “confidence region”around this average where all the values that contributeto the average fall. Thus, the errorbars in our figuressimply mark the boundaries of this region and shouldnot be confused with statistical errorbars.We study the superconducting properties of the modelat several values of the interaction strength and on afine grid of temperature and chemical potential. Thelatter allows us to study the calculated quantities atconstant electronic densities after numerical conversion.As with previous studies of Hubbard models using theNLCE, we employ the site expansion in which the order to which each cluster belongs is determined bythe number of sites it has. In order l , we considerall the open boundary clusters of various shapes andtopologies on the square lattice that have l sites, anduse exact diagonalization to solve for their properties.For the pairing correlations, the Hamiltonian matricesare block-diagonalized in each particle number sector.So, we are able to carry out the expansion to the ninthorder. For the pairing susceptibility, on the other hand,we can only carry out the expansion to the seventhorder since not only particle number is not conservedduring the time-dependent measurements (see Eq. 6), butalso the majority of the computational time is spent onobtaining the off-diagonal expectation values, which, likethe diagonalization, scales like O ( N ). DQMC simulations are performed on a 10 ×
10 lattice,which is large enough to have only small finite size effectsat the temperatures studied here. Results representaverages of at least 8 independent runs with 10,000sweeps each. To fix the density, n , away from half fillingat each temperature and U value, the chemical potentialneeds to be tuned starting from an estimate providedby the NLCE. Therefore, we repeat the calculationsfor several values of µ to achieve an accuracy of about0.01% for the density. For the structure factor, weextrapolate our results to the continuous imaginary timelimit using the outcome of two separate simulations witha discretization of the inverse temperature β = L ∆ τ corresponding to ∆ τ = 1 /
16 and 1/12. In the caseof the susceptibility, we choose an even smaller ∆ τ =1 /
50, in order to perform the imaginary time integrationaccurately. This value leads to Trotter errors that arenegligible in comparison to the statistical ones.One of the quantities we calculate is the equal-timepairing structure factor, S α ( q ) = X r e i q · r P α ( r ) , (2)where P α ( r ij ) = h ∆ α † i (0)∆ αj (0) + ∆ αi (0)∆ α † j (0) i (3)is the equal-time pair-pair correlation function. Here, thepairing operator for the symmetry α is defined as∆ αi ( τ ) = 12 X j f αij e τH ( c i ↑ c j ↓ − c i ↓ c j ↑ ) e − τH . (4)We consider three pairing symmetries in this study;(local) s -wave, d -wave, and s ∗ -wave. For the s -wavesymmetry, f sij = δ ij . In the case of s ∗ -wave, f s ∗ ij is+1 if i and j are nearest neighbors and j > i (to avoiddouble counting) and zero otherwise. f dij for the d -wavesymmetry is the same as f s ∗ ij except it takes the value − i and j is along the y axis. Here,we consider only the uniform pairing structure factor, S α ( q = 0). The correlated structure factors, S α corr , isobtained by first subtracting off the uncorrelated partsof the expressions in Eq. 3. T S c o rrs ( q = ) NLCEDQMC T ξ s c o rr T P c o rrs ( r ) r=(0,0)r=(1,0)r=(1,1)r=(2,0) T / S c o rrs ( q = ) U= -8, n=0.85 T c =0.11 (a) (b)(c) FIG. 1. (a) Temperature dependence of the s -wave pairingstructure factor at n = 0 .
85 for the attractive model with U = −
8. The line is from NLCE and symbols are from DQMC.The inset shows the inverse of the same function vs T anda low-temperature fit to A exp( B/ √ T − T c ) with T c = 0 . T . In the main panel of (a) and in(b), bare NLCE results before resummations for the last twoorders, 8 th and 9 th , are shown as thin dotted and dashed lines,respectively. Having the uniform structure factor ( S α or S α corr ),the corresponding correlation length, ξ , can also becalculated using, e.g.,( ξ α corr ) = 12 dS α corr ( q = 0) X i | r i | P α corr ( r i ) , (5)where d = 2 is the dimension.The other quantity of interest for superconductivity isthe uniform pairing susceptibility, which is defined as χ α = 1 N Z β dτ h O α ( τ ) O α † (0) i , (6)where O α ( τ ) = P i ∆ αi ( τ ). III. RESULTS
We start with the attractive Hubbard model, forwhich we know there exists a finite-temperatureKosterlitz-Thouless (KT) phase transition to an s -wavesuperconducting state away from half filling. InFig. 1(a), we show the correlated part of the s -wavepairing structure factor from the NLCE for U = − n = 0 .
85, where the superconducting transitiontemperature is expected to be maximal. Results arein excellent agreement with the corresponding DQMCresults, plotted as empty circles in that figure. As can beinferred from previous DQMC simulations with a smaller U , finite-size effects in DQMC will not play a role hereat temperatures as low as T = 0 .
25. Whereas the raw S d ( s * ) c o rr s*-waved-wave T S d ( s * ) c o rr T T ξ d ( s * ) T ξ d ( s * ) T ξ d ( s * ) T ξ d ( s * ) (a) n=1.00 (b) n=0.95 (c) n=0.90 (d) n=0.85U=8 FIG. 2. The uniform d -wave and extended s -wave pairingstructure factors for the repulsive model with U = 8 atdensities 1.00, 0.95, 0.90, and 0.85 vs temperature. Linesare from the NLCE and symbols are from the DQMC. BareNLCE results before resummations for the 8 th and 9 th orders,are shown as thin dotted and dashed lines, respectively.The insets show the correlation lengths from the NLCE vstemperature for each case. NLCE results (before resummations) converge only to T ∼ .
4, the averaged value after resummations suggesta divergent behavior for S s corr at lower temperatures.They lead us to a regime where we can take advantageof extrapolations in temperature in order to obtain anestimate for the critical temperature. We find that afit to the KT form [see the inset of Fig. 1(a)], leadsto T c ∼ .
11, which is in good agreement with resultsof past DQMC simulations. The correlation length,which shows an exponential growth, is also plotted inFig. 1(b). Its behavior is consistent with the trend seenin Fig 1(c) for the pairing correlations growing faster atlonger length scales as the temperature is decreased.We now turn our focus to the main subject of thisstudy; pairing in the repulsive Hubbard model. Weknow that if a similar finite-temperature transition toa superconducting phase takes place in the latter model,the pairing symmetry has to be nonlocal because of theonsite Coulomb repulsion. Therefore, in this case, weonly explore the d -wave and the s ∗ -wave symmetries. Wealso expect the corresponding temperature scales to bemuch smaller than those for the attractive model sincewe are looking for attraction in a repulsive model.In Fig. 2, we show the correlated part of the uniformstructure factor for the two pairing symmetries when U = 8 and at various average densities. At half filling,the series converges to a low enough temperature to makeclear that S α corr eventually saturates as we decrease thetemperature. In the absence of the ‘sign problem’ atthis filling, DQMC can easily access lower temperatures.We see in Fig. 2(a) that, while agreeing excellentlywith NLCE at high temperatures, results from DQMCsimulations confirm the saturation at lower T . As wemove away from half filling into the hole-doped region T S c o rr d U=4U=6U=8U=12 T / S d c o rr U S c o rr d n=0.90 T=0.43U=8 FIG. 3. Temperature dependence of the d -wave pairingstructure factor at n = 0 .
90 for U = 4, 6, 8, and 12. Symbolsare the DQMC results. Top inset shows the same structurefactor vs U at a fixed temperature T = 0 .
43. The bottominset shows the inverse of the structure factor vs T , alongwith a fit to the function A exp( B/ √ T − T c ) for T < . T c = 0 . ( n < . s ∗ -wave structure factor is seen to take place athigher temperatures whereas the d -wave structure factorcontinues to grow at the lowest temperatures accessibleto us, although its over all values decrease as we increasethe doping. Hence, if there is an instability to pairingaway from half filling in this model, it would be in the d -wave and not the s ∗ -wave channel. Interestingly, atsmall dopings near half filling, NLCE results are morereliable at generally lower temperatures than those of theDQMC because of the restrictions imposed by a severesign problem in this region [see Fig. 2(b)]. Nevertheless,results from the two methods match within the errorbarsat the available temperatures for all the dopings studied.The favorability of d -wave over s ∗ -wave pairing isalso evidenced by the behavior of the correspondingcorrelation lengths, shown in the insets of Fig. 2.For example, even though the low-temperature s ∗ -wavestructure factor is larger than the d -wave one away fromhalf filling, its correlation length clearly saturates whilethat of the d -wave keeps rising and becomes larger. Thelatter can explain the higher convergence temperature of S s ∗ corr in comparison to S d corr in Fig. 2(b).Focusing on d -wave pairing at a moderately dopedsystem with n = 0 .
9, we find that at temperatures belowone, the structure factor is largest at U ∼
8, which isequal to the non-interacting bandwidth. This can beseen in Fig. 3, where we show S d corr vs temperature for U = 4, 6, 8, and 12. For U = 4, the DQMC resultsare available at lower temperatures than the NLCE andshow a relatively slow increase of this quantity as thetemperature is decreased. In the top inset of Fig. 3,we see that the structure factor at T = 0 .
43 quicklyrises as U increases from 4, reaches a maximum at U = 8, and then slowly decreases. Beyond U = 12, T χ d U=4U=6U=8U=12 T / χ d n=0.90 FIG. 4. The d -wave pairing susceptibility at n = 0 .
90 vstemperature for U = 4, 6, 8, and 12. Bare NLCE resultsbefore resummations for the last two orders, 6 th and 7 th , areshown as thin dotted and dashed lines, respectively. Symbolsare the DQMC results. The inset shows the inverse ofthe susceptibility vs temperature for the same values of theinteraction strength. we expect this quantity to scale as 1 /U as, in thestrong-coupling regime, the only relevant energy scalewill be the exchange interaction of the correspondinglow-energy t − J model, J = 4 t /U . The bottominset in Fig. 3 shows the inverse of the structure factorat U = 8. Unfortunately, we are not close enoughto a transition temperature to be able to make anyquantitative statement about its value. However, thebest estimate from the DCA for a close value of theinteraction ( U = 7), puts T c around 0.05, which isconsistent with a KT fit to our results for T < . χ d vs temperature at n = 0 . U values andfor larger U values when the temperature is not too low.This includes the susceptibility at U = 4 and n = 0 . (not shown). In all cases, there is a rapid increase in thesusceptibility at low temperatures. However, more termsare needed for the susceptibility to capture the sharprise at low temperatures, and to determine how T c maydepend on U . In future, it would be important to extendthe results for the susceptibility to higher orders and alsocalculate pairing susceptibilities at non-zero momenta.In summary, we have employed two unbiased methods,the NLCE and the DQMC to study finite-temperaturesuperconducting properties of the square lattice Fermi-Hubbard model. To benchmark our NLCE approach, wefirst explore the s -wave pairing in the attractive modelaway from half filling. By fitting our low-temperaturepairing structure factor to known forms, we obtain a T c that is consistent with the best estimate from large-scale QMC simulations. We then investigate the nonlocal s ∗ -wave and d -wave pairing instabilities in the repulsivemodel at various dopings and for several interactionstrengths. We find that the d -wave symmetry has thetendency to be dominant at low temperatures and thatits structure factor has a maximum at U ∼
8. We alsocalculate the pairing susceptibility, which shows a similardivergent behavior in the d -wave channel and a sharpupturn at low temperatures for large interactions.An important potential application of the resultsdescribed here is to ongoing emulation of modelHamiltonians which describe fermionic atoms in opticallattices. NLCEs allow the rapid evaluation of physicalproperties on a dense mesh of Hamiltonian parameters,a requirement for accurate modeling of optical latticeexperiments where the confining potential leadsto spatially varying chemical potential, interactionstrength, and hopping matrix elements. Here we have shown the potential importance of NLCEs as a tool toanalyze pairing in these systems.
ACKNOWLEDGMENTS
This work was supported by the Department ofEnergy under DE-NA0001842-0 (EK and RTS), andby the National Science Foundation (NSF) grantnumber DMR-1306048 (EK and RRPS). This workused the Extreme Science and Engineering DiscoveryEnvironment (XSEDE) under project number TG-DMR130143, which is supported by NSF grant numberOCI-1053575. D. Scalapino,
Does the Hubbard Model Have the RightStuff? in Proceedings of the International School ofPhysics , edited by R. A. Broglia and J. R. Schrieffer(North-Holland, New York, 1994). E. Dagotto, Rev. Mod. Phys. , 763 (1994). S. Kivelson, I. Bindloss, E. Fradkin, V. Oganesyan,J. Tranquada, A. Kapitulnik, and C. Howard,Rev. Mod. Phys. , 1201 (2003). N. Lee, P.A. amd Nagaosa and X. Wen,Rev. Mod. Phys. , 17 (2006). Handbook of High-Temperature Superconductivity , editedby J. Robert Schrieffer and James S. Brooks (Springer,New York, 2007). D. Scalapino, Rev. Mod. Phys. , 1383 (2012). S. Raghu, S. A. Kivelson, and D. J. Scalapino,Phys. Rev. B , 224505 (2010). S. R. White, Phys. Rev. Lett. , 2863 (1992). U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005). N. V. Prokof’ev and B. V. Svistunov,Phys. Rev. Lett. , 2514 (1998). E. Kozik, K. V. Houcke, E. Gull, L. Pollet,N. Prokof’ev, B. Svistunov, and M. Troyer,EPL (Europhysics Letters) , 10004 (2010). S. R. White, D. J. Scalapino, R. L. Sugar, E. Y.Loh, J. E. Gubernatis, and R. T. Scalettar,Phys. Rev. B , 506 (1989). C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell,and R. T. Scalettar, Phys. Rev. B , 075116 (2009). T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler,Rev. Mod. Phys. , 1027 (2005). G. Kotliar, S. Y. Savrasov, G. P´alsson, and G. Biroli,Phys. Rev. Lett. , 186401 (2001). E. Y. Loh, J. E. Gubernatis, R. T. Scalettar,S. R. White, D. J. Scalapino, and R. L. Sugar,Phys. Rev. B , 9301 (1990). V. I. Iglovikov, E. Khatami, and R. T. Scalettar,Phys. Rev. B , 045110 (2015). M. Rigol, T. Bryant, and R. R. P. Singh,Phys. Rev. Lett. , 187202 (2006). M. Rigol, T. Bryant, and R. R. P. Singh,Phys. Rev. E , 061118 (2007). M. Rigol, T. Bryant, and R. R. P. Singh,Phys. Rev. E , 061119 (2007). B. Tang, E. Khatami, and M. Rigol,Comput. Phys. Commun. , 557 (2013). R. Hart, P. Duarte, T. Yang, X. Liu, T. Paiva, E. Khatami,R. Scalettar, N. Trivedi, D. Huse, and R. Hulet,Nature (London) , 211 (2015). P. Duarte, R. Hart, T. Yang, X. Liu, T. Paiva,E. Khatami, R. Scalettar, N. Trivedi, and R. Hulet,Phys. Rev. Lett. , 070403 (2015). R. T. Scalettar, E. Y. Loh, J. E. Gubernatis, A. Moreo,S. R. White, D. J. Scalapino, R. L. Sugar, and E. Dagotto,Phys. Rev. Lett. , 1407 (1989). A. Moreo and D. J. Scalapino,Phys. Rev. Lett. , 946 (1991). T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H.Denteneer, Phys. Rev. B , 184501 (2004). F. Assaad, W. Hanke, and D. Scalapino,Phys. Rev. B. , 4327 (1994). P. Staar, T. Maier, and T. C. Schulthess,Phys. Rev. B , 195133 (2014). H. W. Press, B. P. Flannery, S. A. Teukolsky, and W. T.Vetterling,
Numerical Recipes in Fortran (CambridgeUniversity Press, Cambridge, England, 1999) Chap. 5.1. A. J. Guttmann,
Phase Transitions and CriticalPhenomena , edited by C. Domb and J. Lebowitz,Vol. 13 (Academic Press, London, 1989). E. Khatami and M. Rigol,Phys. Rev. A , 053611 (2011). E. Khatami and M. Rigol,Phys. Rev. A , 023633 (2012). B. Tang, T. Paiva, E. Khatami, and M. Rigol,Phys. Rev. Lett. , 205301 (2012). In the ninth order alone, we need to distinguish all 1285clusters that may have the same topology (the sameHamiltonian matrix), but are not related by point groupsymmetry. The largest matrices to be diagonalized have alinear size of N = 15876. In the seventh order, there are108 symmetrically distinct clusters. The largest matricesto be diagonalized here have a linear size of N = 3432. For example, for the s-wave channel, we have P s corr ( r ij ) = P s ( r ij ) − h c † iσ c jσ i − δ ij h − h c † iσ c jσ i i . S. R. White, D. J. Scalapino, R. L. Sugar, N. E. Bickers,and R. T. Scalettar, Phys. Rev. B , 839 (1989). A. G. Truscott, K. E. Strecker, W. I. McAlexander, G. B.Partridge, and R. G. Hulet, Science , 2570 (2001). M. K¨ohl, H. Moritz, T. St¨oferle, K. G¨unter, andT. Esslinger, Phys. Rev. Lett. , 080403 (2005). U. Schneider, L. Hackermller, S. Will, T. Best, I. Bloch,T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science , 1520 (2008). R. J¨ordens, N. Strohmaier, K. G¨unter, H. Moritz, andT. Esslinger, Nature , 204 (2008). C. J. M. Mathy, D. A. Huse, and R. G. Hulet,Phys. Rev. A86