# Numerically exact Green functions from Hirsch-Fye quantum Monte Carlo simulations

## Abstract

We present a new method for extracting numerically exact imaginary-time Green functions from standard Hirsch-Fye quantum Monte Carlo (HF-QMC) simulations within dynamical mean-field theory. By analytic continuation, angular resolved spectra are obtained without the discretization bias previously associated with HF-QMC results. The method is shown to be accurate even at very low temperatures ( for bandwidth ) in the strongly correlated regime.

###### pacs:

71.10.Fd, 71.27.+a, 78.20.Bh, 02.70.SsPhotoemission spectroscopy (PES) and related techniques (inverse PES, X-ray absorption spectroscopy) are among the most useful experimental probes of electronic properties of solids (1). Already the angle- and spin-integrated variants can characterize materials as weakly or strongly correlated metals or as insulators. In addition, angle-resolved photoemission spectroscopy (ARPES) gives a nearly direct view on the occupied part of the electronic band structure. However, the interpretation of such experimental data and the separation of surface from bulk contributions is by no means trivial. Thus, both a reliable analysis of the experiments and an understanding of the underlying physics require comparisons with theoretical predictions. Often, good agreement is found with spectra obtained within density functional theory (DFT), e.g., within the local density approximation (LDA) (2). For strongly correlated materials, however, the LDA and/or the interpretation of the DFT parameters as many-body dispersion break down; then, true many-body approaches, usually based on Hubbard-type models, are required.

A nonperturbative treatment of Hubbard-type models for correlated electron systems is possible by iterative solution of the DMFT self-consistency equations (3) using quantum Monte Carlo (QMC) methods. In the Hirsch-Fye QMC algorithm (4), the imaginary-time path integral is discretized into time slices of uniform width ; a Hubbard-Stratonovich (HS) transformation replaces the electron-electron interaction at each time step by a binary auxiliary field which is sampled by standard Markov Monte Carlo techniques. After convergency, estimates of the local imaginary-time Green function are obtained (only) on the grid with (). For the computation of spectral functions, this data has to be continued to the real axis, typically using maximum entropy methods (MEM) (5); these regularize the ill-conditioned inversion of the equation

(1) |

by finding a spectrum that is as smooth (or similar to a chosen default model) as possible within the constraints set by the data and its error bars. If needed, the real part of the local Green function is then obtained by Kramers-Kronig transformation of Im. Finally, momentum resolved spectra corresponding to ARPES measurements may be obtained via the real-frequency self-energy .

While impressive results have been obtained using the QMC/MEM procedure as outlined above, e.g., in the context of “kink” anomalies in ARPES spectra (6), it has one fundamental problem: the Hirsch-Fye QMC estimates of the imaginary-time Green function contain systematic Trotter errors which are typically much larger than the statistical errors, often by orders of magnitude. However, the MEM takes only the statistical errors into account, partially by quite elaborate formalism, while the larger systematic errors in the data are neglected. Consequently, all resulting spectra are biased, to an unknown extent, by the Trotter discretization.

In this Letter, a method is proposed which essentially eliminates these problems: using a novel extrapolation scheme, we extract continuous estimates of the imaginary-time Green function without significant Trotter error from conventional HF-QMC data; these are then proper bases for analytic continuation techniques. The method works well even at very low temperatures and for comparatively coarse imaginary-time discretizations. These properties make it very attractive for future calculations of spectra, e.g., in ab initio LDA+DMFT studies; they also establish that – contrary to common belief – reliable HF-QMC results do not depend on a good resolution of the rapid initial decay of .

Raw HF-QMC results – In the following, the method will be defined and illustrated using a quite ambitious example: the half-filled single-band Hubbard model (semi-elliptic density of states with bandwidth ) at the very low temperature for , i.e., in the strongly correlated metallic regime. As HF-QMC results had so far only been reported for temperatures , these parameters have been believed to be out of reach of the Hirsch-Fye QMC method (7).

Discrete HF-QMC estimates of the imaginary-time Green function are shown as symbols in the main panel of Fig. 1

for relatively large values of the discretization . Clearly, the Trotter errors in these data sets are very significant at least for , much larger than the statistical and convergency errors (of about ). However, extrapolation schemes that have successfully been used for observables such as quasiparticle weight , energies, and double occupancies (8); (9) cannot directly be applied here, since the grid is different for each data set. In addition, none of the data sets really covers the initial rapid drop of : even for the finest discretization , is already halved at the first nontrivial data point. In fact, according to the assumption (7) that a good resolution of the highly curved initial region was necessary, HF-QMC results for such coarse grids should not yield useful information at all; instead one would have to resort to much finer grids of (7).

These considerations, however, overlook an important physical point: The behavior of at small is restricted by second-order weak-coupling perturbation theory (which is exact for the curvature of at in the half-filled case; see below and Fig. 3). Consequently, reasonable accuracy at small can be expected from a model Green function based on a (here particle-hole symmetric) two-pole approximation to the self-energy, ; since the weight is fixed by the leading large-frequency asymptotics (10), only the position of the poles (at ) remains as a free parameter. As indicated by the solid line in Fig. 1, this simple model (for ) has all the features that we expect from the true Green function, including the rapid initial drop; in fact, a visual inspection of the trends suggests that it better approximates the true Green function (for small ) than any of the QMC data sets. Still, as we will show in the following, the latter contain the necessary information for computing very precise estimates of the true Green function at all ; the model will only be needed for interpolating the raw QMC data.

Extrapolation procedure – As a starting point for the extrapolation, we need accurate QMC data with reliable error bars. While arithmetic averages are clearly appropriate for combining the Green functions obtained in parallel QMC runs for the same impurity model, errors can be minimized (most notably in insulating phases) by using geometric averages for combining estimates of from different DMFT iterations. In practice, we transform all data to before taking arithmetic averages; the appropriateness of this logarithmic scale will become apparent below (cf. Fig. 4).

In a second step, the averaged QMC data sets (for different ) need to be transformed to a common grid. For this purpose, optimized two-pole approximations (as specified above) are obtained for each data set. Each difference (defined only on the specific grid ) is interpolated by a natural cubic spline which is evaluated on a fixed, much finer grid; finally, is added to these intermediate results (11). As seen in Fig. 2,

these interpolated curves (dashed and dotted lines in the main panel) are smooth and vary systematically, with a discretization dependence which vanishes at small .

In the third step and final step, least squares fits are performed for at each value of (on the fine grid), taking quadratic and quartic contributions in into account (13). The result of this procedure is shown as thick solid line in Fig. 2. In the inset, the approximate low- asymptotics of the Green function,

have been subtracted. At this scale, the discretization error in the QMC data is clearly visible; in contrast, our extrapolated result (thick solid line) is hardly distinguishable from a continuous-time QMC estimate (7) (thin solid line). The fluctuations (on a much larger scale for HF-QMC than for CT-QMC) suggest that both results have similar precision. The competitiveness of HF-QMC (plus extrapolation) at this extremely low temperature, i.e., for rather coarse grids, may appear surprising. However, as shown in the inset of Fig. 2, the strongest absolute deviations of the raw HF-QMC data occur at (while the relative errors are largest at ), i.e., at a position independent of . This establishes that the rapid initial decay of does not set the scale for useful values of in HF-QMC and explains the good performance of our extrapolation procedure.

The uniform convergence of the (interpolated) HF-QMC Green functions, in turn, may be traced back to the fact that their curvatures are asymptotically exact for . Indeed, the second derivatives are visibly dependent in the main panel of Fig. 3

only for ; they all approach the correct asymptotic limit and agree at large within statistical errors. In contrast, as shown in the inset of Fig. 3, the Trotter errors of the first derivatives do not vanish at ; concurrently, the asymptotic exact value is nonuniversal, i.e., dependent on temperature and phase (metallic or insulating). Note that the extrapolated HF-QMC results in Fig. 3 (thick solid lines) are hardly distinguishable from the corresponding CT-QMC results (thin solid lines).

The specific form of the Green function extrapolation procedure detailed above is based on the insight that the Green function varies on a logarithmic scale. This is clearly seen in Fig. 4

for a moderately low temperature and similarly strong interaction (close to the thermodynamic Mott transition): While all curves become indistinguishable for , both differences between metal and insulator and the dependence in the insulating phase involve many orders of magnitude; in the latter case, even the error bars span nearly an order of magnitude. Obviously, a direct extrapolation of the insulating Green function (instead of ) would be hopeless for . Note that our choice guarantees the correct sign for the extrapolated Green function.

As seen in Fig. 5,

the absolute Trotter errors (i.e., the deviations at finite from the numerically exact results) of all data sets peak, again, at the same position (here , indicated by vertical lines). The remarkable constancy of this peak position, even across phase transitions, implies that the extrapolation technique can be trusted to capture the low- features of the Green function without bias, even for relatively coarse grids . In comparison, the region is more challenging (in the insulating phase): since HF-QMC results for (dotted line in Fig. 4) are too metallic, a good resolution of the exponential decay of (see lower solid line in Fig. 4) can only be expected from extrapolations which include small discretizations .

Spectra on the real axis – So far, we have only specified how to extract Green functions in the imaginary-time domain. The analytic continuation of this data to the real axis is still a highly nontrivial task; in fact, the extrapolated imaginary-time data arguably deserves even more careful and sophisticated continuation procedures than regular HF-QMC results, since – for the first time – the Green functions are exact within error bars. Adapted maximum entropy procedures and methods for efficient error analysis will be discussed elsewhere (14). In the case of very precise data (as in our examples), however, one may also neglect the small statistical errors and use Padé methods. Specifically, the local spectral function shown in the main panel of Fig. 6

has been obtained via a Padé approximant (15) to the imaginary-frequency self-energy using the first 199 positive Matsubara frequencies. This procedure is stable [as has been checked via sum rules for and and by varying the Padé parameters] and asymptotically exact in the noninteracting limit. In addition, it gives direct access also to momentum-resolved spectra (corresponding to ARPES measurements), as visualized in the left inset of Fig. 6.

The expected pinning of the quasiparticle peak height in to its noninteracting value (cf. dashed line in main panel) is accurately observed; the spectrum is also reasonably smooth without any indications of artefacts. clearly resolves the dispersions of the heavy quasiparticles (with ) and of the incoherent Hubbard bands. So the results appear reasonable – but are they significantly better that those attainable using conventional HF-QMC/MEM procedures? This is indeed the case, as demonstrated in the right hand inset of Fig. 6: here the thin solid line, computed from the final spectrum via Eq. 1, is hardly distinguishable from the numerically exact (thick line; cf. Fig. 2) while the deviations of finite- data (on which spectra would be based in conventional methods) are larger by several orders of magnitude. Thus, the spectrum of Fig. 6 is seen to be unbiased and numerically exact – for a temperature which had been deemed out of reach of HF-QMC.

Discussion – We have presented a method that, finally, allows to obtain numerically exact Green functions and ( resolved) spectral functions from regular Hirsch-Fye QMC calculations. Due to the simplifications in DMFT, transport properties such as the optical conductivity could be computed without additional effort. Since the method has been successfully tested also for doped systems (14), it should greatly improve the attainable quality of spectra, in particular in the context of LDA+DMFT calculations. In principle, the recently developed continuous-time QMC methods should yield spectra of similar quality with comparable effort; however, this might require regularization procedures which fail in insulating phases (16). Numerical renormalization group methods appear hampered by their coarse high-frequency resolution and systematic errors. It contrast, detailed comparisons of numerically exact HF-QMC spectra with ground state estimates from dynamical density-matrix renormalization group are expected to lead to new physical insight. Specificly, one may hope to thereby resolve the still controversial question of how the structure of the Hubbard bands changes across the Mott transition.

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 A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003) 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).
- M. Jarrell and J. E. Gubernatis, Phys. Rep. 269, 134 (1996).
- K. Byczuk, M. Kollar, K. Held, Y.-F. Yang, I. A. Nekrasov, Th. Pruschke, and D. Vollhardt, Nature Physics 3, 168 (2007).
- P. Werner, A. Comanac, L. de Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
- 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).
- M. Potthoff, T. Wegner, and W. Nolting, Phys. Rev. B 55, 16132 (1997).
- A similar procedure is used for Fourier transformations within our DMFT-QMC code (12); (8); (9).
- C. Knecht, N. Blümer, and P.G.J. van Dongen, Phys. Rev. B 72, 081103(R) (2005).
- For stabilization of the least-squares fits, lower- data is overweighted by rescaling errors by .
- N. Blümer, in preparation.
- H. J. Vidberg and J. W. Serene, JLTP 29, 179 (1977).
- I. S. Krivenko and N. A. Rubtsov, arXiv:cond-mat/ 0612233v1.