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.
Numerous studies of the Hubbard model  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) , 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 , 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) . 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)  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 , 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.
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  the Trotter errors can be eliminated with high precision by piecewise extrapolations of (circles) . 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.
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 
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.
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”  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 . 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 ).
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 ), 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  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.
- Hubbard J 1959 \PRL3 77–78
- Schollwöck U 2005 \RMP77 259
- Metzner W and Vollhardt D 1989 \PRL62 324–327
- Georges A, Kotliar G, Krauth W and Rozenberg M J 1996 \RMP68 13–125
- Lee P A and Wen X-G 2006 \RMP78 17
- Armitage N, Fournier P and Greene R 2010 \RMP82 2421
- Rost D, Gorelik E V, Assaad F F and Blümer N 2012 Phys. Rev. B 86 155109
- Mott N F 1968 \RMP40 677
- Gebhard F 1997 The Mott Metal-Insulator Transition (Springer Berlin)
- Staudt R, Dzierzawa M and Muramatsu A 2000 Eur. Phys. J. B 17 411
- Park H, Haule K and Kotliar G 2008 \PRL101 186403
- Gorelik E V, Rost D, Paiva T, Scalettar R, Klümper A and Blümer N 2012 Phys. Rev. A 85 061602
- Chang C-C, Scalettar R T, Gorelik E V and Blümer N 2013 Phys. Rev. B 88 195121
- Blankenbecler R, Scalapino D J, and Sugar R L 1981 Phys. Rev. D 24 2278
- Fuchs S, Gull E, Troyer M, Jarrell M and Pruschke T 2011 Phys. Rev. B 83 235113
- Blümer N 2007 Phys. Rev. B 76 205120
- Blümer N 2008 Preprint arXiv:08011222
- Gorelik E V and Blümer N 2009 Phys. Rev. A 80 051602(R)
- Rost D, Assaad F F and Blümer N 2013 Phys. Rev. E 87 053305
- This nontrivial task is achieved using reference models and spline interpolation of differences between measured and reference Green functions [30, 27, 16, 17, 7].
- Fetter A L and Walecka J D 1971 Quantum Theory of Many-Particle Systems McGraw-Hill
- 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 , in analogy to orbital-selective Mott phases [27, 28, 29] observed in systems with inequivalent bands.
- Toschi A, Katanin A A and Held K 2007 Phys. Rev. B 75 045118
- Held K, Katanin A A and Toschi A 2008 Prog. Theor. Phys. Suppl. 176 117
- 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
- Ferrero M, Cornaglia P, De Leo L, Parcollet O, Kotliar G and Georges A 2009 Phys. Rev. B 80 064501
- Knecht C, Blümer N and van Dongen P G V 2005 Phys. Rev. B 72 081103(R)
- Jakobi E, Blümer N and van Dongen P G V 2009 Phys. Rev. B 80 115109
- Jakobi E, Blümer N and van Dongen P G V 2013 Phys. Rev. B 87 205135
- Blümer N and Kalinowski E 2005 Phys. Rev. B 71 195102