Evaluation of High Order Terms for the Hubbard Model in the Strong-coupling Limit

Evaluation of High Order T erms for the Hubbard Model in the Str ong-coupling Limit Eva Kalinowski ∗ Academy of Computer Science , 43-300 Bielsko-Biała, P ola nd Władysław Gluza Silesian University of T echno logy , 44-100 Gli wice, P oland (Dated: 24 th June 2011) The ground-state energy of the Hubba rd mode l on a Bethe lattice with infinite connecti vity at half filling is calculated for the insulating phase . Using K ohn’ s transformation to deri ve an ef fectiv e Hamiltonian for the strong-coupling limit, the resulting class of diag rams is determined. W e dev elop an algorithm for an algebraic e va luation o f the contrib utions of high-order terms and ch eck it by applying i t to the Falicov–Kimball model that is exactly solv able. For the Hubbard model, the ground-state energy is ex actly calculated up to order t 12 / U 11 . The resu lts of the strong-co upling expansion deviate from numerical calculations as quantum Monte Carlo (or density-matrix renormalization-gro up) b y less than 0 . 13 % (0 . 32 % respecti vely) for U > 4 . 76. P A CS numbers: 71.30.+h, 71.10.Fd, 71.27.+a, 02.10.Ox INTRODUCTION T H E Hubb ard mod el 1 captures the essential elements of the complex beh avior of stro ngly-co rrelated fermion ic sys- tems with short-range repulsive interactio n. Particular ly in- teresting is th e exploration of the transition from the par am- agnetic metallic phase to a paramagnetic Mott–Hubb ard in- sulator in the limit of in finite d imensions 2–6 ; there , the inter- acting lattice problem can be mapp ed onto effective single- impurity models an d so lved within the framew ork of the dy- namical mean-field theory ( DMFT). Since this ph enomen on was discussed controversly 7 –9 , hig h accuracy in determining the g round -state energy and doub le occupancy per lattice site near the tran sition region is n ecessary for resolvin g doubts as to the nature of the transition, for minimizing quantitativ e un- certainties in the p hase diagr am an d establishing an essential benchm ark for o ther, in par ticular numer ical meth ods. T here were very many attempts to study the mo del in the strong- coupling limit (cf., e. g., r efs. 10–12). Ho we ver , it a ppears to be rather difficult to g o beyond the lowest ord ers. Th erefore, we developed a com puter-algebraic app roach. In this work we pr esent a detailed d escription of a “divide- and-co nquer” algorith m used for an exact calcula tion of all co- efficients in the asympto tic expansion of the gro und-state en- ergy of the Mott insulator including order t 12 / U 11 . Results of such an algorithm up to t 10 / U 9 were alread y q uoted in ref. 13, where an extrapo lation scheme to infin ite order (extrap olated Perturbation The ory , “ePT” ) was introduced. Th is showed ex- cellent agre ement with a quantu m Monte Carlo (QMC) tec h- nique, improved the state of the art by 1 –2 orders of m ag- nitude and lead to a well controlled e valuation of the criti- cal exponen t. Quite recently , our metho d was a pplied to the Bose–Hubba rd model 14,15 . The o utline o f this paper is as fo llows: In sec. I, we ∗ Electroni c address: kalinows@m athematik .uni- marburg.de show h ow the effecti ve Hamiltonian is deri ved following K ohn’ s 16 and Kato’ s 16 and T akahashi’ s 17 treatment of the strong-co upling limit. Then, in sec. II, the class of diag rams defined by the resulting ef fecti ve Hamiltonian is d iscussed for the Bethe lattice with infinite conne ctivity , and the algorith m for the ev aluation of electro nic tran sfer pro cesses o n it is d e- scribed. The co ncept of this algorithm is ideally suited for parallelization that will be done in further work. The results are first gi ven for the Falicov–Kimball model that is exactly solvable and serves as a test of ou r treatment (sec. III. A). Our main resu lt is gi ven in eqs. (19) and (20) in sec. III.B. W e apply our metho d “ePT” ( see ref. 13) a nd co mpare our results with results from DMFT -QMC and DMFT -DDMRG (Dynamica l Density-Matrix Reno rmalization Group) 9 for the groun d-state insulating p hase o f the Hubb ard mo del. Fina lly , flowcharts that present the essential parts of the algo rithm are giv en in the Appendix. I. PER TURBA TION EXP ANSION FOR THE STRONG-COUPLING LIMIT W e in vestigate spin-1 / 2 electro ns o n a lattice represented by the Hubbard model H = T + U D , (1) where T = − t ∑ ( i , j ) , σ c † i σ c j σ is the k inetic energy operator de- scribing electron hops between nearest neighbour sites i and j with the tran sfer amplitu de t , U D = U ∑ i n i ↑ n i ↓ is the interac- tion part inclu ding only local contributions n i σ = c † i σ c i σ . c † i σ and c i σ are the creation and ann ihilation operato rs fo r elec- trons with spin σ = ↓ , ↑ on site i . In the f ollowing, we sketch the calcu lation of an effectiv e Hamiltonian in a stron g coup ling expansio n, as was dev eloped in ref. 1 8 and ref. 17. There it is sho wn how this expansion in 1 / U is done systematically . The aim is th e transfo rmation to n ew particles with an ef fectiv e Hamiltonia n that does n ot change the num ber of doubly occu pied sites. Th e consider a- tions ar e valid for any lattice in a ny dime nsion. The op erator 2 for the kinetic energy T in eq. (1) couples states with a dif fer- ent n umber of doubly occup ied sites. In deriving this effec- ti ve Hamiltonian, a decouplin g can be achieved b y introdu c- ing s uitable linear combinations. Rotating to s uch a ne w basis is p erform ed by a un itary tran sformation U ≡ e S developed by K ohn 16 for the stron g-cou pling limit. Th is transform ation introdu ces ne w particles created by ˜ c † so that c † i σ = e S ( ˜ c ) ˜ c † i σ e − S ( ˜ c ) , (2) and therewith H ( c ) = e S ( ˜ c ) H ( ˜ c ) e − S ( ˜ c ) ≡ ˜ H ( ˜ c ) . (3) The generato r S is construc ted in such a way that the hopping of the n ew p articles does n ot ch ange the effective numb er of doubly occupied sites for the new particles ( ˜ c ), [ ˜ H ( ˜ c ) , D ( ˜ c )] = 0 . (4) This requires S (and therefore ˜ H ( ˜ c ) ) to be an operato r series in 1 / U S ( ˜ c ) = ∞ ∑ i = 1 S i ( ˜ c ) U i , (5) and the un itarity of the tran sformation implies S † = − S . Ob - viously , the wavefunctions can be expressed in ter ms of the new particles, and it f ollows for the eigenenergies h ψ m ( c ) | H ( c ) | ψ m ( c ) i = h ˜ ψ m ( ˜ c ) | ˜ H ( ˜ c ) | ˜ ψ m ( ˜ c ) i = E m . (6) The g round state ˜ ψ 0 ( ˜ c ) of ˜ H ( ˜ c ) at half ban d filling will be determined at the end of this section. The lo w orders in 1 / U ar e con veniantly found by sub stitut- ing the expansion (5) in (3); to second order in 1 / U : ˜ H ( ˜ c ) = T ( ˜ c ) + U D ( ˜ c ) + 1 U [ S 1 ( ˜ c ) , H ( ˜ c )]+ + 1 U 2 [ S 2 ( ˜ c ) , H ( ˜ c )] + 1 2 U 2 [ S 1 ( ˜ c ) , [ S 1 ( ˜ c ) , H ( ˜ c )]] + · · · . (7) From the condition ( 4), the coefficients S n ( ˜ c ) are determined order by order as shown now . The kinetic energy operato r can be separated in three parts, each of which increases or decreases the numb er of double occupancies b y one, o r l eaves it unchan ged, T ( ˜ c ) = T U + T − U + T 0 , (8) where T U = − t ∑ ( i , j ) , σ ˜ n i , − σ ( 1 − ˜ n j , − σ ) ˜ c † i σ ˜ c j σ , T − U = − t ∑ ( i , j ) , σ ˜ n j , − σ ( 1 − ˜ n i , − σ ) ˜ c † i σ ˜ c j σ , T 0 = − t ∑ ( i , j ) , σ ( 1 − ˜ n i , − σ − ˜ n j , − σ + 2 ˜ n i , − σ ˜ n j , − σ ) ˜ c † i σ ˜ c j σ . Because [ T U , D ( ˜ c )] = − T U and [ T − U , D ( ˜ c )] = T − U , (9) (4) is fulfilled when the cross terms T U and T − U are cancelled by [ S ( ˜ c ) , U D ( ˜ c )] in th e lowest order in U − 1 . This is ach iev ed by choosing S 1 ( ˜ c ) = T U − T − U . (10) Inserting in ( 7) and d emanding ( 4), on e ob tains the cond ition for the next order , i. e., S 2 that leads to S 2 ( ˜ c ) = [ T U + T − U , T 0 ] . (11) Follo wing this procedu re, o ne de termines S ( ˜ c ) ord er by or- der . Since S does n ot create or annih ilate bare particles, the vacuum state is eq ual for both old and new p articles. The lowest order term s of the resulting 1 / U -exp ansion of the Hamiltonian ˜ H ( ˜ c ) = U D ( ˜ c ) + ∑ ∞ i = 0 U − i ˜ h i are (here and in the following, T = T ( ˜ c ) ) ˜ h 0 = ∑ j P j T P j , (12a) ˜ h 1 = ∑ j P j T S 1 j T P j , (12b) ˜ h 2 = ∑ j P j T S 1 j T S 1 j T P j − − 1 2 ∑ j  P j T S 2 j T P j T P j + P j T P j T S 2 j T P j  , (12c) . . . where P j projects on to th e sub space with j ≥ 0 double occu - pancies and S k j is defined as ( D j = j ) S k j = ∑ l 6 = j P l ( D j − D l ) k , k > 0 . (13) Calculation of th e next or ders in this way needs an increas- ing compu tational ef fort. Therefore, a comp uter pr ogram has been developed that e valuates the gen eral f ormula of Kato, ref. 18 (cf. eq. 22), up to a given or der . ˜ H ( ˜ c ) , the first terms of which are given by (1 2) is the de- sired effecti ve Hamiltonian. It is valid in any subspa ce with fixed number of j double occup ancies (of p articles co rre- sponding to ˜ c † ). For large U , the grou nd state ˜ ψ 0 ( ˜ c ) of ˜ H ( ˜ c ) must have the lowest value o f U D ( ˜ c ) , and because ˜ H ( ˜ c ) lea ves the num ber of d oubly occupied sites unchang ed, see (4), this state does not co ntain a ny d ouble oc cupancies at half b and filling, so we put j = 0 in the following. So far , the consid- erations are v alid for arbitrar y d imension. Now , we focus on the case of half-filling in infinite dimensions. Then, all glo bal singlets are ground states of ˜ H ( ˜ c ) , cf. ref. 19. Therefore, each lattice sit e is equally lik ely to be occupied by an electron with spin ↑ or ↓ , irrespective of th e spin on any other lattice site. This enables u s to p erform the g roun d-state exp ectation v al- ues in h ˜ ψ 0 ( ˜ c ) | ˜ h n | ˜ ψ 0 ( ˜ c ) i in the case of a h alf-filled band in a second computer program . Both these c omputer pro grams are of algeb raic nature (work w ith integers only ) a nd thus giv e exact results for any giv en order in 1 / U . 3 II. GRAPHS AND ALGORITHM In this section , we gi ve m ore details regard ing the evaluation of ˜ H ( ˜ c ) . In fixed order i , the operator ˜ h i containes all possible electron hops resulting f rom i + 1 applications o f T ; it con- tributes to the en ergy in the o rder 1 / U i . In the following, we will deal on ly with states with zero dou bly oc cupied sites (in the new particles), so we drop the index 0 and denote S k 0 ≡ S k and P 0 ≡ P . A sequen ce of elec tron hops d escribed throu gh an op erator chain PT S k 1 · · · S k i T P is called pr ocess 19 . Only processes tha t resto re the initial state con tribute to the ground- state energy . Con sider now a giv en process and perfo rm the sum on lattice sites in the op erators T . The individual ter ms in this sum are called “seq uences” or “diagra ms”. Becau se of the linked-cluster th eorem, we n eed to keep o nly connected diagrams. Each of these co ntains n = 2 + i − 1 2 sites, connected throug h i + 1 jumps. Diagr ams of odd o rder in t do not con - tribute at h alf fillin g f or any lattice type. Now , we specialize our considerations to the Bethe lattice. There, all clo sed paths are self-r etracing. This can be seen in fig. 1 where the possible sequences (of hops) for the lo w orders i = 1 ( n = 2) and i = 3 ( n = 3) are d isplayed o n a Bethe lattice of conne ctivity 3. In the f ollowing, we put t = t ∗ / √ Z ( Z is th e nu mber of nearest neighbo urs) and con sider the limit Z → ∞ for fixed t ∗ ≡ 1 . This limit im plies that in the en ergy , in each order in U − i , the leading order t i + 1 is taken into account. Thus, d iagrams with more than two transitions between any two given sites are suppressed at least as 1 / Z : every additional connection of already doubly c onnected sites is smaller by 1 / Z comp ared to those with only two ju mps between any two sites. Since the paths are self-retracing , they can be collapsed into ‘Butcher trees’ 20 as also indicated in the righ t of fig. 1. The position of the first h op in the sequen ce (the index j o f c j , σ at the r ightmost T ) defines th e “r oot” of the tree, cf. fig. 1. Because th ere are exactly two hops between neighb oring s ites ( Z → ∞ ), th ere is a one-to -one corresponde nce between sequences and Butcher trees. The n umber of Butcher trees built with n vertices, A [ n ] , is still moder ate for moderate n ; it is gi ven, see ref. 21, by the following recursi ve defin ition A [ n ] : = 1 n − 1 n − 1 ∑ j = 1 A [ n − j ] k max ∑ k = 1  d k ( j ) A [ d k ( j )]  , (14) with A [ 1 ] : = 1 and A [ 2 ] : = 1, and d k ( j ) is the k th element of the set of di visors of j . In table I, we g iv e the number of trees, the n umber of initial spin con figuration s, and the number of sequences for given or der of perturbation theory . Illustrating the co mplexity of the problem, figu re 2 shows the Butch er trees up to seven vertices (representing the co n- nected sites) and th us all graphs c ontributing up to ele venth order in the perturbatio n series. In ord er to illustrate our p rocedu re, we show all intermedi- ate states for all sequen ces fo r all processes in third orde r in fig. 3. I n the main part of the figu re, these states are d isplayed on a Bethe lattice with co nnectivity 3. Note th at only th ree sites a re af fected b y the hops, the two in the center and the up- per rig ht o ne o n this segment of the lattice. All other sites of the lattice are unaffected by the h ops and we can restrict o ur W ( n D 2 ) C W C ( n D 3 ) FIG. 1: Correspondence between sequenc es of hops on the l attice (left) a nd Butcher trees (right): S ho wn are the sequ ences in first ( n = 2) and third order ( n = 3) and the related Butcher trees. The first hop defines the root of the tree (encirceled). i n A [ n ] B [ n ] C [ n ] 1 2 1 2 2 3 3 2 8 20 5 4 4 32 648 7 5 9 136 45472 9 6 20 596 564488 0 11 7 48 2712 1099056000 . . . . . . . . . . . . . . . T ABL E I: Number of trees with n = ( i + 3 ) / 2 sites or order i in per- turbation theory , A [ n ] , see text. B [ n ] i s the numbe r of the i nitial spin configurations; C [ n ] is the number of all hopp ing sequences, i. e. t he number of differen t ways to realize all processes of the given order with initial spin configurations. consideratio ns to the “Butcher trees” shown in the left part o f the figure. The re, the sites of the first hop, the roots of the trees, are encircled. The processes are generated by follow- ing parts of th e hamiltonian (term s containing the sequ ence . . . P T P . . . vanish at half band filling (hf)): ˜ h hf 3 = PT S 1 T S 1 T S 1 T P − PT S 2 T PT S 1 T P . (15) This exp ression coincid es with the fourth order (h ighest avail- able) in ref. 10, eq. (6). The probab ilities for the occurr ence of the p rocesses in th e par amagnetic ph ase, mu ltiplied by the prefactors o f the pro cesses in ˜ h hf 3 (1 and − 1 here), are gi ven in the righ t colum n of figure 3. Their sum yields the co ntribution to the ground -state e nergy . The nu merical algorithm to calcu late the expectation value of the ope rators ˜ h i is based on this diagrammatic appro ach: After co nstructing all i th order Butch er trees fo r the lattice, all po ssible seque nces o n them resulting from different terms of the ˜ h i are generated thr ough a recur si ve pro cedure, within which the condition s for th e re alisation of the e lectron tra ns- fers, as the fulfillmen t of the Pauli principle and the con sider- ation of p receding hopping steps on the b ranches, are defin ed. The first electron ho p starts fr om the roo t of the gr aph (encir- cled sites in figur es 2 and 1) to a neig hbour site. In a single loop of the prog ram, possible following jumps are test ed by a subroutin e an d executed wher e ap plicable. Th ereby , th e sec- ond and la st electron h op on a bran ch has to be per formed by the same spin species as the first o ne. This r equiremen t guar- antees the restoration of th e in itial spin configuration . The 4 n = 2 n = 3 n = 4 n = 5 n = 6 n = 7 n = 5 FIG. 2: Butcher trees up t o the order i = 11. T he order expressed by the number of si tes is i = 2 n − 3. T he encirceled sites denote the root of the trees. actual number o f double o ccupancie s that enters the op erators S k , (13) , is stored and used for the computation of the factor for the given pr ocess. T he final spin configu ration determines the factor’ s sign, ( − 1 ) P , where P is the num ber of permu ted spin p airs. Su mmation yields the contribution of given order to the grou nd-state e nergy of the Mott insulator . As sh own now , th is algo rithm has b een successfully tested by com puting the gro und-state energy of the exactly solvable Falicov–Kimball mod el, a simplified Hub bard model with one immobile spin species. III. RESUL TS A. Falico v–Kimball Model For the Hamiltonian of the half-filled Falicov–Kimball mo del we refer to van Dongen’ s fundamen tal work 22 . W e ca lculated with o ur procedure the groun d-state ene rgy on the Bethe lat- tice with infinite conn ectivity (b andwidth W = 2 √ 2 t ∗ ) up to O ( t 12 / U − 11 ) . T aking t ∗ ≡ 1 as o ur energy unit in the follow- ing, we find that all con tributions in the series in 1 / U vanish, except the first: E FK 0 ( U ) = − 1 4 U + O  1 U 13  . (16) Next, w e verify 23 our resu lt (16) using the exact solution in ref. 2 2. W e start from the expression of the kin etic energy eq. 4.7; we denote th at by E FK, T 0 ( U ) . All grou nd-state en- ergies ar e gi ven as den sities (intensiv e the rmody namic vari- able). W e use their spectral rep resentation in o rder to expre ss the Green fu nction in eq . 4.7 in terms of th e density o f states z ( ε , U ) , eq. 7.5 in ref . 22. Thus E FK, T 0 ( U ) = − 2 ∞ Z 0 d ε ∞ Z 0 d ε ′ z ( ε , U ) z ( ε ′ , U ) 1 ε + ε ′ . (17) Finally , we show by nu merical integration for different choices of U between 2 and 10 that E FK, T 0 ( U ) = − 1 2 U  1 + O ( 10 − 16 )  (18) and that confirms ou r result, eq. 16. (Here , 10 − 16 is th e nu- merical accuracy .) 5 P T S 1 T S 1 T S 1 T P Probability :      !      !      !      !  1 2U 3      !      !      !      !  1 4U 3 :      !      !      !      !  1 4U 3      !      !      !      !  1 2U 3  P T S 2 T P T S 1 T P :      !      !      !      ! 1 4U 3      !      !      !      ! 1 4U 3 :      !      !      !      ! 1 4U 3      !      !      !      ! 1 4U 3 FIG. 3: The intermediate states for the two processes and all seq uences contributing to the ground-state ener gy in third order in 1 U . The arro ws correspond to the application of T . The symbol ↑ ↓  denotes a doubly occupied lattice site; the symbols  and ·  denote a hole and a singly occupied lattice site, respecti vely . In the right column, the contribution of each sequence to the ground -state energy is indicated. They sum up to − 1 2 U 3 . W e have to conclude that a ll high er order h oppin g co ntri- butions to the groun d state energy cancel. The reason may be th at on ly o ne spin species can hop in the Falicov–Kimball model. B. Hubbard Model The calculation of the groun d-state energy of the Hub bard model to the 11 th order yields E H 0 ( U ) = − 1 2 1 U − 1 2 1 U 3 − 19 8 1 U 5 − 593 32 1 U 7 − − 23877 128 1 U 9 − 44962 45 2048 1 U 11 + O  1 U 13  . (19 ) Consequently , the doub le occ upancy of the original particles 1 L D ( U ) = d E ( U ) / d U is gi ven by 1 L D H ( U ) = 1 2 1 U 2 + 3 2 1 U 4 + 95 8 1 U 6 + 4151 32 1 U 8 + + 21489 3 128 1 U 10 + 49458 695 2048 1 U 12 + O  1 U 14  . (20 ) These perturb ation-theo retical (PT) results are sho wn as solid lines in figure 4. T he comparison s of the first (second) and third ( fourth ) order PT (do tted/dashed lines) d emonstrate a fast conv ergence for th e values of U shown. The agreement with QMC (circles) and DDMRG (crosses) r esults extrap o- lated to zer o temper ature is excellent for U > 5 (smaller than the line width in fig. 4 ). As U decreases, dev ations from these (nu merical) DMFT data increase noticea bly , since re- sults of finite order PT rapidly bec ome inaccu rate as U ap - proach es U c 1 , the critical interaction. In the fo llowing, we de- scribe a meth od of how to estimate the critical coup ling U c 1 . W e assume that th e radius o f con vergence of the 1 / U ex- pansion of the energy coincide s with the critical coupling U c 1 , beyond which the insulating phase becom es stable. W e perfor m an extrapo lation o f th e co mputer-aided h igh-o rder ev aluation to infinite ord er (ePT 13 ) that exceeds former a c- curacy in ref. 13. W ith E H ( U ) = ∑ ∞ s = 1 a 2 s U 1 − 2 s , we have U c 1 = lim s → ∞ p a 2 s + 2 / a 2 s . I n figu re 5, we plot p a 2 s + 2 / a 2 s against 1 / s . As seen in the figure, the data points are fitted by a nearly straigh t line as a function of 1 / s . T aking a slight 6  0:12  0:11  0:10  0:09  0:08 E .U / 1 st order PT 3 rd order PT 11 th order PT QMC DDMRG 4 :5 5:0 5:5 6:0 0:014 0:018 0:022 0:026 0:030 0:034 U D.U / 2 nd order PT 4 th order PT 12 th order PT QMC DDMRG FIG. 4: QMC and DDMRG 9 results for the groun d-state energy E (top) and the double occupanc y D (bottom) in comparison wi th PT , see eq. (19) and (20). curvature into account in a least-squares fit, U c 1 ( s ) = r a 2 s + 2 a 2 s ≈ U c 1 + u 1 2 s + u 2 ( 2 s ) 2 , (21) one finds U c 1 = 4 . 7 6, u 1 = − 1 6 . 4712 57, and u 2 = 5 . 7 14707 2. The cr itical expo nent (fo r de tails see ref. 13) defined by E crit ( U ) ∝ ( U − U c 1 ) τ − 1 is o btained with τ = 3 . 46, th at gives support to our assumption 13 for τ = 7 / 2. The ePT estimates for the energy are strongly supported by QMC results 13 : E PT is con verged within O ( 10 − 5 ) for U ≥ 6, while the ePT p rovides an estimate for E with a precision of the same ord er above the stability edg e of the insulator . These ePT results for E have b een repro duced at U = 4 . 8, 5 . 5, an d 6 within O ( 10 − 6 ) using the self-energy f unction al approa ch/dyna mical impurity appro ach (SFT/DIA) 24 . Other methods based on DDM RG gi ve U c 1 ≃ 4 . 4 5 and τ = 2 . 5 see ref. 7,9. As seen fro m figu re 4, a h igh accuracy of d ata is indispensable for a correct analysis of the transition. CONCLUSIONS Using K ohn’ s un itary transformation, an 1 / U -expansion fo r the Hubb ard model was derived up to order ( 1 / U ) 11 at zero temperatur e. The expansion has been fo rmulated in terms of diagramm atic rules fo r th e calculation of the groun d-state energy and the resulting doub le o ccupancy . These rules re- duce the calculation of finite-order contributions to an alge- bra wh ich becom es increasin gly complex for h igher or ders. 0 0:05 0:10 0:15 0:20 0:25 0:30 0 1 2 3 4 5 1=s U c .s / FIG. 5: Construction of ePT : PT v alues for U c ( s ) = p a 2 s + 2 / a 2 s (squares) are extrap olated to 1 / s → 0 using a quadratic least squares fit (solid line). E v aluations at smaller 1 / s (circles) define ePT coef fi- cients to all orders, cf. ref. 13. Any step of th e ru les is car ried out exactly b y our computer progr am. Explicit results were ob tained for the grou nd-state energy up to 11 th order in 1 / U . An inspection of th e contributions of the diagrams sho ws that th ere are group s of dominan t ones, namely the widespread diagrams (first ones in each order in fig. 2), an d they a re sig- nificant in v iew of the metal-insu lator transition . Th is sho uld be analysed quantitatively in the large order limit. Then, ev en an exact determination of, e. g ., U c 1 might be possible. Acknowledgments W e thank the r eferee for m any helpful sug gestions that greatly improved the p resentation of this work and E. Jeckelmann and—in particular—W . and V . Apel for m any discussions and help in revising this manuscript. Appendix: FLO WCHAR T OF THE ALGORITHM The kerne l of the program is the recursi ve pr ocedur e ho pping which is calling the proced ure hopp ing_f wd and vice versa; their flowcharts ar e shown in figure 6. ☛ ✡ ✟ ✠ BEGIN ✲ hopping ✲ ☛ ✡ ✟ ✠ END ♠ ❄ ✻ hopping_fwd The main task o f these “ (non- linear) indirect- recursive” proced ures is to find all electron hopping processes generated by the Hamiltonian, which occur on a set of graphs, as sho wn in the pr eceding sectio n. The grap hs are defined in su ch a way that the set of all p ossible pro cesses of the g iv en order is identical to th e set of processes starting fro m the root of the graphs. All possible spin config urations o n the g raphs are gen- erated and stor ed in arrays. Due to sym metry , it is suf ficient 7 hopping(hop,bran ch,doubleocc,pf, pt,ps, spinconfig,jumph ist,direction,do ublehist) ☛ ✡ ✟ ✠ BEGIN ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ branch>0 No Yes ✲ ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ jumphist[branch, not direction] in [ps,hole] No Yes ❄ ☛ ✡ ✟ ✠ END ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ direction=true No Yes ❄ Refresh of the ’spinconfig’ data structure ❄ Refresh of the ’spinconfig’ data structure ❄ jumphist[branch,direct ion]:=ps doublehist:=doublehist + char(48+doubleocc) ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ nr=2*branch_count No Yes ✛ ✲ Computation of the number of permutations ✛ Return of result to the evaluation procedure ❄ ❄ i:=1 ❄ ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ i<=branch_count No Yes ✛ ✲ ✲    ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ nr>0 No Yes ❄ ❄ hopping_fwd(i,nr+1,fal se) ✻ ✲    ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ nr>0 or ((nr=0) and (tree[branch _nr].p=0)) No Yes ❄ ❄ hopping_fwd(i,nr+1,tru e) ✛ ✛ i:=i+1 ✛ ✻ hopping_fwd(i,nr +1,direction) ☛ ✡ ✟ ✠ BEGIN ❄     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ jumphist[branch, direction]= hole Yes No ❄ ☛ ✡ ✟ ✠ END ✛     ❅ ❅ ❅ ❅     ❅ ❅ ❅ ❅ direction=true Yes No ❄ Initialisation of ’_from’ and ’_to’ ✻ Initialisation of ’_from’ and ’_to’ ❄    ❅ ❅ ❅    ❅ ❅ ❅ _from=up Yes No ❄    ❅ ❅ ❅    ❅ ❅ ❅ _to=hole Yes No ✲ hopping(hop,branch,dou bleocc,hole,up,up, spinconfig,jumphist,tr ue,doublehist) ✲ ✲   ❅ ❅ ❅    ❅ ❅ ❅ _to=down Yes No ✲ hopping(hop,branch,dou bleocc+1,hole,pair,up, spinconfig,jumphist,tr ue,doublehist) ✲ ✲ ❄    ❅ ❅ ❅    ❅ ❅ ❅ _from=down Yes No ❄    ❅ ❅ ❅    ❅ ❅ ❅ _to=hole Yes No ✲ hopping(hop,branch,dou bleocc,hole,down,down, spinconfig,jumphist,tr ue,doublehist) ✲ ✲   ❅ ❅ ❅    ❅ ❅ ❅ _to=up Yes No ✲ hopping(hop,branch,dou bleocc+1,hole,pair,down, spinconfig,jumphist,tr ue,doublehist) ✲ ✲ ❄    ❅ ❅ ❅    ❅ ❅ ❅ _from=pair Yes No ❄    ❅ ❅ ❅    ❅ ❅ ❅ _to=up Yes No ✲ hopping(hop,branch,dou bleocc,up,pair,down, spinconfig,jumphist,tr ue,doublehist) ✲ ✲   ❅ ❅ ❅    ❅ ❅ ❅ _to=down Yes No ✲ hopping(hop,branch,dou bleocc,down,pair,up, spinconfig,jumphist,tr ue,doublehist) ✲ ✲   ❅ ❅ ❅    ❅ ❅ ❅ _to=hole Yes No ✲ hopping(hop,branch,dou bleocc-1,down,up,up, spinconfig,jumphist,tr ue,doublehist) ❄ ❄ hopping(hop,branch,dou bleocc-1,up,down,down, spinconfig,jumphist,tr ue,doublehist) ✲ ✲ ✲ FIG. 6: Fl o wcharts of the proced ures hopping and hopping_fwd . 8 to occupy the root always by an u p-spin lead ing to 2 n − 1 con- figuration s, wher e n is the numb er of vertices of the graphs. For the descriptio n o f th e procedu res the following variables are used, cf. fig. 6: hop : nu mber of the current hop branc h : nu mber of the branch, on which the hop occurs doubl eocc : n umber of double occupanc ies pf , pt , ps : spin states ( hole , up , d own , pai r ): p f spin state of the site o ut of which the jump occurs, pt spin state of the site to which the jump occurs, ps jump ing spin spinc onfig : structure describ ing the current spin configura- tion on the graph jumph ist : table con taining the jump history on the branches direc tion : d irection of th e jump, defined by d escription o f the graph: true jump forward, false jump back ward doubl ehist : sequence o f digits repr esenting the h istory of the number of double occupan cies branc h_cou nt : total number of branches in the graph The initial call of the proc edure h oppin g is done with the following par ameters: hop= 0 , bra nch=0 , double occ=0 , pf=up , pt=down , ps=up , s pinco nfig , jump hist , t rue , ’0’ . The pr ocedur e executes jump s on all branche s in both direction s; the possibility of a jump is tested throug h the pro cedure hopping_ fwd . Its ca ll is done with hopping _fwd( branch_nr,hop+1,true) (fo rwards) or hopp ing_f wd(br anch_nr,hop+1,false) (back wards) and it checks if on the curr ent bran ch a hop has already oc- cured in the given direction . I f not, d epend ing on the spin state of the inv olved s ites, the proced ure hopping is called, th e ar- ray of spin states is refresh ed, an d a n ext br anch is tested. The second jump on a branch h as always to be perfo rmed b y the same spin species as the first one. The recu rsion has the prope rty that, in case of exiting the proced ure when the jump was not p ossible, the values of variables resulting from prec eding steps a re autom atically re- stored. This construction of the algo rithm gu arantees th at all po s- sible variants o f elec tron jum ps are tested in accor dance with the assum ptions, and that the final spin configuration equals the initial spin configura tion; therefo re this condition d oes n ot need not to be tested. Af ter the last step (carr ying the number 2*bra nch_c ount ), the characteristic factor f proc. = ( − 1 ) P N spin [ proc. ] 2 n − 1 ∑ m f m i ∏ j = 1  1 d j  k j (A.1) for a process appearing in ord er 1 / U i is computed. I n (A.1) P denotes the n umber of permuted spin pairs, N spin is the n um- ber of spin con figurations n ot chang ed b y the process. Th e sum runs over all terms of the h amiltonian which ge nerate the p rocess an d f m is th e related factor o btained fro m eq ua- tion (7) and calculated by a sepa rate algorithm. d j are the number s o f double occupancies and k j are also obtained from equation (7). The sum of the factors yields the final result. 1 A. Montorsi, ed., T he Hubbar d Model: A R eprint V olume (W orld Scientific, Singapore, 1992). 2 F . Gebhard, The Mott Metal-Insulator T r ansition (Springer , Berlin, 1997). 3 J. Schlipf, M. Jarrell, P . G. J. v an Dongen, N. Blümer , S. K ehrein, T . Pruschke, and D. V ollhardt, Phys. Rev . Lett. 82 , 4890 (1999). 4 M. J. Rozenber g, R. Chitra, and G. K otliar , Phys. Rev . Lett. 83 , 3498 (1999). 5 W . Krauth, P hys. Re v . B 62 , 6860 (2000). 6 G. Kotliar , E. Lange, and M. J. Rozenberg, Phys. Rev . Lett. 84 , 5180 (2000). 7 M. P . Eastwood, F . Gebhard, E. Kalino wski, S. Nishimoto, and R. M. Noack, Eur . Phys. J. B 3 5 , 155 (2003). 8 M. Feldbacher , K. Held, and F . F . Assaad, P hys. Rev . Lett. 93 , 136405 (2004). 9 S. Nishimoto, F . Gebhardt, and E. Jeckelmann, J. Phys.: Con- dens. Matter 16 , 7063 (2004). 10 A. L. Cherny she v , D. Galanakis, P . Phill ips, A. V . Rozhko v , and A.-M. S. Tremb lay , Phys. Rev . B 70 , 235111 (2004). 11 J.-Y . P . Delanno y , M. J. P . Gingras, P . C. W . Holdsworth, and A. - M. S. T remblay , P hys. Re v . B 72 , 115114 (2005). 12 P . Phillips, D. Galanakis, and T . D. St anescu, Phys. Re v . Lett. 93 , 267004 (2004). 13 N. Blümer and E. Kalino wski, Phys. Rev . B 71 , 195102 (2005). 14 N. T eichmann, D. Hinrichs, M. Holthaus, and A. E ckardt, Phys. Rev . B 79 , 100503 (2009). 15 A. Eckardt, Phys. Rev . B 79 , 195131 (2009). 16 W . Koh n, P hys. Re v . 133 , A171 (1964). 17 M. T akahashi, J. Phys. C: Solid State Phys. 10 , 1289 (1977). 18 T . Kato, Prog. Theor . Phys. 4 , 514 (19 49). 19 E. Kalinowsk i, Ph.D. t hesis, Fachbe reich P hysik, Philipps- Univ ersität Marbu rg (2002 ). 20 J. C. Butcher , The Numerical Analysis of Ordin ary Differ ential Equations (W iley , Chichester , 1976). 21 Number of tr ees w ith n unlabeled nodes , The On-Line Encyclop edia of Integer Seque nces, research.att.com/ njas/sequence s/A000081 . 22 P . G. J. van Don gen, Phys. Rev . B 45 , 2267 (1992). 23 W e thank W . Apel for a discussion of t his point. 24 K. Požgaj ˇ ci ´ c, arXi v:cond-mat/0407 172v1 (2004).