Thermodynamics of the metal-insulator transition in the extended Hubbard model

In contrast to the Hubbard model, the extended Hubbard model, which additionally accounts for non-local interactions, lacks systemic studies and reference data of thermodynamic properties especially across the metal-insulator transition. We perform such a systematic investigation using an extended variational principle and give a detailed description of how non-local interactions screen local correlations differently in the Fermi-liquid and in the insulator. Furthermore, we study the thermodynamic properties of a discontinuous transition induced by non-local interactions and discuss its relevance in context of metal-insulator transitions in real materials. We substantiate the mechanism of non-local interactions being at least in parts responsible for first-order metal-insulator transition in real materials.


I. INTRODUCTION
The Hubbard model 1-5 of itinerant electrons with a local interaction is arguably the simplest model describing the competition between kinetic energy and interaction. It is of enormous importance for the understanding of strongly correlated materials and it is thus not surprising that much effort has been spent on unveiling its properties through analytical, numerical, and experimental methods. Since the advent of experiments on cold atoms, which closely simulate the Hubbard model, simulations can even be benchmarked against actual experimental data 6 for high temperatures.
Energetics and entropy are, for instance, well known for a broad parameter regime from cold atom experiments 7,8 and complementary numerical simulations using the numerical linked cluster expansion (NLCE) 9 , dynamical cluster approximation (DCA) 10 , determinant quantum Monte Carlo (DQMC) 11 , and recently variational cluster approximation (VCA) 12 and cellular dynamical mean field theory 13 .
Despite its usefulness for understanding the principles of strongly correlated electron physics, the Hubbard model is quite far from describing real materials: Electrons in real materials interact by long-range Coulomb interaction which in the Hubbard model is approximated by neglecting all but on-site interactions.
A more realistic model Hamiltonian in that sense is the extended Hubbard model that incorporates interactions between electrons on different lattice sites. In contrast to the Hubbard model, the properties of the extended Hubbard model are far less well understood. It is known that strong non-local interactions can induce charge density waves [14][15][16][17] , influence superconducting properties [18][19][20][21] , screen local correlations [22][23][24] lead to band-widening effects 25,26 , first-order metalinsulator transitions 27 , and renormalization of Fermi velocities 28,29 . However, systematic studies of thermodynamic properties are rare and reference data for detailed comparisons and benchmarks between methods especially in context of the metal-insulator transition is missing.
Here, we systematically study the half-filled square lattice in the thermodynamic limit (Sec. III). The study includes the transition from Fermi-liquid to insulator in the Slater regime (compare Fig. 1a). We describe two distinct mechanisms of how non-local interactions suppress correlation effects: First, in the Fermi-liquid by effectively reducing the local interaction and second, in an insulating state by effectively increasing the hopping amplitude (Sec. III A). In the presence of non-local interactions, we find that the competition of these two mechanisms drives a first-order metal-insulator transition which we characterize by its latent heat (Sec. III B). We discuss the transition in terms of its relevance in real materials (Sec. III C). We now start by introducing the model and method of approximation in the following Secs. II A and II B.

A. Extended Hubbard model
The extended Hubbard model with non-local interactions reads, where c iσ is the annihilation (creation) operator for electrons on site i and spin σ, t is the nearest-neighbor hopping amplitude, U is the local interaction, and V ij is the nonlocal interaction between electrons at sites i and j. n iσ and n i are the spin-resolved and total occupation operators, respectively. Here, we consider two cases: (a) nearest-neighbor interaction only V 01 ≡ V and (b) long range Coulomb interaction V ij = 1/|r i − r j |. For the latter we fix V 01 by choosing the lattice constant a = 1/V 01 .

B. Variational principle
We investigate the U -V -T phase diagram of the extended Hubbard model by approximating its thermodynamic ground state using the Peierls-Feynman-Bogoliubov variational principle 31-33 with a Hubbard model as the effective system 22 . In former applications only the local interaction was varied, capturing the competition between kinetic and potential energy. Here, we extend that scheme by treating both the hopping and local interaction as variational parameters, which captures the competition between energy and entropy. It turns out, that combining both parameters is necessary to properly study the thermodynamics of the first-order transition. For finite systems, we explore the quality of the variational solution for varying only the local interaction, only the hopping, and both in Ref. 26. The effective Hubbard model with two variational parameters, reads We varyH via the effective hopping parametert and the effective local interactionŨ in order to minimize the free energy functional where the expectation value is performed with the density matrix corresponding toH. FH is the effective free energy. We perform Determinantal Quantum Monte Carlo (DQMC) 34 simulations for the effective Hubbard model on a (t,Ũ ) grid, performing finite-size and Trotter discretization extrapolation, as described in Ref. 27. A two-dimensional Savitzky-Golay filter is used to smooth the data as detailed in Appendix B. We calculate the free energy by integrating from a non-interacting starting point (see Appendix A). For any value of V, we subsequently evaluate the variational free energy function (3) as a function of (t,Ũ ) and determine its minimum.
Here we use a bivariate spline interpolation to lett and U take continuous values. More details of the computational procedure and error estimates can be found in Appendix B.

III. RESULTS
We start by analyzing how non-local interactions influence observables and thermodynamic properties: Fig. 2 shows the temperature dependence of the double occupancy for different non-local interaction strengths in the cases of nearest-neighbor and long-range Coulomb interaction. First, we focus on the double occupancy for the Hubbard model (V = 0). The temperature dependence of the double occupancy, shown in Fig. 2, is rather unusual due to a maximum at finite T . The maximum exists only in the Slater regime at low temperatures, as previously found in Refs. 12 and 30 and can be seen more clearly and in the context of its behavior at larger temperatures in Fig. 1b. We briefly review the nature of this maximum to set the stage for the discussion of effects of non-local interactions.
Generally, disregarding long-range anti-ferromagnetic fluctuations (e.g., in DMFT 35 or in frustrated models 36 ) the double occupancy first decreases with temperature, since the formation of unordered local moments gains entropy over itinerant electrons (discussed extensively in context of adiabatic or Pomeranchuk cooling in optical lattices 37 ). For larger temperatures, the double occupancy eventually increases to meet its non-interacting value in the high-temperature limit. Put together, this gives a minimum close to T /t = 1 for parameters investigated here (compare, e.g., Fig. 11 in Ref. 12 for a larger range of temperatures and interaction strengths). Now, the introduction of long-range anti-ferromagnetic (AF) fluctuations (in more than two dimensions, this corresponds to long range AF ordering) has a drastic effect on the double occupancy: For small interactions (Slater regime), AF fluctuations decrease the potential energy U D by decreasing the double occupancy D. For large interactions (Heisenberg regime) AF fluctuations decrease the kinetic energy through increased hopping probabilities for spins aligned antiparallel 30 . This comes with paying potential energy, i.e., a larger double occupancy D.
In both regimes, AF fluctuations diminish for larger temperatures (roughly the DMFT Néel temperature), such that we expect D to increase with T in the Slater regime and decrease with T in the Heisenberg regime.
Combining AF fluctuations at small temperatures and local moment formation we find two distinct situations: a) In the Heisenberg regime, both effects cooperate such that we still simply find a minimum in D(T ). b) For the Slater regime, AF fluctuations counteract the local moment formation and we find a maximum in D(T ) close to the CDMFT Néel temperature (compare Fig. 1b) and a minimum at larger temperatures. The case b) is exactly the parameter range we study in this work.
To summarize, two mechanisms are responsible for the maximum at V = 0 in Fig. 2. Going to lower temperatures, the double occupancy is lowered to decrease the potential energy (Slater effect). Going to higher temperatures, the double occupancy is lowered to increase the spin entropy (Pomeranchuk effect).
On top of competing formation of local moments and long-ranged AF fluctuations in the Hubbard model, Fig. 2 reveals that non-local interactions introduce multiple effects: First, V generally increases the double occupancy, i.e., decreases local correlation. Second, the position of the maximum is shifted -in the case of nearestneighbor interaction to smaller temperatures and in the case of long range interaction to slightly larger temperatures. Third and most strikingly, a large enough V in both cases introduces a discontinuity. On the high temperature side of this discontinuity, the correlation reducing effect of the non-local interactions is stronger, hinting to a more efficient screening where mobile electrons are present. As highlighted by dotted lines, the discontinuity comes with the introduction of meta-stable states (shown as dotted lines). The possibility of hysteresis and the discontinuous behavior of observables indicate a first-order transition, here.
In the Hubbard model the U dependence of the entropy, S(U ), is connected to D(T ) through a Maxwell nearest neighbor long range relation, ∂S/∂U = −∂D/∂T . Thus, the maximum of S(U ) visible in Fig. 3 has a common origin with maximum in D(T ). Here, it marks the border between the entropy gain through formation of local moments for small U and the Slater regime where U decreases the entropy through AF fluctuations. For even larger U , the entropy passes through a minimum to rise again by local moment formation and finally reach ln 2 in the atomic limit. The introduction of non-local interactions decreases the entropy for most U . This decrease is larger on the large-U side of the discontinuity, i.e., in the regime where AF ordering is favored. For the case of nearestneighbor interactions we observe an increase of the entropy through V for a small regime directly on the low-U side of the discontinuity. For the nearest-neighbor case the maximum is shifted to larger values of U .
For the long-range case, the location of the maximum is hardly influenced by non-local interactions for only observing the region of U values above the discontinuity. However, through the introduction of the discontinuity, we find a maximal S at smaller U than in the case of the pure Hubbard model.
Interestingly, in both cases the discontinuity takes place at rather different values of U . The discontinuity comes with a larger jump in the case of nearest-neighbor interactions.
A. Mechanism of screening local correlations.
So far we have gained insight into the effects non-local interactions have on the properties of the Hubbard model by examining the double occupancy and entropy: we have learned that local correlation effects are reduced and under certain circumstances a first-order transition takes place. Now we study the underlying mechanisms of both effects in terms of the free energy of the effective system.
We essentially find three scenarios for how non-local interactions decrease correlations: 1) In the Fermi-liquid, non-local interactions push the system deeper into the Fermi-liquid by decreasing the effective interaction at fixed hopping. 2) Deep in the correlated regime, nonlocal interactions push the system to effectively lower temperatures by increasing the hopping at fixed interaction strength. 3) In a correlated system close to the metal-insulator transition non-local interactions push the system into the Fermi-liquid regime while increasing the entropy. Under certain circumstances this happens discontinuously. In the following, we will give a detailed analysis of these three regimes, where we from now on focus on the case of nearest-neighbor interaction only, since the case of long-range interaction is qualitatively similar. Figure 4 shows examples of the variational free energy for all mentioned scenarios. A representative example of the Fermi-liquid regime is shown in Fig. 4 a), where we identify a minimum close to thet = t line. Here, a variation ofŨ alone gives basically the same result as varying both parameters: Correlations are reduced by a quenching of the effective interaction. In the strongly correlated regime, shown in Fig. 4 b), we find the opposite case: The minimum is located close to theŨ = U line and the correlations are quenched by increasing the effective hopping.
We observe these two distinct behaviors by tracking how increasing V changes the effective parameters in a T /t-Ũ /t phase diagram, shown in Fig. 5  We identify the different regimes by expanding the free energy functional (Eq. 3) parabolically in V . One can then explicitly show whether the variation oft or that ofŨ leads to a smaller free energy for infinitesimal V as shown in Ref. 26. One finds that the variation oft minimizes the free energy if We finally come to the last and most interesting case of non-local interactions inducing a (discontinuous) metalinsulator transition. We show the corresponding free energy surface in Fig. 4 c), where we find two local minima. Interestingly, both minima are found with similarŨ /t ratio such that the competing configurations differ only in effective temperature T /t. Here the global minimum is realized for the higher effective temperature. By decreasing the non-local interaction V slightly the free energy landscape changes such that the minimum at low temperature (hight) is the global minimum: the non-local interaction induces a first-order phase transition.
We also investigate this third scenario in terms of how the effective parameters in the T /t-Ũ /t phase diagram behave. Both the red and magenta track start for V = 0 at the square markers in the correlated regime such that non-local interactions lead to an increase of the hopping. In the first case we find a smooth transition to a Fermiliquid. In the second case, we find a critical V for which the track jumps at constantŨ /t to larger effective temperature. The entropy, which is color coded in the same figure, increases during the jump with the effective temperature.
This general description of how non-local interactions screen local correlation effects extends findings from GW+DMFT in Ref. 38 which show that the partially (only by non-local interactions) screened local interaction at zero frequency strongly depends on V in the Fermiliquid regime, while it is nearly constant in the Mottinsulator. However, in contrast to our findings the transition from strong to weak screening happens continuously in GW+DMFT.

B. Critical non-local interaction
In the following, we identify the critical non-local interaction V c , which we identify as the smallest V for which we find a discontinuous transition. We investigate the Udependence of V c and the accompanying jump in entropy ∆S (which is proportional to the latent heat) for a fixed hopping of t = 1.26, as shown for T = 0.1 in Fig. 6.
At this temperature, the smallest V c (U ) is ∼ 0.1t, which is significantly smaller than the strong coupling estimate for the charge density wave phase (V U/4). This rules out any kind of masking by a charge-densitywave phase. The jump in entropy is on the order of 10 −3 , on which we base the assessment of observability in experimental setups in the next section. Importantly, finding V c (U ) < U/4 also leads to the possibility of simulating the models containing the prospective discontinuity by numerically exact simulations, e.g., bond-centered auxiliary field Monte Carlo 39 , since V < 4U is the parameter regime, where no sign problem occurs in the square lattice.
We note that for fewer variational degrees of freedom (by fixing eithert = t orŨ = U ) we do not find a firstorder phase transition for the parameters shown in Fig. 6: The increased variational freedom is crucial, here.

C. Comparison to experimentally observed first-order transitions
Experimentally observed discontinuous metalinsulator transitions can be classified into two groups: The first involves transitions where electronic and lattice degrees of freedom participate, e.g., in the case of the metal-insulator transition in VO 2 . These typically come with a large change of entropy (∆S ∼ 1.2 per V atom).
The second group involves purely electronic transitions with a typically smaller jump in entropy. Some examples of the second group are the quantum magnet TiOCl with an entropy change smaller than ∆S = 0.1 per unit cell 40 . In Sm 2 IrIn 8 where an antiferromagnetic transition leads to ∆S = 0.08 per unit cell. 41 Even smaller jumps are observed for an electric-field induced transition in the ferroelectric ceramic Ba 0.73 Sr 0.27 TiO 3 with ∆S = 0.01 per Ti atom 42 and for the metal-insulator transition in NiS 2−x Se x with ∆S between 0.003 and 0.08 per unit cell. 43,44 Since we do not consider any lattice degrees of freedom, we should compare the jumps in entropy found here to the second group. The largest jumps found (which are not masked by a charge density wave) are on the order of ∆S ∼ 0.005, which is well inside the range associated with purely electronic transitions in real materials.

IV. CONCLUSION
We have investigated thermodynamic properties of the metal-insulator transition in the extended Hubbard model, established a reference for future numerical studies, and found discontinuous transitions accompanied with latent heat which is in the range of what is experimentally observed in real materials whenever mainly electronic transitions are involved.
Generally, we have found distinct mechanisms for the suppression of correlation effects in the Fermi-liquid and Mott insulating regime: in the metallic phase the influence of non-local interactions is effectively described by a reduction of the local interaction U . Conversely, in the insulating phase it is described by an increase of the hopping (t). The discontinuous transition finally is described by simultaneously increasing t and decreasing U at a fixed ratio, which comes with an increase of the effective temperature and therefore the entropy.
The regime of parameters for which we find discontinuous Mott transitions rules out masking by charge density waves.
By comparing jumps of the entropy to experimentally observed first-order transition, we are confident, that the mechanism for a first-order transition is of experimental relevance and that it should be a significant contributor to discontinuous metal-insulator transitions in real materials, particularly in situations where lattice distortions are not relevant.
Acknowledgments. We acknowledge financial support of DFG via RTG QM 3 . MIK acknowledges financial support of NWO via Spinoza Prize. The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing computing resources via project number hbp00046 that have contributed to the research results reported in this paper.

Appendix A: Calculation of the free energy
In the DQMC method, we have no direct access to the free energy. However, we do have access to its first derivative the double occupancy: We obtain the free energy at anyŨ by integrating from U 0 = 0 where the free energy in the analytically solvable noninteracting case is given by where ε k are the single-particle energies of the noninteracting system. We calculate the integral over the double occupancy numerically using the composite trapezoidal rule. In a recent work, Walsh et al. 13 extract the entropy of the Hubbard model from the Gibbs-Duhem relation. Their approach is based on an integration over chemical potential, that is, it requires simulations away from half-filling. In contrast, the coupling-constant integration presented here can be performed using only half-filled simulations. For our purpose here, that is a substantial advantage since DQMC suffers from a sign problem away from half-filling.

Appendix B: Error estimation
Performing an error estimate for the results of the variational principle is difficult due to a multitude of error sources: a) statistical (finite sampling) and systematic (finite size and finite Trotter discretization) DQMC errors, b) uncertainties in the evaluation of the free energy, and c) errors from the limits of the variational principle.
In Secs. B 1 and B 2 we explain how we deal with the statistical and systematic DQMC errors. In the remaining sections B 3 and B 4 we discuss errors in the free energy and convergence of sums of charge correlation functions over neighbors throughout the lattice. The most difficult error to assess is the intrinsic error introduced by the limitations of the variational principle. We have performed benchmarks in Ref. 24 for the case of only varying U and we expect qualitatively similar results here. I.e., that the variational principle performs naturally well in describing properties concerning spin fluctuations (e.g., the Mott transition) but bad at properties concerning charge fluctuations (e.g., the charge density wave).

Statistical DQMC Errors
We use the determinant quantum Monte Carlo method (DQMC) 34 implemented in the quest code 45 . We solve Hubbard models in steps of ∆t = 0.02 and ∆Ũ = 0.1 and use the fact that in the case of a finite system all observables behave smoothly withŨ andt. We drastically reduce the statistical error by using a Savitzky-Golay approach 46 : We fit a two-dimensional polynomial to the data in a box of widths wt and wŨ around the point (t 0 ,Ũ 0 ). We use wt = 1.0, wŨ = 1.0, and third order polynomials. We additionally give larger weight to data close to (t 0 ,Ũ 0 ) as to those further away by using a tricubic weighting function (1 − d 3 ) 3 where the distance is defined as d = max{|t −t 0 |/wt; |Ũ −Ũ 0 |/wŨ }. For all parameters we have performed 500 warm up sweeps and 50000 measurement sweeps.

Systematic DQMC errors
Next to statistical errors stemming from finite sampling times in the Monte Carlo procedure which we deal with by smoothing (c.f. Sec. B 1), all observables acquire a systematic error by the Trotter-decoupling for finite imaginary time differences ∆τ . In many cases one achieves qualitatively correct results by choosing small enough values of ∆τ (e.g., by choosing ∆τ ∼ 0.125/U , as suggested in Ref. 47). However, in the case of calculating the free energy (and the entropy) by integrating the double occupancy overŨ we magnify the systematic error (while averaging out the statistical error), which leads to even qualitatively wrong results. We show this for simulations 48 for a larger set ofŨ /t. In Fig. 7 we show the entropy calculated from observables for finite ∆τ and from those which are extrapolated to ∆τ = 0. The entropy for the case of finite ∆τ is negative for large U which breaks basic thermodynamic laws. This is the case, even though both values of ∆τ shown here are smaller than required by the rule of thumb given above for all values ofŨ . Interestingly, this effect worsens for smaller temperatures. Only for the data extrapolated to ∆τ = 0, do we find strictly positive entropy. In the end, we alleviate finite-size and Trotter errors by extrapolating from finite Trotter discretizations of ∆τ = 0.2, 0.1, and 0.05 and linear lattice sizes of L = 8, 10, and 12 for the square lattice following Refs. 49 and 27.
Determinantal Quantum Monte Carlo maps the Hubbard model at finite temperature to an auxiliary d + 1 dimensional Ising model. In this transformation, the temperature is discretized and Suzuki-Trotter discretization errors appear. The observation of negative entropy at ∆τ > 0 in Fig. 7 shows that the discretized system is not thermodynamic, in the sense that basic laws of thermodynamics are violated. Only after the extrapolation to ∆τ = 0 are they restored. Validating these thermodynamic laws serves as a useful check on the extrapolation process and the quality of the Monte Carlo data.
The free energy also provides a useful validation of the Monte Carlo data. The free energy at constant temperature is a function oft andŨ . The nearest-neighbor Green's function and the double occupancy are the observables conjugate to these parameters, i.e., they could be obtained as first derivatives of the free energy. (In practice, we actually obtain the observables from DQMC and then get the free energy by integration.) Thus, the derivative of the double occupancy w.r.t.t is a second derivative of the free energy and by Clairaut's theorem it is equal to the derivative of the Green's function w.r.t. U , modulo prefactors from spin factors. Furthermore, the free energy is a state function and integrating the free energy from t 1 , U 1 to t 2 , U 2 should be independent of the path taken through the t − U plane. These equalities hold for the exact solution of the Hubbard model, for Monte Carlo data it is subject to statistical errors and systematic ∆τ errors and provides a useful estimate of the quality of the data.

Errors related to the free energy
To determine the parametersŨ 0 andt 0 which minimize the free energy functional F ν defined in Eq. 3, we calculate the functional for a fixed mesh ofŨ andt and perform a minimization using interpolating splines of third order. The accuracy of this approach crucially depends on the accuracy of the observables entering F ν , i.e., the double occupancy, nearest-neighbor correlation function and Green's function, and the free energy (from Eq. A2).
We assess the accuracy of the functional by checking if its minimum is located atŨ = U andt = t for the trivial case ofH = H, (no non-local interactions, V = 0). In the upper panels of Fig. 8 we show the deviations from the expected result in terms of the variational parameters, i.e., ∆U =Ũ 0 − U and ∆t =t 0 − t. Both values show considerable deviations from the ideal behavior (∆U = ∆t = 0), which worsens for larger U . The errors are strongly correlated with each other, i.e., we find large positive ∆U where we do so for ∆t. The jumps int corresponding to the first order transition, which can be seen in Fig. 5, are about 0.15 -twice the size of the error estimate shown here.
SinceŨ andt are only variational parameters without any direct physical meaning, it is more interesting to investigate the relative deviations of observables defined as ∆O = (Õ − O)/Õ. We find the relative errors of the double occupancy and nearest-neighbor Green's function to be below one per mille. The error in the free energy is even two orders of magnitude smaller. This roots first in a generally quite flat free energy surface and secondsince the errors inŨ andt are positively correlated -in the weak dependence of the free energy in constantŨ /t direction (c.f. surfaces shown in Fig. 4).
Since the entropy does not only depend onŨ /t but also strongly ont via the effective increase of the temperature T /t, we find more considerable errors, here. The deviations of the entropy from the expected result for the parameters studied in this work are between 1% and 2%.

Finite-size convergence for the Coulomb case
For the case of long-range Coulomb interaction, the calculation of the free energy function (Eq. 3), involves a sum over neighbors: i V 0j n 0 n j . Although we perform a finite size extrapolation of the observables n 0 n j , we can only do this for observables which are actually present in the smallest lattice considered (which are 14 non-equivalent neighbors in the case of an 8 × 8 lattice). We therefore have to check how the sum and thereby the effective parameters converge. We plot the convergence oft andŨ in Fig. 9. Both basically converge for 8 neighbors, such that finite size errors in the case of long-range Coulomb interaction are negligible. The fast convergence relates to the combination of the 1/r behavior of V (r), the exponential decay of spatial correlations for finiteŨ and a checkerboard like oscillation of the sign of ∂/∂Ũ n 0 n j at smallŨ 50 .