Deciding the fate of the false Mott transition in two dimensions by exact quantum Monte Carlo methods

Deciding the fate of the false Mott transition in two dimensions by exact quantum Monte Carlo methods


We present an algorithm for the computation of unbiased Green functions and self-energies for quantum lattice models, free from systematic errors and valid in the thermodynamic limit. The method combines direct lattice simulations using the Blankenbecler-Scalapino-Sugar quantum Monte Carlo (BSS-QMC) approach with controlled multigrid extrapolation techniques. We show that the half-filled Hubbard model is insulating at low temperatures even in the weak-coupling regime; the previously claimed Mott transition at intermediate coupling does not exist.

1 Introduction

Numerous studies of the Hubbard model [1] have greatly enhanced our understanding of strongly correlated electron systems within the last four decades. Specialized methods, namely the semi-analytic Bethe ansatz and the density matrix renormalization group (DMRG) [2], yield reliable high-precision results (only) in the case of one spatial dimension. Conversely, the dynamical mean-field theory (DMFT) [3, 4] provides deep insight in the limit of high dimensionality. However, the intermediate regime of two (and three) dimensions is less well understood, e.g., with respect to pseudogap physics [5, 6, 7] and high- superconductivity. In this paper, we will address another open question: the nature of the Mott metal-insulator transition (MIT) [8, 9] of the half-filled Hubbard model in two dimensions ().

As indicated by shaded regions in Fig. 2, the Hubbard model is insulating at half filling (1 electron per site) and sufficiently strong coupling, in any dimension, while it is metallic at weak coupling – as long as the translational symmetry is not broken, e.g., by magnetic order. For , one finds antiferromagnetic (AF) ordering [10], which implies insulating behavior, below the Néel temperature (gray dash-dotted line). Such long-range order is ruled out by the Mermin-Wagner theorem in (at temperature ). Single-site DMFT predicts a MIT with critical interactions (pink dotted line). This transition line is shifted to weaker interactions within cluster DMFT (CDMFT), which includes short-range correlations on small clusters (violet solid line) [11]. However, the low- metallic phase implied by the CDMFT result is incompatible with several earlier studies highlighting effects of “short-range antiferromagnetism” [12, 7, 13]. To clarify the situation, we apply the Blankenbecler-Scalapino-Sugar quantum Monte Carlo (BSS-QMC) [14] approach with controlled multigrid extrapolation techniques to selected points in the phase diagram (green crosses). We will show that these points are separated by a MIT line (or narrow crossover) and conclude that the Mott physics is quite similar to the case [15], in spite of the Mermin-Wagner theorem.

2 Model and Methods

The single-band Hubbard Hamiltonian is given by


with number operators , next-neighbor hopping , and local interactions . In this paper, we set and consider finite-size clusters with periodic boundary conditions.

Figure 1: Schematic representation of the MIT in the Hubbard model at half filling on a square lattice. Green crosses: parameters selected for our unbiased QMC study. Gray dash-dotted line: Néel temperature for the cubic lattice (rescaled).


Figure 2: Extrapolation of BSS-QMC Green functions to . Analytic continuation (MEM, black line) regularizes the results of pointwise extrapolations (symbols) and yields stable results (gray line) even in the case of noisy raw data (thin color lines).

The BSS-QMC algorithm is based on a Trotter-Suzuki decomposition of the partition function


where () corresponds to the interaction (kinetic) term in (1), is the inverse temperature () and denotes the number of time slices (). A discrete Hubbard-Stratonovich transformation simplifies the interaction term, which is quartic in the fermionic operators, to a quadratic form and a coupling to an auxiliary Ising field . As usual for quadratic Hamiltonians, the fermionic degrees of freedom can be integrated out and the partition function is expressed by determinants of matrices:


The sum in (3), and finally the computation of all observables, is performed by Monte Carlo methods. The main results of the BSS-QMC algorithm are Green functions on a discrete imaginary-time grid for each pair of real-space coordinates on the given lattice. The involved matrix operations for the calculation of the determinants in (3) lead to a scaling of , which restricts the algorithm to relatively small clusters ().

To arrive at reliable conclusions regarding the Mott transition using BSS-QMC, one has first to eliminate all systematic errors (finite-size as well as Trotter errors) and, secondly, calculate the self-energy on the Matsubara axis, , from the imaginary-time Green function in a stable manner. These steps are described in the following.

3 Elimination of systematic errors

As the elimination of Trotter errors is well established [16, 17, 18, 19], we will sketch the scheme only briefly and focus on the problem specific and physically more relevant finite-size analysis and present a stable method to calculate unbiased Green functions, even if the raw data are noisy.

Trotter bias – The raw BSS-QMC Green functions (thick colored lines in Fig. 2) do not only live on different imaginary-time grids , depending on the chosen discretization , but are also shifted with respect to each other (and to the exact solution). After aligning them on a common fine grid [20] the Trotter errors can be eliminated with high precision by piecewise extrapolations of (circles) [19]. However, some fluctuations remain inevitably, especially when using noisy raw data (using a factor of 10 less QMC sweeps; thin colored lines). These can be greatly reduced by regularization via a maximum entropy method (MEM), which computes corresponding spectral functions via the inversion of


Note that the intermediate spectra are used only for producing continuous and smooth Green functions via (4); in this combination, the procedure is stable. After this step, the results based on high- or low-precision data (thick black line vs. gray line) are hardly distinguishable and can both be used for stable computations of self-energies.

Figure 3: (a) Influence of statistical noise in the raw data on the self-energy and stabilization by MEM regularization. Thin lines: noisy data, thick lines and symbols: high precision results. (b) Impact of Trotter discretization on the Matsubara self-energy.

Self-energy – These quasi-continuous can be reliably Fourier transformed to the Matsubara axis:


A quantity of great interest for the analysis of the Mott transition is the imaginary part of the momentum-resolved self-energy, which is connected to the BSS-QMC Green function and the non-interacting Green function via a Dyson equation [21]


where is the chemical potential and the dispersion of the noninteracting problem. As shown in Fig. 3 (b), the MEM regularizing procedure leads to estimates (crosses and black line) of the imaginary part of the self-energy (the real part vanishes at due to particle-hole symmetry) which are smooth even at higher frequencies. In contrast, the direct results (green circles) fluctuate significantly even when based on high-precision raw data (and even more so for noisy raw data; thin dashed green lines). Note that both methods agree very well at the lowest Matsubara frequencies, where the QMC predictions are most reliable.

In Fig. 3 (b) the influence of the Trotter discretization is shown. The Trotter biased self-energies (broken colored lines) deviate from the exact result (black solid line) for large values of , in particular at the lowest frequencies. However, in the parameter regime explored in this study, we can afford small enough values of so that the elimination of Trotter errors is less essential than the finite-size extrapolation to be discussed below.

Figure 4: Finite-size scaling of Matsubara self-energy at , . Finite-size BSS-QMC data (open symbols and broken colored lines), extrapolated BSS-QMC results in the thermodynamic limit (circles and black bold solid line), and DA data (gray line) versus Matsubara frequency at (a) and (c); also shown are momentum-independent single-site DMFT results (thin black lines). (b)+(d) Finite-size BSS-QMC (symbols) data for the first three Matsubara frequencies versus inverse system size plus extrapolations in linear order in (thin lines) and quadratic order (thick lines).

Finite-size extrapolation of the self-energy – Using the methods described above, unbiased and stable estimates of the Matsubara self-energy have been obtained for square clusters with sizes ranging from to (with a computational effort that varies as , i.e., by a factor of in this range). Results at the “anti-nodal” [22] momentum point [throughout the paper, unit lattice spacing is assumed] are shown (colored symbols and broken lines) in Fig. 4(a) for the low temperature and in Fig. 4(c) for the elevated temperature , respectively. Evidently, the finite-size effects are enormous: while the smallest systems (, triangles) show insulating behavior, i.e., a strong enhancement of towards small frequencies, at both temperatures, this tendency is reduced with increasing at [Fig. 4(a)] and completely eliminated at [Fig. 4(c)]. Obviously, careful extrapolations are needed for reliable predictions in the thermodynamic limit.

These are, indeed, possible, as shown in Fig. 4(b) for and in Fig. 4(d) for , respectively, for the three lowest Matsubara frequencies (where the finite-size effects are largest): as a function of , the finite-size results (symbols) can be fitted with second-order polynomials (thick lines) with high precision; as these fit functions have small curvatures, linear extrapolations (thin lines) to the thermodynamic limit (i.e., ) deviate only slightly from the quadratic ones. We use these deviations as error bars and the arithmetic average of both extrapolation procedures as final result, as indicated for by a black symbol in Fig. 4(b).

4 Mott physics and correlation lengths

These final QMC results [black symbols and lines in Fig. 4(a) and Fig. 4(c)] show that the character of the system changes drastically between the selected phase points (crosses in Fig. 2): while the QMC self-energy indicates a metallic phase at , very similar to the DMFT solution [thin black line in Fig. 4(c)], the low-temperature phase is clearly insulating - and completely unlike the paramagnetic DMFT solution [thin black line in Fig. 4(a)]. In contrast, the QMC results are strikingly similar to the DA predictions (gray dash-dotted lines), especially at low . The dynamical vertex approximation () [23, 24] is a diagrammatic extension of DMFT which works directly in the thermodynamic limit (without any need for finite-size extrapolations); for the full set of results see [25]. QMC and DA together show convincingly that the suspected phase transition (dotted line in Fig. 2) is, indeed, present.

The question remains how exactly this metal-insulator transition (or crossover) is related to magnetic properties. After all, long-range AF order is excluded, in , by the Mermin-Wagner theorem. However, as shown in Fig. 5(a), the spin correlations decay much more slowly at (circles) than at (squares). More precisely, the correlation length at the elevated temperature can be estimated (using QMC data) to while it is clearly too large () at to allow precise determination using QMC. For this observable, DA predictions, shown in Fig. 5(b) and Fig. 5(c) are more reliable, yielding correlation lengths of at the higher temperature, in excellent agreement with QMC, and at , also consistent with QMC. Using DA we could also show that the temperature dependence of changes its character at the Mott transition (see [25]).

Figure 5: Real-space spin correlation function. (a) The correlation length for (red squares) from BSS-QMC calculations for a lattice is estimated by an exponential fit (red solid line). For (green circles) is much larger than the extend of the lattice. Right panel: The corresponding estimations from DA data [25] for (b) and (c) confirm the BSS-QMC results.

Conclusion and acknowledgements

In spite of its simplifications, the Hubbard model remains a challenging target for theorists, especially in the interesting case of two spatial dimensions. All available methods have certain limitations, many of which are of particular concern in the nonperturbative intermediate-coupling regime at low temperatures and in the presence of longer-range correlations. Therefore, a single method can hardly yield authoritative results; some predictions may even be completely off (such as the CDMFT transition line in Fig. 2).

In the study presented in this work (with focus on the QMC methodology; for the full story, see [25]), we have applied two methods (an unbiased variant of BSS-QMC as well as the dynamical vertex approximation) with completely different characteristics. Only the near-perfect agreement between both sets of results makes the predictions truly authoritative (and validates technical choices on either side). We have settled one important question at half filling, namely the character of the Mott metal-insulator transition on the square lattice: it is driven by exponentially long-ranged [25] antiferromagnetic correlations, which act similarly to the AF long-range order in the cubic case, and is not connected to a quantum-critical point at .

The same methodology may be useful in other parameter ranges, e.g., for frustrated or anisotropic Hubbard models, doped systems, or multi-band models. It might also be worthwhile to compute transport properties, within bubble approximation or beyond, in order to observe the metal-insulator transition even more directly (than via the self-energy).

We acknowledge support from the research unit FOR 1346 of the German Research Foundation (DFG) and the graduate school GSC 266.



  1. Hubbard J 1959 \PRL3 77–78
  2. Schollwöck U 2005 \RMP77 259
  3. Metzner W and Vollhardt D 1989 \PRL62 324–327
  4. Georges A, Kotliar G, Krauth W and Rozenberg M J 1996 \RMP68 13–125
  5. Lee P A and Wen X-G 2006 \RMP78 17
  6. Armitage N, Fournier P and Greene R 2010 \RMP82 2421
  7. Rost D, Gorelik E V, Assaad F F and Blümer N 2012 Phys. Rev. B 86 155109
  8. Mott N F 1968 \RMP40 677
  9. Gebhard F 1997 The Mott Metal-Insulator Transition (Springer Berlin)
  10. Staudt R, Dzierzawa M and Muramatsu A 2000 Eur. Phys. J. B 17 411
  11. Park H, Haule K and Kotliar G 2008 \PRL101 186403
  12. Gorelik E V, Rost D, Paiva T, Scalettar R, Klümper A and Blümer N 2012 Phys. Rev. A 85 061602
  13. Chang C-C, Scalettar R T, Gorelik E V and Blümer N 2013 Phys. Rev. B 88 195121
  14. Blankenbecler R, Scalapino D J, and Sugar R L 1981 Phys. Rev. D 24 2278
  15. Fuchs S, Gull E, Troyer M, Jarrell M and Pruschke T 2011 Phys. Rev. B 83 235113
  16. Blümer N 2007 Phys. Rev. B 76 205120
  17. Blümer N 2008 Preprint arXiv:08011222
  18. Gorelik E V and Blümer N 2009 Phys. Rev. A 80 051602(R)
  19. Rost D, Assaad F F and Blümer N 2013 Phys. Rev. E 87 053305
  20. This nontrivial task is achieved using reference models and spline interpolation of differences between measured and reference Green functions [30, 27, 16, 17, 7].
  21. Fetter A L and Walecka J D 1971 Quantum Theory of Many-Particle Systems McGraw-Hill
  22. In the context of high- and pseudogap physics, one distinguishes “nodal” from “anti-nodal” points on the noninteracting Fermi surface. The observed dichotomy in the pseudogap phase can be viewed as “momentum-selective” Mott physics [26], in analogy to orbital-selective Mott phases [27, 28, 29] observed in systems with inequivalent bands.
  23. Toschi A, Katanin A A and Held K 2007 Phys. Rev. B 75 045118
  24. Held K, Katanin A A and Toschi A 2008 Prog. Theor. Phys. Suppl. 176 117
  25. Schäfer T, Geles F, Rost D, Rohringer G, Arrigoni E, Held K, Blümer N, Aichhorn M and Toschi A 2014 Phys. Rev. B 91 125109
  26. Ferrero M, Cornaglia P, De Leo L, Parcollet O, Kotliar G and Georges A 2009 Phys. Rev. B 80 064501
  27. Knecht C, Blümer N and van Dongen P G V 2005 Phys. Rev. B 72 081103(R)
  28. Jakobi E, Blümer N and van Dongen P G V 2009 Phys. Rev. B 80 115109
  29. Jakobi E, Blümer N and van Dongen P G V 2013 Phys. Rev. B 87 205135
  30. Blümer N and Kalinowski E 2005 Phys. Rev. B 71 195102
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description