# Multigrid Hirsch-Fye quantum Monte Carlo method for dynamical mean-field theory

## Abstract

We present a new algorithm which allows for direct numerically exact solutions within dynamical mean-field theory (DMFT). It is based on the established Hirsch-Fye quantum Monte Carlo (HF-QMC) method. However, the DMFT impurity model is solved not at fixed imaginary-time discretization , but for a range of discretization grids; by extrapolation, unbiased Green functions are obtained in each DMFT iteration. In contrast to conventional HF-QMC, the multigrid algorithm converges to the exact DMFT fixed points. It extends the useful range of , is precise and reliable even in the immediate vicinity of phase transitions and is more efficient, also in comparison to continuous-time methods. Using this algorithm, we show that the spectral weight transfer at the Mott transition has been overestimated in a recent density matrix renormalization group study.

###### pacs:

71.10.Fd, 71.27.+a, 71.30.+h, 02.70.SsMott metal-insulator transitions and other effects of strong electronic correlations are among the most intriguing phenomena in solid state physics (1). They occur when the effective electronic bandwidths of or shell electrons become comparable with the local Coulomb interactions. This most interesting regime is also most challenging: here, perturbative expansions (both at weak and at strong coupling) and effective independent-electron methods such as density functional theory (2) within local density approximation (LDA) break down.

Unfortunately, nonperturbative methods for a direct treatment of correlated lattice electrons are either restricted to one-dimensional (i.e., chain- or ladder-like) systems or suffer from severe finite-size and/or sign problems. However, a significant reduction of complexity in higher-dimensional cases is achieved by the dynamical mean-field theory (DMFT) which maps the electronic lattice problem onto an Anderson impurity model with a self-consistent bath; this mapping becomes exact in the limit of infinite lattice coordination (3). Within the last 15 years, much insight into strongly correlated systems and phenomena has been gained using the DMFT. Many of these studies have relied on the Hirsch-Fye quantum Monte Carlo (HF-QMC) method (4) for solving the DMFT impurity problem at finite temperatures.

The HF-QMC method discretizes the imaginary-time path integral into slices of uniform width and employs a Trotter decoupling; this modified impurity problem is solved via Monte Carlo sampling of a binary Hubbard-Stratonovich field. In principle, arbitrary precision can be achieved using HF-QMC in the combined limit of infinitely many updates of the Hubbard-Stratonovich field and of vanishing discretization ; in this broader sense, HF-QMC is numerically exact. However, practical lower limits for the discretization exist (primarily due to a scaling of the numerical effort with for fixed temperature ), so that raw HF-QMC results usually contain systematic Trotter errors which are much larger than the statistical Monte Carlo errors. This implies, in the DMFT context, that all observables and phase boundaries depend on the auxiliary parameter chosen in the self-consistency cycle. While high accuracy and efficiency have been demonstrated for a posteriori extrapolations of certain static observables (5); (6), these procedures require special care and experience and may fail close to phase transitions. Up to recently, the Trotter bias of HF-QMC derived spectral functions could not be reduced at all.

For a long time, Trotter errors could only be avoided using fundamentally different impurity solvers that either introduce a logarithmic frequency discretization directly (numerical renormalization group, NRG) or represent impurity plus bath as a finite cluster or chain [exact diagonalization (ED), density-matrix renormalization group (DMRG)]. Consequently, their results, too, contain systematic finite-size/discretization errors; in particular, these alternatives to QMC yield continuous spectra only after numerical broadening. Only very recently, two new quantum Monte Carlo impurity solvers have become available (7); (8) which are numerically exact in the stricter sense that their results are correct within statistical error bars, without the need for explicit extrapolations. These methods avoid the imaginary-time discretization of HF-QMC and rely, instead, on perturbative expansions in the interaction and hybridization, respectively, which are statistically sampled to arbitrary order. However, the continuous-time quantum Monte Carlo (CT-QMC) methods are often less efficient than HF-QMC (6).

In this Letter, we propose a new algorithm for solving the DMFT self-consistency equations in quasi continuous imaginary time, based on a multigrid implementation of the HF-QMC method. This implementation removes many of the limitations of HF-QMC; in particular, it yields precise and reliable results even at phase boundaries. As a first application, we test DMRG predictions (9) of the spectral weight transfer at the Mott transition – a study that would not have been possible using conventional HF-QMC.

Multigrid HF-QMC method – Before we formulate the new multigrid approach, let us discuss the flow diagram Fig. 1(a) of the conventional HF-QMC method.

Starting with some guess for the self-energy, the lattice Green function is obtained via the lattice Dyson equation (upper box, here written for the 1-band Hubbard model); as a next step, the impurity Dyson equation (middle box) is used to determine the bath Green function which defines the DMFT impurity problem (lower box). Conventionally, this is solved using HF-QMC at a fixed discretization; the resulting Green function closes the DMFT cycle via the impurity Dyson equation. Obviously, the Trotter error in the impurity solver affects the whole DMFT cycle; thus, at self-consistency, , , and all observables deviate from their physical values, ideally with corrections in powers of . Only then, numerically exact estimates of observables can be extrapolated a posteriori from results of independent HF-QMC DMFT runs, as illustrated in Fig. 1(b).

The new multigrid method, visualized in Fig. 1(c), splits the impurity-solving part of the DMFT self-consistency cycle into two steps: first, the impurity problem is solved – at fixed bath Green function – for a suitable range of discretizations using HF-QMC. In a second step, a numerically exact estimate of the true Green function of the given impurity problem is obtained by extrapolation of the multiple resulting Green functions . This highly nontrivial task is accomplished using a scheme developed recently (in the context of a posteriori extrapolation) (10): (i) each of the discrete HF-QMC estimates is interpolated with the help of a suitably chosen reference model onto a common fine grid; (ii) quadratic (in ) least-squares extrapolations are performed on this grid. In practice, the raw Green functions are averaged over 10–20 impurity solutions. Thus, for about 10 discretizations , the multigrid method parallelizes to about 100 CPU cores. For optimal accuracy, a hierarchy of discretization scales and high-frequency cut-offs needs to be maintained (11); details will be presented elsewhere (12). Since the multigrid algorithm closes the DMFT self-consistently with an impurity Green function that is (for correctly chosen parameters, see below) unbiased, all state variables , , and converge to their numerically exact forms.

This fundamental advantage of the multigrid HF-QMC method is illustrated in the scheme Fig. 2:

Here, the solid line depicts the exact Ginzburg-Landau (GL) free energy in the multi-dimensional space of bath Green functions (for fixed physical parameters, , ); its stationary points mark the DMFT fixed points (13). However, conventional HF-QMC implementations converge to the stationary points (squares, diamonds) of generalized GL functionals (dotted lines) with additional dependence on . These fixed points can vary irregularly or even discontinuously as a function of which complicates or prevents extrapolations . In contrast, the multigrid HF-QMC procedure drives the DMFT iteration to the numerically exact fixed point: at self-consistency, each finite- HF-QMC evaluation is performed for the bath Green function at which is stationary (full circle) while the generalized functionals are not (empty circles).

Benchmark results – For a quantitative discussion, let us now consider the one-band Hubbard model with semi-elliptic density of states (bandwidth ) at low temperatures in the strongly correlated regime. Specifically, we will concentrate on the double occupancy , an observable which is well-defined for the impurity model irrespective of DMFT self-consistency and which is best computed at fixed in both algorithms. Figure 3

shows results for the coexisting metallic and insulating phases versus squared discretization at and interaction . Conventional HF-QMC (diamonds) finds an insulating solution only for relatively small . While a metallic solution can be stabilized at arbitrarily large , the dependence of is highly nonuniform which restricts the useful range for extrapolations (dashed lines), again, to small values . In contrast, the multigrid HF-QMC raw data (circles) shows regular dependence even at large discretizations which allows for very precise extrapolations (solid lines) based on cheap large- HF-QMC data (14). Evidently, the multigrid method is superior.

However, such extrapolations, as well as all observable estimates that can be derived from the Green function and self-energy without explicit extrapolation, are only reliable if the multigrid algorithm has really converged to the exact DMFT fixed point. Thus, we have to check for any bias that might survive from the intrinsic parameters, most notably the range used in the internal Green function extrapolation. In the inset of Fig. 4,

multigrid HF-QMC results (at , i.e, in the metallic phase) are plotted for various grids; at this scale, both the raw data and the fit curves fall on top of each other. For a fine-scale analysis, the approximate leading Trotter errors have been subtracted in the main panel of Fig. 4. Even for this extremely precise raw data (with sweeps for each in 10 iterations after convergency), no impact of the grid range can be detected for . A slight bias of the order is visible for the grid range ; even the coarsest range of discretizations shifts only by about . Similar studies indicate even slightly smaller grid range effects in the insulating phase (at ) (12). Thus, the multigrid HF-QMC method can be considered numerically exact for grid ranges up to in almost all cases; finer grids or a detailed analysis are only needed when extreme precision is sought. This should be contrasted with the conventional HF-QMC method for which biases in Green functions and self-energies should remain detectable down to – if reasonable statistical accuracy could be maintained.

Spectral weight transfer at the Mott transition – While many aspects of the Mott metal-insulator transition are well-established by now, including the phase diagram of the half-filled frustrated one-band model within DMFT, some fundamental questions have remained open. In particular, it is still not clear how exactly the spectra change across the first-order Mott transition, where the quasiparticle peak disappears. Recently, a new scenario has been suggested on the basis of dynamical density matrix renormalization group (DDMRG) calculations (9). As seen in the inset of Fig. 5,

this study finds sharp features at the inner edges of the Hubbard bands in the strongly correlated metallic phase (solid lines), which are separated from the quasiparticle peak by quite broad pseudo gaps; in contrast, the Hubbard bands appear rather featureless in the insulating phase (dashed lines) for the same interaction , with inner edges that are shifted to much smaller frequencies. Obviously this behavior, which has been interpreted in terms of collective magnetic excitations in the metal (9), should be verified. One may also ask whether it extends to finite temperatures, i.e., governs the spectral weight transfer at the Mott transition, which, for , occurs at .

Estimates of the spectral weight transfer, the difference of metallic and insulating spectra, are shown in the main panel of Fig. 5. At this scale, the DDMRG data (solid line) clearly resolves the features discussed above, in particular the Hubbard subpeak at . The corresponding multigrid HF-QMC results (broken lines) have been obtained, for minimal bias, via Padé analytic continuation of Matsubara difference Green functions. They are much smoother even at very low : apart from the quasiparticle peak, with a reduced weight, only a shallow negative region at is visible, which is evidence against the existence of subpeaks. We obviously cannot exclude that this smoothness is, at least partially, an artifact of the ill-conditioned analytic continuation; it is also not clear whether the change of shape at the lowest temperature , with a closer approach of the DDMRG line shape, is significant. In contrast, the QMC difference Green functions shown in Fig. 6 as broken lines (basis for the spectral data of Fig. 5) are precise within linewidth. In the low- regime that we are focusing on, even the extrapolation to is reliable (see inset of Fig. 6). However, the transformed DDMRG data (thin solid line in main panel) deviates markedly from this QMC ground state result (thick solid line, with a precision of about linewidth): the DDMRG overestimates the differences between metal and insulator, in this measure, by about , with an even larger discrepancy compared to the finite-temperature result at , where the metal-insulator transition takes place. We conclude that the DDMRG does not capture the spectral-weight transfer at the Mott transition quantitatively correctly; possible explanations include the bias towards metallic/insulating solutions for odd/even chains in the DMRG calculation and, in particular, the deconvolution technique. In terms of spectra, a large part of the DDMRG error can be traced back to its overestimation of the quasiparticle peak (since the QMC data of Fig. 5 is particularly reliable in this range); similar mechanisms might have generated the Hubbard-band subpeaks.

Discussion – We have formulated a new algorithm for solving Hubbard-type models within DMFT, based on a multigrid implementation of the HF-QMC method with internal Green function extrapolation. This technique directly yields numerically exact results even at very low temperatures, similarly to recently developed continuous-time QMC methods. However, our algorithm retains the advantages of HF-QMC, i.e., generates Green functions with minimal fluctuations, now at about halved temperatures for similar precision and effort. The method extends to the general multi-band case, with enormous potential, e.g., in the context of ab initio LDA+DMFT calculations. While it is too early to judge which of the three (quasi) continuous-time methods will ultimately prevail, our multigrid method clearly supersedes conventional HF-QMC for most applications.

Using this method, we have found inconsistencies in the DDMRG scenario of the spectral weight transfer at the Mott transition. A definite judgement whether the DDMRG results are, at least, qualitatively correct, will require additional computing ressources and detailed analyses of analytical continuation procedures.

Stimulating discussions with P.G.J. van Dongen and support by the DFG within the Collaborative Research Centre SFB/TR 49 are gratefully acknowledged.

### References

- See G. Kotliar and D. Vollhardt, Physics Today 57, 53 (2004) and references therein.
- See R. Martin, Electronic structure, Cambridge University Press (2004) and references therein.
- A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986); M. Jarrell, Phys. Rev. Lett. 69, 168 (1992).
- N. Blümer and E. Kalinowski, Phys. Rev. B, 71, 195102 (2005); Physica B 359, 648 (2005).
- N. Blümer, Phys. Rev. B 76, 205120 (2007).
- A.N. Rubtsov, V.V. Savkin, and A.I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
- P. Werner, A. Comanac, L. de Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
- M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 72, 113110 (2005).
- N. Blümer, arXiv:0712.1290v1 (unpublished).
- E.g.: 100–400 time slices in HF-QMC, 1000 Matsubara frequencies for state variables, 10000 interpolated time slices, and 40000 Matsubaras for the reference model.
- N. Blümer, in preparation.
- G. Kotliar, E. Lange, and M. J. Rozenberg, Phys. Rev. Lett. 84, 5180 (2000).
- We have verified that the assumed extrapolation law remains accurate at least up to , when is increased from 3 to about 8.