First-order metal-insulator transitions in the extended Hubbard model due to self-consistent screening of the effective interaction
M. Schüler, E. G. C. P. van Loon, M. I. Katsnelson, T. O. Wehling
aa r X i v : . [ c ond - m a t . s t r- e l ] A p r First-order metal-insulator transitions in the extended Hubbard model due toself-consistent screening of the effective interaction
M. Sch¨uler,
1, 2, ∗ E. G. C. P. van Loon, M. I. Katsnelson, and T. O. Wehling
1, 2 Institut f¨ur Theoretische Physik, Universit¨at Bremen, Otto-Hahn-Allee 1, 28359 Bremen, Germany Bremen Center for Computational Materials Science,Universit¨at Bremen, Am Fallturm 1a, 28359 Bremen, Germany Radboud University, Institute for Molecules and Materials,Heyendaalseweg 135, NL-6525 AJ Nijmegen, The Netherlands (Dated: October 2, 2018)While the Hubbard model is the standard model to study Mott metal-insulator transitions, it isstill unclear to what extent it can describe metal-insulator transitions in real solids, where nonlocalCoulomb interactions are always present. By using a variational principle, we clarify this issue forshort- and long-range nonlocal Coulomb interactions for half-filled systems on bipartite lattices.We find that repulsive nonlocal interactions generally stabilize the Fermi-liquid regime. The metal-insulator phase boundary is shifted to larger interaction strengths to leading order linearly withnonlocal interactions. Importantly, nonlocal interactions can raise the order of the metal-insulatortransition. We present a detailed analysis of how the dimension and geometry of the lattice as wellas the temperature determine the critical nonlocal interaction leading to a first-order transition: forsystems in more than two dimensions with non-zero density of states at the Fermi energy the criticalnonlocal interaction is arbitrarily small; otherwise, it is finite.
PACS numbers: 72.80.Rj; 73.20.Hb; 73.61.Wp
I. INTRODUCTION
The Hubbard model is a central model for under-standing various aspects of strongly correlated electrons.It incorporates the competition between kinetic and in-teraction energies in the most basic way and exhibitsphenomena such as magnetism and metal-insulator tran-sitions with and without magnetic transitions . How-ever, particularly due to neglecting nonlocal interaction,the Hubbard model can be quite far from describingreal materials whenever nonlocal interactions are not effi-ciently screened, e.g., in two-dimensional materials. Oneexample is the plasmon dispersion in metals which dif-fers qualitatively in models with and without nonlocalinteractions . Most obviously, in insulating systems,where screening is by definition incomplete, prominentnonlocal interaction effects should be expected. It isthus unclear whether the Hubbard model can describethe Mott metal-insulator transition (MIT) even qualita-tively correctly.Indeed, the question about the order of the MIT hasbeen controversial for about five decades . In theHubbard model with strictly local interaction, the orderof the MIT depends on the degree of magnetic frustra-tion in the system. If magnetic order is fully suppressed,the transition is of first order below a critical tempera-ture T c as, e.g., dynamical mean-field theory (DMFT) and related quantum cluster theories have demon-strated. Otherwise, the MIT is accompanied by mag-netic (quasi)order and is continuous . Thus, Hubbardmodels on bipartite lattices, like the honeycomb, square,diamond, and cubic lattice, as well as higher dimensionalgeneralizations thereof, feature continuous MITs. Dueto the various simplifications implied by the Hubbard model, it is, however, unclear how well this picture of theMIT relates to the experimentally realized one. Alreadyin his original work, Mott, for instance, argues that thephysically realized MIT should, due to the long-range na-ture of the Coulomb interaction, be of first order , whichis indeed found experimentally in many transition-metaloxides . In this context, the question of how the Hub-bard model’s MIT is connected to that of the extendedHubbard model, which includes nonlocal interactions, ishighly relevant.Here, we show that the MIT in the half-filled extendedHubbard model on bipartite lattices is of first order fornonlocal interactions larger than a critical V c , as depictedschematically in Fig. 1(a). V c depends on the dimensionand lattice topology and can be even arbitrarily small incubic systems in d >
2. The first-order transition canbe masked by a charge density wave [CDW; Fig. 1(b)], asituation which we find, e.g., in the honeycomb lattice.We propose (and later substantiate) the followingmechanism for how nonlocal interactions induce a dis- U V (b) metal insulator CDWU V (a) metal insulator V c (U)CDW Figure 1. (Color online) Schematic phase diagram of the ex-tended Hubbard model: (a) The continuous (solid) and first-order (dotted) metal-insulator transition lines touch at V c ( U ).A coexistence region surrounds the first-order transition. Forlarge V /U a CDW phase occurs. (b) The CDW phase canmask the first-order MIT. continuous MIT: Nonlocal interactions generally decreasecorrelations in half-filled extended Hubbard models .The amount of decrease, due to different screening ,is larger in the metallic regime than in the insulatingregime. Now consider two systems close to the MIT withinitially no nonlocal interactions: one metallic and oneinsulating. Nonlocal interactions V will push the MITof the formerly metallic system to larger local interac-tion U met. c > U ins. c than the formerly insulating system,resulting in a discontinuous (i.e., first-order) MIT at suf-ficiently large V . II. MODELS AND METHODSA. Hubbard model
To begin, we briefly review the MIT on bipartite lat-tices in the Hubbard model [Eq. (2)], i.e., without non-local interactions. For the systems here, in d > U induces a MIT from a paramagneticmetal to an antiferromagnetic insulator . For latticeswith perfect nesting the critical interaction U c vanishesfor zero temperatures . For lattices with a vanishingdensity of states (DOS) at the Fermi energy E F , U c is fi-nite for T = 0 . In two dimensions the Mermin-Wagnertheorem prevents long-range order and therefore a mag-netic transition . However, quasi-long-range antifer-romagnetic fluctuations lead to a similar phase diagramfor the transition from a paramagnetic metal to a quasi-ordered insulator. Again, perfect nesting in the squarelattice leads to U c = 0 for T = 0 , while the vanishingDOS at E F in the honeycomb lattice leads to U c > T = 0 . In all cases the gap appears continuously;that is, the MIT is not of first order. B. Extended Hubbard model
We now turn to the extended Hubbard model, H = − t X h i,j i ,σ c † iσ c jσ + U X i n i ↑ n i ↓ + 12 X i = j V ij n i n j , (1)where c ( † ) iσ is the annihilation (creation) operator for elec-trons on site i and spin σ , t is the nearest-neighbor hop-ping amplitude, U is the local interaction, and V ij is thenonlocal interaction between electrons at sites i and j . n iσ and n i are the spin-resolved and total occupationoperators, respectively. We focus on nearest-neighbor(NN), V j = V δ , and long-range (L), V j = V /r j ,interactions .In this model repulsive nonlocal interactions decreasecorrelations and can lead into a CDW phase .Our focus is solely on how nonlocal interactions influencethe order of the MIT. C. Variational principle
We investigate the U - V - T phase diagram of theextended Hubbard model by approximating its ther-modynamic ground state using the Peierls-Feynman-Bogoliubov variational principle with a Hubbardmodel as the effective system . The effective Hubbardmodel, reading˜ H = − t X h i,j i ,σ c † iσ c jσ + ˜ U X i n i ↑ n i ↓ , (2)is varied via the effective local interaction ˜ U in order tominimize a free energy functional. Therefore ˜ U is˜ U = U + X j =0 V j ∂ ˜ U h n n j i ˜ H ∂ ˜ U h n n i ˜ H . (3)Although the variational principle provides only an upperbound of the free energy, it has been found to give an ac-curate description of the physics present in the effectivereference Hubbard model and even gives exact dou-ble occupancies for infinitesimal nonlocal interactions .This makes the approach appropriate for capturing theMIT, a hallmark of Hubbard model physics, in theextended Hubbard model. We introduce the effectivescreening factors α NN ( ˜ U ) = − ∂ ˜ U h n n i ˜ H ∂ ˜ U h n n i ˜ H and α L ( ˜ U ) = − P j =0 1 r j ∂ ˜ U h n n j i ˜ H ∂ ˜ U h n n i ˜ H with which Eq. (3) simplifies to˜ U = U − V α
NN/L ( ˜ U ) , (4)for the case of NN and L interactions. Here, α ( ˜ U ) is aproperty of the effective Hubbard model and quantifiesthe above mentioned decrease of correlations by V : Non-local interactions shift the transition at ˜ U MIT to leadingorder linearly with a slope of α − ( ˜ U = ˜ U MIT ); that is,a positive α leads to nonlocal interactions stabilizing themetallic phase. Indeed, we find strictly positive α in ournumerical calculations (see Fig. 2).From Eq. (4) we calculate the change in the effectiveinteraction ˜ U with V and U ,( ∂/∂U, ∂/∂V ) ˜ U = (cid:18) V ∂α∂ ˜ U (cid:19) − (1 , − α ) . (5)The derivatives diverge for − V c ∂ ˜ U α . At this point wefind a bifurcation of the solution of Eq. (4); that is, a sin-gle point ( U, V ) is mapped to multiple ˜ U and thus a jumpof the observables of the effective system. Particularly,discontinuities in α ( ˜ U ) lead to arbitrarily small nonlocalCoulomb interactions inducing a first-order phase tran-sition ( V c = 0).In the following we numerically determine the critical V c required to induce such a first-order phase transition.To this end we calculate α for different lattice topolo-gies and dimensions: We choose to investigate the squareand honeycomb lattices and their three-dimensional gen-eralizations the cubic and diamond lattices. We rely on (c) square ̃Ũt (d) square α NN (a) honeycomb ̃Ũt α L (b) honeycombβt Figure 2. (Color online) Effective screening factor α as definedabove Eq. (4) for NN (top panels) and L (bottom panels)interaction. Left (right) panels show the case of the honey-comb (square) lattice. We show results from low (red) to high(cyan) inverse temperatures βt . quantum Monte Carlo simulations of all effective Hub-bard models.For the case of two dimensions we use the determinantquantum Monte Carlo method (DQMC) implementedin the quest code . We alleviate finite-size and Trot-ter errors by extrapolating from finite Trotter discretiza-tions of ∆ τ = 0 .
2, 0 .
1, and 0 .
05 and linear lattice sizes of L = 8, 10, and 12 for the square lattice and L = 6, 9, and12 for the honeycomb lattice . We evaluate derivativesin the definition of α numerically by solving Hubbardmodels in steps of ∆ ˜ U /t = 0 . III. RESULTSA. Honeycomb lattice
First, we discuss the case of the honeycomb lattice forwhich the α ( ˜ U ) are plotted for βt = 2 to βt = 10 inFigs. 2(a) and (b). For βt & E F . Thus, we can draw conclusions for finite tem-peratures and T →
0. For low temperatures α showsa minimum at ˜ U min /t ∼ . U kink /t ∼ .
2. The MIT is in between at˜ U MIT /t ≈ . .Notably, the dependence of the effective screening fac-tor α on ˜ U is rather weak (no steplike features). Hence,there is no ˜ U where the slope ∂α/∂ ˜ U is particularly large,and from Eq. (5) we expect that rather large V c wouldbe needed to push the MIT to first order here. Clarifyingthis, we solve Eq. (4) and calculate the U -dependent dou-ble occupancy of an extended Hubbard model with dif-ferent nearest-neighbor interaction at βt = 10 [Fig. 3(a)]. For increasing V the V = 0 line is shifted towards larger U (i.e., V weakens correlations).Concerning the influence of nonlocal interactions onthe order of the transition, we calculate V c from the slopeof α . We find V c /t ≈ . V c /t ≈ . V /t ∼ . since larger V stabilize aCDW phase which we estimate in strong coupling, aspresented in Appendix D. For the honeycomb lattice V c is always larger than V CDW such that no first-order MITwill be observable.We infer a U - V phase-diagram schematically shown inFig. 1(b). The slope of the transition line at V = 0 isgiven by 1 /α ( U MIT ), with α NN(L) ( U MIT ) = 0 . . reveal a slope compatible to α L . .
55; For nearest-neighbor interactions dynamical clus-ter approximation (DCA) calculations indicate α NN ∼ .
23, whereas QMC calculations lead to α NN . . B. Square lattice
The case of the square lattice turns out to be differ-ent. We show α NN/L in Figs. 2(c) and (d) for βt = 2to βt = 20. Here, α is strongly temperature dependent.Prominently, ˜ U min and ˜ U kink are shifted to smaller ˜ U ,and the slope at ˜ U kink gets steeper with increasing β .The increase in the slope of α ( ˜ U ) traces back to the de-velopment of a soft kink in h n n i ( ˜ U ). Comparing U min and U kink to the temperature dependent critical interac-tion ˜ U MIT from Ref. 21, we find that ˜ U min ≈ ˜ U MIT andthat ˜ U kink approaches ˜ U MIT with lower temperatures.Concerning the resulting phase diagram [Fig. 1(a)],the slope of the U - V phase-transition line at V = 0 is1 /α ( ˜ U MIT ), with α NN(L) ( ˜ U MIT ) ∼ . at βt = 6 and combined GW plus extendedDMFT at βt = 25 are compatible with α ∼ . α ∼ . U/t(b) square, βt = 20 ⟨ n ↑ n ↓ ⟩ U/t(a) ⟩oney⟨omb, βt = 10V/t = 0V/t = 0.6V/t = 1.2V/t = 1.8
Figure 3. (Color online) Double occupancy for (a) the honey-comb ( βt = 10) and (b) square lattice ( βt = 20) for V /t = 0(black), 0 . . . tion V c leading to a first-order transition we calculatemax | ∂ ˜ U α ( ˜ U ) | , which turns out to exhibit a linear β de-pendence, as can be seen in Fig. 4(a). Thus, by ex-trapolating to decreasing temperatures [Fig. 4(b)] we ex-pect a first-order phase transition at decreasing nonlocalinteraction strengths: V c → T →
0. Here, wefind V c /t ≈ .
15 for βt = 20. For NN interaction thefirst-order phase transition is probably not observable for βt = 20 since it is deeply buried in the CDW phase (seeAppendix D). However, since V c scales linearly with T but the nonlocal interaction for the CDW phase scalesas V CDW /t = 4 π [ln(8 t/T )] − , [black line in Fig. 4(b)],low enough temperatures always lead to a favoring of thefirst-order MIT over the CDW. Moreover, long-range in-teractions partially suppress the CDW phase such thatin this case the first-order MIT will be observable at onlyslightly lower temperatures (see Appendix D).
10 15 20 βt m a x | ∂ α / ∂ ̃ U | (a) 0.04βt + 0.04 T/t V c / t (b) (0.04t/T + 0.04) −1 CDW̃weak
Figure 4. (Color online) (a) Dependence of max[ | ∂α/∂ ˜ U | ]on βt for the square lattice with linear fit (red line). (b)Corresponding V c ( T ) = − [ ∂α/∂ ˜ U ] − . Results for NN and Lcoincide within error bars. The fading black curve is a weak-coupling estimate of the critical V for the CDW valid for small V . In Fig. 3(b) we present the double occupancy depen-dent on U at βt = 20 for different nonlocal interactions.The curves are shifted to larger U with increasing V ,where the different amounts of shifting are apparent forthe metallic and insulating regimes. This different ef-fective screening on the metallic and insulating sides ofthe transition eventually lifts the MIT to first order here.The soft kink visible for V = 0 (black solid line) gets asteplike shape for V >
0, which for
V > V c ≈ . t even-tually leads to unphysical (see below) loops, as shown indetail in Fig. 3(d) for V = 1 . t . The real double occu-pancy in the coexistence region is obtained by Maxwellconstruction and shown as a dashed line. This coexis-tence region is shown schematically in Fig. 1(a).The double occupancy D = h n ↑ n ↓ i and the localHubbard interaction U are conjugate variables in thethermodynamic sense, i.e., D = N ∂F/∂U , where F isthe free energy and N is the number of lattice sites inthe system. Since F is not only extremal but actuallyminimal in a stable thermodynamic state, a small devia-tion from the thermodynamic ground state must increasethe free energy. The resulting thermodynamic inequal-ity ∂D/∂U < D ( U ) curve inside the hysteresis re-gion [Fig. 3(d)], which thus corresponds to thermody-namically unstable states. This behavior is characteris-tic of a first-order transition and signals that the metallicand the insulating sides of the transition are not linkedcontinuously through a series of thermodynamically sta-ble states. C. Cubic and diamond lattices
We now turn to higher-dimensional systems withcubic and diamond lattices, which generalize thesquare and honeycomb lattice to three dimensions. Whilethe diamond lattice preserves the linearly vanishing DOSat E F , the cubic lattice loses the van Hove singularity at E F but keeps a nonzero DOS at E F . The main differencewith d = 2 is the absence of the Mermin-Wagner theoremand the presence of finite-temperature antiferromagneticlong-range order. We solve the Hubbard models in d = 3in DMFT using triqs . We allow for antiferromag-netic long-range order to study the thermodynamicallyrelevant transition from a (semi)metal to an antiferro-magnetic insulator . We provide raw data, details onthe simulations, and results for the case of infinite dimen-sions in Appendix B. From DQMC simulations in d = 3we find that the discontinuous behavior of α is essentiallydetermined by ∂ ˜ U h n ↑ n ↓ i (see Appendix C for details)and thus search for discontinuities in ∂ ˜ U h n n i directly,circumventing the calculation of ∂ ˜ U h n n j i for j > m cubic diamond ̃Ũt −0.06−0.04−0.02 ∂ ̃ ∂ ̃ U ⟨ n n ⟩ βt = ⟩βt = 5 βt = ⟨0βt = 20 ⟩ 4 5 ̃Ũt βt = 5βt = ⟨0 βt = 20βt = 40 Figure 5. (Color online) DMFT results for the cubic (left pan-els) and diamond (right panels) lattices. From top to bottomthe panels show staggered magnetization m and ∂ ˜ U h n n i . The results for different temperatures are presentedin Fig. 5. The onset of a finite staggered magnetiza-tion m determines the critical effective interaction ˜ U c atthe MIT. For the diamond lattice, ˜ U c becomes temper-ature independent for small enough temperatures likein d = 2. For the cubic lattice U c approaches zero for T → . We find clear discontinuities in the doubleoccupancies’ derivative and thus find discontinuities in α for all temperatures only for the cubic lattice. Thesize of the discontinuity grows as the system approacheslarge N´eel temperatures . From this discontinuity weconclude that infinitesimal positive nonlocal interactionsinduce a first-order phase transition in three or more di-mensions; that is, we expect V c = 0 for d > E F in the caseof the diamond lattice leads to no discontinuities at lowtemperatures. For large enough temperatures ( βt . ∂ ˜ U h n n i appears. We conclude that forlow temperatures only finite nonlocal interactions inducefirst-order phase transitions in diamondlike lattices in ar-bitrary dimensions.This dimensional dependence of the MIT in the cubicsystems can be understood from the nature of the anti-ferromagnetism. An antiferromagnetic phase transitiontranslates to a kink in the double occupancy since the lat-ter is related to the magnetic moment as m = n − n ↑ n ↓ .In d = 2, the Mermin-Wagner theorem forbids conven-tional (i.e., second order in the Ehrenfest sense) antiferro-magnetic phase transitions at finite temperature, whichleads to a smooth double occupancy and a finite V c . For d >
2, on the other hand, an antiferromagnetic phasetransition leads to a kink in the double occupancy anda first-order phase transition at infinitesimal V . The van-ishing DOS for the diamond lattice, on the other hand,leads to an unusual critical behavior , and thus no kinkin the double occupancy and a finite V c . IV. CONCLUSION
We showed that nonlocal interactions in bipartite ex-tended Hubbard models can lead to a first-order MIT.This result is highly relevant in the context of the ques-tion of whether Hubbard models describe discontinuousMITs occurring in realistic materials. The underlyingmechanism is governed by nonlocal interactions screen-ing correlations differently in the insulating and metallicphases, with the metallic phase being generally stabilizedby the nonlocal interactions. Interestingly, this is in con-trast to the mechanism envisioned by Mott , which isbased on nonlocal interactions stabilizing the insulating phase. We found first-order transitions for nonlocal in-teractions larger than a critical V c ( U ). Our calculationsindicate V c = 0 for cubic systems in d >
2, whereas sys-tems with vanishing DOS at E F (e.g., diamond) and two-dimensional systems in general show V c > , explaining how the continuous MITin Hubbard models is reconvened with the discontinuousMIT in real materials. Acknowledgments.
We acknowledge the kind help ofR. Staudt. We thank S. Haas, S. Wessel, A. Rosch,and F. Gebhard for comments. E.G.C.P.v.L. and M.I.K.acknowledge support from ERC Advanced Grant No.338957 FEMTO/NANO. Computer time at the HLRNis acknowledged.
Appendix A: Details of DQMC simulations
We perform simulations at a fixed value of ˜ U and athalf filling, i.e., by setting µ = ˜ U /
2. As discussed inRef. 57, simulations of the Hubbard model for
U/t & U/t .
6, we include global moves as a precau-tionary measure and indeed find no “sticking” behavior ofthe occupancies described in Ref. 57. With 500 warm-upsweeps we perform between 10000 and 1 million measure-ment sweeps depending on the temperature, lattice size,and Trotter discretization. We provide our raw data to-gether with more equal time measurements provided bythe quest software for all lattice sizes, Trotter discretiza-tions, and temperatures for both the square and honey-comb lattices together with error estimates and a com-plete set of input parameters on the Zenodo platform .We deal with finite-size and finite Trotter errors byextrapolating schemes: We simulate effective Hubbardmodels in d = 2 with imaginary time discretizations∆ τ = 0 .
05, 0 .
1, and 0 . τ . We present an example of this pro-cedure in Fig. 6. The errors for finite ∆ τ do not existfor ˜ U = 0 and increase for larger ˜ U . From linear systemsizes L = 8, 10, and 12 for the square lattice and L = 6,9, and 12 for the honeycomb lattice we extrapolate a lin-ear dependence on L − . We show an example in Fig. 7.The finite-size errors in the charge correlation functionget smaller for larger ˜ U since the electrons localize. Δτ ⟨ n n ⟩ ̃U⟨t = 1.0 rawlinear fitextrapol. ̃Ũt ⟨ n n ⟩ Δτ = 0.2Δτ = 0.1Δτ = 0.05Δτ → 0
Figure 6. (Color online) Finite-∆ τ extrapolation for the near-est neighbor charge correlator on the square lattice for βt = 10and L = 10. Left: linear fit (blue line) to raw data on ∆ τ − (black dots) resulting in extrapolated value at ∆ τ → U/t = 1 .
0. Right: ˜ U -dependent results for the ex-trapolation. In order to reduce the inherent noise in the MonteCarlo data, which poses a serious problem when calcu-lating derivatives with respect to ˜ U , we use a Savitzky-Golay approach ; that is, we analytically take deriva-tives of polynomials which are locally fitted in a win-dow of width w to the numerical values of h n n j i ( ˜ U ).We show an example of this procedure for two differ-ent cases: smooth dependence on ˜ U with little noise forhigh temperatures ( βt = 4 .
0; Fig. 8) and rather large L −2 ⟨ n n ⟩ ̃U⟩t = 1.0 rawlinear fitextrapol. ̃Ũt ⟨ n n ⟨ L = 8L = ⟨0L = ⟨2L → ∞
Figure 7. (Color online) Finite- L extrapolation for the nearestneighbor charge correlator on the square lattice for βt = 10and ∆ τ = 0 .
2. Left: linear fit (blue line) to raw data on L − (black dots) resulting in the extrapolated value at ∆ τ → U/t = 1 .
0. Right panel: ˜ U -dependent resultsfor the extrapolation. dependence on ˜ U with large noise for low temperatures( βt = 20 .
0; Fig. 9). We show results for different fit win-dows ( w = 0 . w = 1 .
0) for cubic polynomials. Inall cases, the raw data and the smoothed data are hardto distinguish. However, the derivative with respect to˜ U vastly increases noise for the raw data. Taking thederivative analytically for the smoothed data avoids this.The high-temperature data show that the larger windowleads to smoother results. The case of low temperature,however, exemplifies the drawback of too large windows:The steep feature in ∂ ˜ U h n n i at ˜ U /t ∼ .
3, which canbe clearly seen in the raw data, is washed out for w = 1 . w = 0 . U /t but alsoshows larger noise for larger ˜
U /t . Since the derivative of ∂ ˜ U h n n i (via that of α ) determines the critical V c , wehave extracted them with w = 0 . U w /t = 1 . U /t does not obstruct the trends visible in Fig. 2.A word on the error bars shown in Figs. 8 and 9: Thequantities (with error bars) actually measured in theDQMC algorithm are h n ↑ n i ↓ i and h n ↑ n i ↑ i , such thatwe obtain h n n i i by summing over the two observables.Since the two constituent observables are correlated (for i = 0), the error bar on their sum cannot simply beobtained by Pythagorean addition. Visual inspection ofthe raw data at neighboring values of ˜ U shows that er-ror cancellation happens in the determination of h n n i i (equivalently, the statistical errors in h S z S zi i are largerthan the Pythagorean sum of the error bars of h n ↑ n i ↓ i and h n ↑ n i ↑ i ). In the figures, we have done Pythagoreanaddition and scaled the errors, such that the error barsenclose 70% of the smoothed data, which leads to visu-ally reasonable results but overall too large error bars for˜ U /t &
3. Since the error estimates do not influence thecalculations, this is, however, not a crucial point. ⟨ n n ⟩ ⟨ raww = 1.0w = 0.4 ⟨ n n ⟩ ⟨ raww = 1.0w = 0.4 ̃U⟩t −0.0̃−0.02−0.01 ∂ ̃ U ⟨ n n ⟩ ⟨ central⟨diff.w = 1.0w = 0.4 ̃U⟩t ∂ ̃ U ⟨ n n ⟩ ⟨ central⟨diff.w = 1.0w = 0.4 Figure 8. (Color online) Raw data and smoothed data for on-site and nearest-neighbor charge correlators ( h n n i , h n n i ;top panels) and their derivatives with respect to ˜ U (bottompanels) for the square lattice with L = 12, ∆ τ = 0 .
2, and βt = 4 .
0. We show smoothing results of different fit windows, w = 1 . w = 0 .
4, with cubic polynomials. Where no errorbars are visible, they are hidden behind the markers. ⟨ n n ⟩ ⟨ raww = 1.0w = 0.4 ⟨ n n ⟩ ⟨ raww = 1.0w = 0.4 ̃̃⟩t −0.06−0.04−0.020.00 ∂ ̃ ̃ ⟨ n n ⟩ ⟨ central⟨diff.w = 1.0w = 0.4 ̃̃⟩t ∂ ̃ ̃ ⟨ n n ⟩ ⟨ central⟨diff.w = 1.0w = 0.4 Figure 9. (Color online) Same as Fig. 8, but for L = 12,∆ τ = 0 .
05, and βt = 20 . Appendix B: Details of DMFT simulations
For our calculations we use the triqs package and the accompanying continuous-time quantum MonteCarlo hybridization expansion solver . We find a crit-ical slow down of the DMFT convergence on the para-magnetic side of the antiferromagnetic transition. Forsome cases we use more than 400 DMFT cycles to obtainreasonable convergence. We subsequently perform someiterations with increased statistics with up to 4 millionsweeps on 80 cores each to obtain data with little noise.We cannot rely on a smoothing algorithm as in the caseof finite-size DQMC data since the kink in the doubleoccupancy would vanish with any smoothing algorithm.We calculate the derivative ∂ ˜ U h n n i by performing for-ward and backward finite differences. If both (forwardand backward differences) are equal within a toleranceof 2%, we take the mean value (i.e., we perform centraldifference); if they are not, we assume the derivative isnot defined at that U [i.e., there is a kink in h n n i ( ˜ U )]. m hypercubic hyperdiamond ̃Ũt −0.14−0.10−0.06 ∂ ̃ ∂ ̃ U ⟨ n n ⟩ βt ⟩ 10βt ⟩ 20 βt ⟩ 40 ̃Ũt βt ⟩ 5βt ⟩ 10 βt ⟩ 20βt ⟩ 40 Figure 10. (Color online) DMFT results for the d = ∞ hyper-cubic (left panels) and hyperdiamond (right panels) lattices.From top to bottom the panels show staggered magnetization m and ∂ ˜ U h n n i . We present DMFT results similar to the DMFT sim-ulations in d = 3 presented in Fig. 5 for U values closeto the antiferromagnetic phase transition for the corre-sponding lattices in d = ∞ , i.e., the hypercubic andhyperdiamond lattices. The case of d = ∞ is interest-ing since DMFT provides an exact solution. The resultsare presented in Fig. 10 in the same way as for d = 3above. We have performed calculations for βt = 10, 20,and 40 and βt = 5, 10, 20, and 40 for the hypercubicand hyperdiamond lattices, respectively. The results arequalitatively very similar to the case of d = 3; that is, wefind discontinuities in ∂ ˜ U h n n i for all temperatures forthe hypercubic case. For the hyperdiamond case we finda discontinuity only at the highest temperature ( βt = 5).The different values for the critical U in d = 3 and d = ∞ can be understood in terms of different bandwidths. Us-ing effective half bandwidths of 2 . t and 2 . t for the hy-percubic and hyperdiamond lattice, respectively, the val-ues of U c /w ( w is the bandwidth) for d = 3 and d = ∞ match nearly perfectly.We provide raw data (Greens function and self-energyof Matsubara frequencies, occupancy, and double oc-cupancy) for the last DMFT iteration for all temper-atures and all four lattices [ cubic, diamond, hypercu-bic ( d = ∞ ), and hyperdiamond ( d = ∞ )] togetherwith a complete set of input parameters on the Zenodoplatform . Appendix C: DQMC results for the cubic lattice andcomparison with DMFT
A DQMC treatment of the Hubbard model on the cu-bic lattice suffers from the scaling of computational timewith the linear lattice size, which is ∝ L and limits thecalculations to L ≤
10. To assess the finite-size scalingwe have performed simulations for βt = 10 at fixed Trot-ter discretization of ∆ τ = 0 .
1. We do not perform anextrapolation to ∆ τ = 0 since our results for the squarelattice and test calculations at βt = 4 for the cubic latticeshow that results for ∆ τ = 0 . L = 4, 6, 8, and 10 for interaction strengthsbetween ˜ U /t = 2 and ˜
U /t = 3 .
78 in steps of 0 .
02 in termsof derivatives of the local and nearest-neighbor chargecorrelators with respect to ˜ U as well as α NN . We smooththe data with w = 0 . ∂ ˜ U h n n i translate into adiscontinuity in α NN , or is it canceled by a respective dis-continuity in ∂ ˜ U h n n i ? Second, what do we learn fromthe comparison of DMFT and DQMC data?To answer the first question, we investigate the finite-size scaling of ∂ ˜ U h n n i and ∂ ˜ U h n n i . As can be seenfrom, e.g., the position of the minimum of ∂ ˜ U h n n i , thefinite-size behavior is nonmonotonous, which makes anextrapolation to L → ∞ impossible based on these data.However, we can observe that the step-like feature visiblein ∂ ˜ U h n n i gets monotonously sharper for larger L . Forthe nearest-neighbor case, no steplike feature is observedfor any L . Consequently, we find steplike features in α NN which get sharper in the same way as they do for thederivative of the local correlator. If we take for grantedthat this behavior holds for L = ∞ , i.e., discontinuitiesexist only for the on-site case and not for the nearest-neighbor case, a discontinuity in ∂ ˜ U h n n i directly trans-lates to one in α NN and thus signals a first-order phasetransition for arbitrarily small nonlocal interactions.The result of the DMFT solution of the cubic latticeis presented in Fig. 11 (a). Although DMFT providesonly an approximate solution for the cubic lattice, it doesgive a result in the thermodynamic limit and can thus beseen as a crude finite-size extrapolation of the DQMCdata. From a superficial inspection of the data the finite-size DQMC data seem to converge against the DMFTresult. A detailed comparison of DMFT with DCA ,dynamical vertex approximation (DΓA) , and DQMC for the cubic Hubbard model suggest that for small inter-action strengths ( U/t . .
5) the DMFT result coincidesrather well with the results obtained with more sophisti-cated methods. This is in line with the finding that thesecond-order correlation energies in d = 3 and d = ∞ donot differ strongly . Finally, results in the thermody-namic limit for the U -dependent double occupancy usingthe numerical linked-cluster expansion show a kink in thedouble occupation in line with our DMFT results .In summary, the DQMC data suggest that, first, a dis- ̃U − . − . (a)̃∂ ̃U < n n > ̃U . . (b)̃∂ ̃U < n n > L = 4L = 6L = 8L = 10DMFT ̃U . . . . . (c)̃α NN Figure 11. (Color online) Results for the derivative of the (a)local and (b) nearest-neighbor charge correlator with respectto ˜ U and (c) α NN as defined in the main text from DQMCsimulation of the cubic lattice at βt = 10 for lattices withlinear sizes L = 4 (solid blue line), L = 6 (dashed yellowline), L = 8 (dash-dotted green line), and L = 10 (dotted redline). Results for the local correlator in the thermodynamiclimit L → ∞ in the DMFT approximation are presented aspurple crosses in (a). continuity in ∂ ˜ U h n n i sufficiently signals a first-orderphase transition at arbitrary small V and, second, theDMFT approximation leads to reasonable results in thethree-dimensional case, especially for small interactionstrengths. Appendix D: Strong coupling calculation of CDWphase
We calculate the CDW phase with a strong-couplingapproach and identify a CDW instability by a negativeFourier component of the Coulomb interaction. For theCDW transition line in the case of the honeycomb latticewe find U = 3 V and U ≈ . V for nearest-neighbor andlong-range interactions, respectively. For the case of thesquare lattice we find U = 4 V and U ≈ . V for nearest-neighbor and long-range interactions, respectively. Forthe interaction strengths of interest ( U/t ∼ U/t ∼ . Appendix E: Proof of thermodynamic inequality
Let H = H + U P i D i − µ P i n i , where D i is the dou-ble occupancy of site i , n i is the occupancy of site i , and H contains all other terms in the Hamiltonian. Then − U and µ can be interpreted as the Lagrange multipliersfixing the average double occupancy, D = 1 /N P i D i ,and the average particle number, n = 1 /N P i n i , re-spectively. The free energy per site of a state with den-sity matrix ρ is given by f ( ρ ) = f ( ρ ) + U D − µn , where f ( ρ ) = 1 /N [ E ( ρ ) − T S ( ρ )]. The thermody-namic ground state ρ minimizes the free energy, so devi-ations δρ from ρ increase the free energy: δf >
0. If weparametrize the density matrix via the double occupancyand the particle number, deviations from the thermody-namic ground state lead to the following changes in thefree energy: δf = ∂f ∂D δD + U δD + 12 ∂ f ∂D δD + ∂f ∂n δn − µδn + 12 ∂ f ∂n δn + ∂ f ∂D∂n δDδn. (E1)The condition that f is at an extremum demands thatthe first-order terms vanish, i.e., ∂f /∂D = − U and ∂f /∂n = + µ . The second-order term can be writtenin matrix form as δf = 12 (cid:0) δD δn (cid:1) ∂ f ∂D ∂ f ∂D∂n∂ f ∂D∂n ∂ f ∂n ! (cid:18) δDδn (cid:19) = 12 (cid:0) δD δn (cid:1) − ∂U∂D − ∂U∂n∂µ∂D ∂µ∂n ! (cid:18) δDδn (cid:19) (E2)Now, the condition δf > < − ∂U∂D , (E3)0 < ∂µ∂n , (E4)0 < ∂U∂n ∂µ∂D − ∂U∂D ∂µ∂n . (E5)Equation (E4) tells us that the compressibility κ = ∂n/∂µ is positive (at constant double occupancy), andEq. (E3) indicates that the double occupancy decreasesas a function of U (at constant density). The symme-try of the matrix implies the Maxwell relation ∂µ/∂D = − ∂U/∂n = A .We consider the relation between ∂D/∂U at constantchemical potential and at constant n . We define the im-plicit function µ ( U, n ) to give the chemical potential cor-responding to U and n , via n ( U, µ ( U, n )) = n . We find ∂n∂U (cid:12)(cid:12)(cid:12)(cid:12) µ + ∂n∂µ ∂µ ( U, n ) ∂U = 0 . (E6)Using this, we obtain ∂D∂U (cid:12)(cid:12)(cid:12)(cid:12) n − ∂D∂U (cid:12)(cid:12)(cid:12)(cid:12) µ = ∂D∂µ ∂µ ( U, n ) ∂U (E6) = − ∂D∂µ ∂n/∂U∂n/∂µ = κ − A − ≥ , (E7)where κ is the compressibility and A is the off-diagonalelement in Eq. (E2). The positivity follows since κ has tobe positive for thermodynamic stability and A appears as a square. Together with Eq. (E3) this gives0 < − ∂D∂U (cid:12)(cid:12)(cid:12)(cid:12) n ≤ − ∂D∂U (cid:12)(cid:12)(cid:12)(cid:12) µ . (E8) ∗ [email protected] J. Hubbard, Proc. R. Soc. A. , 238 (1963). M. C. Gutzwiller, Phys. Rev. Lett. , 159 (1963). J. Kanamori, Prog. Theor. Phys. , 275 (1963). J. Hubbard, Proc. R. Soc. A. , 401 (1964). M. C. Gutzwiller, Phys. Rev. , A923 (1964). J. C. Slater, Physical Review , 538 (1951). J. C. Slater,
The self-consistent field for molecules andsolids , Vol. 4 (McGraw-Hill New York, 1974). N. F. Mott,
Metal-insulator transitions , 2nd ed. (Taylor &Francis, London ; New York, 1990). F. Gebhard,
The Mott metal-insulator transition (Springer, Berlin, 1997). H. Hafermann, E. G. C. P. van Loon, M. I.Katsnelson, A. I. Lichtenstein, and O. Parcollet,Phys. Rev. B , 235105 (2014). E. van Loon, H. Hafermann, A. Lichtenstein, A. Rubtsov,and M. Katsnelson, Phys. Rev. Lett. , 246407 (2014). N. F. Mott, Proc. Phys. Soc., Sect A , 416 (1949). N. F. Mott, Rev. Mod. Phys. , 677 (1968). M. Imada, A. Fujimori, and Y. Tokura,Rev. Mod. Phys. , 1039 (1998). N. Bl¨umer,
Mott-Hubbard Metal-Insulator Transition andOptical Conductivity in High Dimensions , Ph.D. thesis,University of Augsburg (2003). H. Park, K. Haule, and G. Kotliar,Phys. Rev. Lett. , 186403 (2008). M. Balzer, B. Kyung, D. S´en´echal, A.-M. S. Tremblay, andM. Potthoff, Europhys. Lett. , 17002 (2009). J. Merino and O. Gunnarsson,Phys. Rev. B , 245130 (2014). G. Santoro, M. Airoldi, S. Sorella, and E. Tosatti,Phys. Rev. B , 16216 (1993). M. Ulmke, V. Janis, and D. Vollhardt,Phys. Rev. B , 10411 (1995). T. Sch¨afer, F. Geles, D. Rost, G. Rohringer, E. Ar-rigoni, K. Held, N. Bl¨umer, M. Aichhorn, and A. Toschi,Phys. Rev. B , 125109 (2015). F. J. Morin, Phys. Rev. Lett. , 34 (1959). D. B. McWhan and J. P. Remeika,Phys. Rev. B , 3734 (1970). S. A. Shivashankar and J. M. Honig,Phys. Rev. B , 5695 (1983). M. Sch¨uler, M. R¨osner, T. O. Wehling,A. I. Lichtenstein, and M. I. Katsnelson,Phys. Rev. Lett. , 036601 (2013). E. G. C. P. van Loon, M. Sch¨uler, M. I. Katsnelson, andT. O. Wehling, Phys. Rev. B , 165141 (2016). T. Ayral, S. Biermann, and P. Werner,Phys. Rev. B , 125149 (2013). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). M. B. Walker and T. W. Ruijgrok,Physical Review , 513 (1968). D. K. Ghosh, Phys. Rev. Lett. , 330 (1972). T. Koma and H. Tasaki, Phys. Rev. Lett. , 3248 (1992). S. Sorella, Y. Otsuka, and S. Yunoki,Sci. Rep. , 992 (2012). F. F. Assaad and I. F. Herbut,Phys. Rev. X , 031010 (2013). The r j are scaled such that V = V . M. Hohenadler, F. Parisen Toldin, I. F. Herbut, and F. F.Assaad, Phys. Rev. B , 085146 (2014). W. Wu and A.-M. S. Tremblay,Phys. Rev. B , 205128 (2014). P. Buividovich, D. Smith, M. Ulybyshev, and L. vonSmekal, in
PoS (Lattice2016) 244 (2016). H. Terletska, T. Chen, and E. Gull,Phys. Rev. B , 115149 (2017). T. Ayral, S. Biermann, P. Werner, and L. Boehnke,Phys. Rev. B , 245130 (2017). J. E. Hirsch, Phys. Rev. Lett. , 2327 (1984). R. Peierls, Phys. Rev. , 918 (1938). N. N. Bogoliubov, Dokl. Akad. Nauk SSSR , 244(1958). R. P. Feynman,
Statistical Mechanics (Benjamin, ReadingMass., 1972). R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,Phys. Rev. D , 2278 (1981). “QUantum Electron Simulation Toolbox” quest http://quest.ucdavis.edu/ ). D. Rost, E. V. Gorelik, F. Assaad, and N. Bl¨umer,Phys. Rev. B , 155109 (2012). W. Wu, Y.-H. Chen, H.-S. Tao, N.-H. Tong, and W.-M.Liu, Phys. Rev. B , 245102 (2010). P. Fazekas,
Lecture Notes on Electron Correlation and Magnetism. ,Series in Modern Condensed Matter Physics No. v. 5(World Scientific Publishing Company, 1999). M. Jarrell, Phys. Rev. Lett. , 168 (1992). O. Parcollet, M. Ferrero, T. Ayral, H. Hafer-mann, I. Krivenko, L. Messio, and P. Seth,Comput. Phys. Commun. , 398 (2015). P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet,Comput. Phys. Commun. , 274 (2016). G. Rohringer, A. Toschi, A. Katanin, and K. Held,Phys. Rev. Lett. , 256402 (2011). S. Sorella and E. Tosatti, Europhys. Lett. , 699 (1992). Y. ¯Ono, M. Potthoff, and R. Bulla,Phys. Rev. B , 035119 (2003). T. Pruschke and R. Bulla, Eur. Phys. J. B , 217 (2005). J. I. Facio, V. Vildosola, D. J. Garc´ıa, and P. S. Cornaglia,Phys. Rev. B , 085119 (2017). R. T. Scalettar, R. M. Noack, and R. R. P. Singh,Phys. Rev. B , 10502 (1991). M. Sch¨uler, “DQMC data for the Hubbard model...”(2017), doi: 10.5281/zenodo.1118152. A. Savitzky and M. J. E. Golay,Analytical Chemistry , 1627 (1964). M. Sch¨uler, “DMFT data for single band Hubbard model,” (2017), doi: 10.5281/zenodo.1095402. S. Fuchs, E. Gull, L. Pollet, E. Burovski,E. Kozik, T. Pruschke, and M. Troyer,Phys. Rev. Lett. , 030401 (2011). R. Staudt, M. Dzierzawa, and A. Muramatsu,Euro. Phys. J. B , 411 (2000). W. Metzner and D. Vollhardt,Phys. Rev. Lett. , 324 (1989). E. Khatami, Phys. Rev. B94