Quantum Monte Carlo solution of the dynamical mean field equations in real time

TL;DR

This paper demonstrates that real-time inchworm quantum Monte Carlo methods can effectively solve dynamical mean field equations, enabling accurate non-equilibrium simulations without analytic continuation issues.

Contribution

The paper introduces a numerically exact real-time Monte Carlo approach for dynamical mean field theory that overcomes the dynamical sign problem and is suitable for time-dependent non-equilibrium problems.

Findings

Successfully applied inchworm Monte Carlo to DMFT on Bethe lattice

Avoided analytic continuation issues in real-time calculations

Enabled simulation of far-from-equilibrium dynamics

Abstract

We present real-time inchworm quantum Monte Carlo results for single-site dynamical mean field theory on an infinite coordination number Bethe lattice. Our numerically exact results are obtained on the L-shaped Keldysh contour and, being evaluated in real-time, avoid the analytic continuation issues typically encountered in Monte Carlo calculations. Our results show that inchworm Monte Carlo methods have now reached a state where they can be used as dynamical mean field impurity solvers and the dynamical sign problem can be overcome. As non-equilibrium problems can be simulated at the same cost, we envisage the main use of these methods as dynamical mean field solvers for time-dependent problems far from equilibrium.

Equations

$$\hat{H}=-\tilde{v}\sum _{⟨ij⟩}{c}_{i\sigma }^{†}{c}_{j\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$G_{\sigma }(t,t^{\prime })\equiv -i⟨{\mathcal{T}}_{\mathcal{C}}{c}_{\sigma }(t)c_{\sigma }^{†}(t^{\prime })⟩,$$

$$G_{\sigma }^{\mathrm{ret}}(t,t^{\prime })=-i\theta (t-t^{\prime })⟨{\left[c_{\sigma }(t),c_{\sigma }^{†}(t^{\prime })\right]}_{+}⟩.$$

$${\Delta }_{\sigma }(t-t^{\prime })=G_{\sigma }(t-t^{\prime })$$

$$\begin{gathered} G^{\mathrm{ret}}(t-t^{\prime })=G_0^{\mathrm{ret}}(t-t^{\prime })+ \\ +{\iint }_0^{t_{\max }}d{t}_{1}d{t}_2{G}_0^{\mathrm{ret}}(t-t_1){\Sigma }^{\mathrm{ret}}(t_1-t_2)G^{\mathrm{ret}}(t_2-t^{\prime }). \end{gathered}$$

$$\begin{gathered} G^{\mathrm{ret}}(t_i-t_j)=G_0^{\mathrm{ret}}(t_i-t_j)+ \\ +(\Delta t{)}^2\sum _{k,l=0}^{N_t-1}{G}_0^{\mathrm{ret}}(t_i-t_k){\Sigma }^{\mathrm{ret}}(t_k-t_l)G^{\mathrm{ret}}(t_l-t_j)w_{kl}. \end{gathered}$$

$${\Sigma }_{ij}^{\mathrm{ret}}{w}_{ij}=\frac{(G_0^{\mathrm{ret}}{)}_{ij}^{-1}-(G^{\mathrm{ret}}{)}_{ij}^{-1}}{(\Delta t{)}^2}.$$

$$w_{i\geq j}=\left\{\begin{array}[]{ll}0,&i=0,\\ 0,&j=N_{t}-1,\\ 0,&i=j,\\ 1/2,&j=i-1,\\ 1,&\mathrm{otherwise}.\end{array}\right.$$

$$\tilde{A}(\omega )={\int }_{-\infty }^{\infty }d{\omega }^{\prime }K(\omega ,{\omega }^{\prime })A({\omega }^{\prime }),$$

$$K(\omega ,{\omega }^{\prime })=\frac{t_{\max }}{2\pi }\mathrm{sinc}\left(\frac{\left(\omega -{\omega }^{\prime }\right)t_{\max }}{2}\right).$$

$${\tilde{x}}_n=-\sum _{i=1}^p{a}_i{x}_{n-i}.$$

$$Q\mathbf{a}=-\mathbf{x},$$

$$R\mathbf{a}=-\mathbf{r}.$$

$$R_{ji}=\sum _n{x}_{n-j}^{\ast }{x}_{n-i};$$

$${\Sigma }^{\mathrm{ret}}(t)={\int }_{-\infty }^{+\infty }\frac{d\omega }{2\pi }{e}^{-i\omega t}{\Sigma }^{\mathrm{ret}}(\omega ).$$

$${\Sigma }^{\mathrm{ret}}(0^{+})=\underset{t\to 0^{+}}{\lim }\frac{1}{t}{\int }_{-\infty }^{+\infty }\frac{dz}{2\pi }{e}^{-iz}{\Sigma }^{\mathrm{ret}}(z/t).$$

$${\Sigma }^{\mathrm{ret}}(\omega )=\underset{\delta \to 0^{+}}{\lim }\sum _{m=1}^{\infty }\frac{C^{(m)}}{(\omega +i\delta {)}^m},$$

$${\Sigma }^{\mathrm{ret}}(0^{+})=\underset{t\to 0^{+}}{\lim }\frac{1}{t}\underset{\delta \to 0^{+}}{\lim }\sum _{m=1}^{\infty }{\int }_{-\infty }^{+\infty }\frac{dz}{2\pi }{e}^{-iz}\frac{t^m{C}^{(m)}}{(z+i\delta t{)}^m}$$

$$\begin{gathered} {\Sigma }^{\mathrm{ret}}(0^{+})=\underset{t\to 0^{+}}{\lim }\underset{\delta \to 0^{+}}{\lim }{\int }_{-\infty }^{+\infty }\frac{dz}{2\pi }{e}^{-iz}\frac{C^{(1)}}{z+i\delta t}= \\ =\underset{t\to 0^{+}}{\lim }\underset{\delta \to 0^{+}}{\lim }(-i{e}^{-\delta t}{C}^{(1)})=-i{C}^{(1)}. \end{gathered}$$

$${\Sigma }^{\mathrm{ret}}(0^{+})=-i{U}^2/4.$$

Full text

I Introduction

The direct solution of quantum systems with many interacting degrees of freedom is believed to be intractable in general, and in order to understand the salient aspects of these systems, suitable approximations have to be employed. The dynamical mean field theoryGeorges et al. (1996) is one such approximation. It is based on the realization that if the momentum dependence of the self-energy can be neglected, as occurs in certain infinite coordination number limits, the solution of a lattice model can be mapped onto that of an auxiliary quantum impurity model with self-consistently determined parameters.Georges and Kotliar (1992); Metzner and Vollhardt (1989)

Quantum impurity models are amenable to numerical study, and the last decade has seen rapid advancement in the development of algorithms for their solution. In particular, remarkable progress was achieved by the continuous-time quantum Monte Carlo methods,Rubtsov et al. (2005); Werner and Millis (2006); Gull et al. (2008); Otsuki et al. (2007); Prokof’ev and Svistunov (2007, 2008) which by now are the standard methods for studying multi-orbital and cluster impurity problems;Gull et al. (2011a) as well as in numerical renormalization group;Bulla et al. (2008) density matrix renormalization group;Wolf et al. (2014) and configuration interaction methods.Knowles and Handy (1984); Zgid et al. (2012)

Modern quantum Monte Carlo methods are numerically exact, in the sense that they can provide results that converge to the exact answer with an uncertainty proportional to the square root of the number of stochastic samples taken. However, because these methods are formulated within an imaginary-time statistical mechanics formalism, real-frequency data such as spectral functions needs to be extracted from an ill-conditioned analytic continuation procedure,Jarrell and Gubernatis (1996); Bryan (1990); Vidberg and Serene (1977); Sandvik (1998); Beach (2004); Mishchenko et al. (2000); Fuchs et al. (2010); Prokof’ev and Svistunov (2013); Otsuki et al. (2017) in which these uncertainties are exponentially amplified. Furthermore, studying systems in general nonequilibrium states or under time-dependent driving beyond linear response is not possible.

The desire to understand nonequilibrium transport in correlated impurity models motivated the development of real-time generalizations of continuous-time quantum Monte Carlo methods.Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009); Werner et al. (2010); Antipov et al. (2016) These early methods require exponentially increasing computer time as a function of the simulated time due to the dynamical sign problem, i.e. the sign problem occurring when real-time dynamics is evaluated. The development of bold-line impurity solver algorithmsGull et al. (2010, 2011b); Cohen et al. (2014a, b) substantially alleviated this sign problem and increased the accessible parameter space. A recent development, the inchworm quantum Monte Carlo method,Cohen et al. (2015) showed a reduction of the computational scaling from exponential to quadratic, effectively eliminating the dynamical sign problem altogether.

For dynamical mean field applicationsEckstein et al. (2009, 2010); Eckstein and Werner (2011); Tsuji et al. (2011) of the inchworm method, two more components are necessary: the ability to obtain two-time response functions, and the extension of the method to an initial thermal ensemble. Both components have recently been implemented,Antipov et al. (2017) allowing us to test the method. As a first application, we address a lattice model for which the dynamical mean field approximation becomes exact: the Hubbard model on the infinite coordination number Bethe lattice. While we see the main application of our method in non-equilibrium,Aoki et al. (2014) we demonstrate its power here for the equilibrium case, where a large variety of well-developed and competitive methods are available and the physics is well understood.

The paper is organized as follows. In Section II, we write down the lattice model, recapitulate the dynamical mean field theory, introduce our impurity solver, discuss how self-energies are extracted, and explain how linear prediction can be used to obtain spectral functions. In Section III we present our results, including real-time Green’s functions, real-time self-energies and real-frequency spectral functions with and without linear prediction. Finally, in Section IV, we discuss our conclusions and outlook. The appendices contain further technical and numerical information.

II Methods

II.1 Lattice model and dynamical mean field theory

We study the repulsive Hubbard model in equilibrium on the infinite coordination number Bethe lattice. The Hamiltonian is given by

[TABLE]

where $\tilde{v}$ denotes the hopping matrix element, and $\sigma$ denotes the spin index. The operators $c^{\dagger}_{i\sigma}$/$c_{i\sigma}$ create/annihilate electrons with the spin $\sigma$ on the corresponding site $i$. $U$ is the on-site Coulomb repulsion between electrons with opposite spins. We restrict our discussion to the paramagnetic solution of a half-filled infinite coordination number Bethe lattice. The hopping matrix element must be properly scaled with the coordination number $Z$ to remain finite in the limiting case, $\lim_{Z\to\infty}(\tilde{v}\sqrt{Z})=v$.Metzner and Vollhardt (1989) Throughout this paper we employ dimensionless units by dividing/multiplying all energy/time scales by the hopping strength $v$. The Bethe lattice is characterized by a semi-elliptical noninteracting density of states $D(\epsilon)=\sqrt{4-\epsilon^{2}}/(2\pi)$. The model’s interacting self-energy $\Sigma_{i\sigma,j\sigma^{\prime}}(t,t^{\prime})$, corresponding to Eq. 1, is purely local (i.e. zero for $i\neq j$), and its local part $\Sigma_{i\sigma,i\sigma^{\prime}}(t,t^{\prime})$ is equal to the self-energy of an auxiliary impurity model. The lattice Green’s function can easily be obtained from this self-energy.Georges et al. (1996)

We simulate the real-time dynamics of the model in equilibrium by means of the real-time dynamical mean field formulated on the $L$-shaped Keldysh contour.Schmidt and Monien (2002); Freericks et al. (2006) We define the impurity Green’s function as the contour-ordered expectation value Aoki et al. (2014)

[TABLE]

where $c_{\sigma}$ and $c_{\sigma}^{\dagger}$ denote impurity operators, such that the retarded Green’s function is given by

[TABLE]

Time-translation invariance of the system in equilibrium implies that the two-time Green’s function only depends on time differences: $G_{\sigma}(t,t^{\prime})=G_{\sigma}(t-t^{\prime})$. The DMFT self-consistency condition for the Bethe lattice thus reads

[TABLE]

for any two times $t,t^{\prime}$ on the L-shaped contour.

We solve the dynamical mean field equations iteratively, starting from a metallic, insulating, or high-temperature guess for the Green’s function.

II.2 Inchworm impurity solver

Each DMFT iteration requires the solution of an impurity problem. We employ the inchworm quantum Monte Carlo solverCohen et al. (2015); Antipov et al. (2017) to obtain a numerically exact solution of the Green’s function of this single impurity Anderson model, up to some finite maximal time.

The inchworm method used here is a continuous-time (or diagrammatic) Monte Carlo methodGull et al. (2011a) based on sampling Feynman diagrams from a hybridization expansion. Unlike standard continuous-time methods, the inchworm method takes advantage of the causal structure of propagation: information from easy-to-evaluate propagation over short times is reused to more efficiently evaluate propagation over longer times. The method has been shown to circumvent the dynamical sign problem in several models, diagrammatic expansions, and physical regimes.Cohen et al. (2015); Antipov et al. (2017); Chen et al. (2017a, b) Its total run time (when applied to the hybridization expansion) grows quadratically with the simulated time. In practice, for the model and initial condition treated in the present work, times about an order of magnitude larger than possible with a bare algorithm can be accessed. Details of the solver implementation, and in particular of the Keldysh Green’s function measurement, are reported in Ref. Antipov et al., 2017.

II.3 Real-time self-energies

The local self energy function is of special interest in the DMFT context as it contains all relevant information about single-particle correlations in the system. Its retarded component $\Sigma^{\mathrm{ret}}(t-t^{\prime})$ is defined as the solution of the Dyson equation (spin indices are omitted for clarity)

[TABLE]

Here, $G_{0}^{\mathrm{ret}}(t-t^{\prime})$ is the retarded Green’s function of the noninteracting lattice, $G^{\mathrm{ret}}(t-t^{\prime})$ is obtained numerically as a result of a DMFT/inchworm calculation, and $t_{\mathrm{max}}$ is a maximum simulation time.

Eq. (5) is a Volterra integral equation of the first kind with respect to $\Sigma^{\mathrm{ret}}$. In a numerical implementation, $G^{\mathrm{ret}}_{0}(t)$ and $G^{\mathrm{ret}}(t)$ are known on some finite time mesh (for example, a uniform grid). It is therefore natural to project the integral equation onto the mesh, and use numerical linear algebra methods to solve the resulting linear system. In principle, the solution of Eq. (5) could be obtained in Fourier space by applying the convolution theorem. However, this approach is not practical here, as $G^{\mathrm{ret}}(t)$ is known only up to a finite time, and does not generalize to the non-equilibrium case without additional modifications.

Given a number of time slices $N_{t}$, we introduce mesh nodes $t_{i}=i\Delta t$ with $i\in\left\{0,1,\ldots,N_{t}-1\right\}$ and $\Delta t=t_{\mathrm{max}}/(N_{t}-1)$. A discretized version of the Dyson equation then reads

[TABLE]

The quadrature weights $w_{kl}$ define the integration method, and at this point we do not specify their precise form. Using matrix notation $F_{ij}=F(t_{i}-t_{j})$ we find

[TABLE]

The numerical matrix inversion used here is stable, as all retarded Green’s functions are represented by lower triangular matrices with $G^{\mathrm{ret}}(0^{+})=-i$ on the main diagonal, such that their condition number is 1. The $(\Delta t)^{2}$ in the denominator of (7) suggests that this procedure is similar to numerical calculation of the second derivative. The number of time slices at which $G^{\mathrm{ret}}(t_{i})$ is known is limited by the computational effort required by the inchworm QMC algorithm. In order to make the numerical differentiation accurate, we choose $N_{t}$ much larger than used in the inchworm simulation and employ a cubic spline interpolation of $G^{\mathrm{ret}}$ between the original nodes.

Only certain choices of quadrature weights give a convergent solution in the small $\Delta t$ limit. In his study of the one-dimensional Volterra equations of the first kind, Linz (1969) showed that trapezoidal and rectangular rules are convergent, whereas higher order quadrature methods in general are not. Using this result as a starting point, we construct $w_{ij}$ as possible direct products $w_{i}w_{j}$, where $w_{i}$ and $w_{j}$ correspond to the rectangular rule with the first/last point excluded and to the trapezoidal rule. Most combinations can be ruled out, as they cause $\Sigma^{\mathrm{ret}}(t_{\mathrm{max}})$ to diverge in the small $\Delta t$ limit. We found that the following choice:

[TABLE]

yields stable and accurate results. The first two lines in (8) show that we choose a rectangle rule approximation excluding the first slice from the $t_{1}$-integral and the last slice from the $t_{2}$-integral. Values on the first sub-diagonal receive half the weight, because retarded functions are proportional to $\theta(t-t^{\prime})$. Finally, $w$ as well as all other matrices entering the equation must be of Toeplitz form in equilibrium. We therefore also set the main diagonal of $w$ to zero.

Equations (7) and (8) allow for extracting all mesh values of the self energy, except $\Sigma^{\mathrm{ret}}_{ii}=\Sigma^{\mathrm{ret}}(0^{+})$. This element is known analytically, see Appendix D.

II.4 Linear Prediction

Our simulations are performed up to a finite maximum real time $t_{\rm max}$ determined by the available computational resources. However, many quantities are best described in the real-frequency domain, and are therefore expressed as Fourier transforms over the entire time axis. When Fourier transforming quantities with a hard time cutoff, the transform can be expressed as a convolution with a sampling kernel proportional to $\text{sinc}(\omega)=\frac{\sin(\omega)}{\omega}$ in the frequency domain, i.e.

[TABLE]

where $\tilde{A}(\omega)$ is the result from a Fourier transform with data up to finite times, $A(\omega)$ is the true result and the convolution kernel is given by

[TABLE]

This convolution introduces broadening and unphysical oscillations into the spectral function. Linear prediction is a technique to remove these artifacts by using a physically motivated extrapolation scheme to extend data beyond the maximum simulated time. Linear prediction has previously been used for this purpose within the framework of $t$-DMRGWhite and Affleck (2008); Barthel et al. (2009); Pereira et al. (2012); Ren and Sirker (2012). We start with the ansatz that the value of the signal (in this case the Green’s function) at the $n$-th real time point is a linear function of its previous $p$ values, i.e.

[TABLE]

It can be shown that this corresponds to fitting the function in time to a superposition of $p$ complex exponential terms. In the case of a Green’s function dominated by a few isolated poles, this approximation is justified. For the data presented here, the validity needs to be assessed by systematically varying $p$.

In order to use this model for extrapolation, the coefficients $a_{i}$ must be fit to the known data. This is done over the region $\left(t_{\rm max}-t_{\rm fit},t_{\rm max}\right)$ in order to exclude spurious short time behavior from the fit results. The linear prediction ansatz leads to a matrix equation

[TABLE]

where $Q_{nk}=x_{n-k}$ is an $N\times p$ complex matrix with $N$ the number of points in $\left(t_{\rm max}-t_{\rm fit},t_{\rm max}\right)$. Solving this in the least squares sense leads to the normal equations

[TABLE]

This is written in terms of the autocorrelations of the data,

[TABLE]

These matrix equations are often unstable and require some form of regularization. Here we choose the simple regularization, $R^{-1}\to\left(R+\epsilon I\right)^{-1}$ and check that the results are not strongly affected by the regularization parameter.Barthel et al. (2009) With this scheme, the coefficients $a_{i}$ are readily obtained and the Green’s function can be extrapolated until it decays to zero. The procedure yields spectral functions free from unphysical finite time oscillations, at the cost of some additional systematic uncertainty due to the assumptions imposed by the linear prediction ansatz. We emphasize that linear prediction is only applied to our converged data as a post-processing routine, in order to interpret the real-time results as functions of frequency. The DMFT / inchworm iteration procedure preceding it is formulated in terms of real, finite times only, and independent of the linear prediction formalism.

III Results

Our impurity solver obtains results for $G(t,t^{\prime})$ as a function of real times. A typical raw simulation output is shown in Fig. 1, which shows the real (solid lines) and the imaginary part (dashed lines) of the greater Green’s function up to a maximum time $tv=2$ and at a low temperature of $T/v=0.02$. The interaction strength is taken to be $U/v=4$, equal to the full bandwidth of the infinite coordination number Bethe lattice. The error bars of the measured Green’s function could, in principle, be estimated as standard deviation from a set of completely independent DMFT/inchworm runs. With our present implementation this approach has proven too computationally expensive.

The inchworm Monte Carlo method is only exact when two numerical parameters are controlled. The first of the parameters is the discretization of the imaginary and the real time branch $\Delta t$, which we chose to be $\Delta tv=0.05$. Second, the maximum order to which diagrams in the inchworm expansion are considered. We find, especially in the metallic low temperature regime, that results converge at an inchworm expansion order of around seven. As the inchworm order is directly related to the number of crossings considered in an $N$-crossing approximation,Cohen et al. (2015) this result implies that results from non-crossing or one-crossing approximations are not valid in this parameter regime. Throughout this paper, all the results presented are converged in both $\Delta t$ and maximum order.

We observe that the dynamical mean field solution converges in a causal manner, in the sense that results are converged within one iteration by time $tv=0.3$, within two iterations by time $tv=1.0$, and results up to $tv=2$ are indistinguishable between iterations $4$ and $5$, indicating that the self-consistency loop has converged. While this causal convergence can be used to avoid the usual dynamical mean field iteration scheme and replace it with a direct solution of the self-consistency equations,Aoki et al. (2014) this scheme has not been employed here.

As the model is tuned from metallic to insulating, the shape of the real-time Green’s functions changes substantially, from slow oscillations at weak coupling to rather rapid, quickly decaying oscillations at large interaction strength. This is shown in Fig. 2, for interaction strengths $U=2/v$, $U/v=4$, $U/v=6$, and $U/v=8$, at a relatively high temperature of $T/v=0.5.$

In contrast, lowering the temperature by a factor of $10$, as shown in Fig. 3, causes relatively little change in the Green’s function, with both oscillation frequency and amplitude staying more-or-less invariant for the time range simulated. Fig. 3 shows a regime in which a quasi-particle peak emerges upon cooling. However, the properties of this peak are not obvious in the data for the greater Green’s function shown here.

Using the procedure described in Sec. II.3, we can directly extract a real-time self-energy. Fig. 4 shows the imaginary parts of the noninteracting and the interacting retarded Green’s functions (left vertical axis) along with the imaginary parts of the computed self-energy (right vertical axis) as a function of time. The orange curve is the Green’s function reconstructed from replacing the bare Green’s function and the self-energy into eq. (6), showing that the scheme for extracting the real-time self-energy, which is numerically demanding, has converged. The numerical noise visible in the real-time self-energy may be used to qualitatively assess the size of the Monte Carlo errors intrinsic to this simulation. However, we note that in the case of inchworm simulations this approach typically underestimates the errors, and a more complicated procedure is required if one is interested in rigorous error estimates.Cohen et al. (2015); Antipov et al. (2017)

Fig. 5 shows the evolution of the imaginary part of the real-time self-energy with interaction strength, for the same parameters as shown in Fig. 2. The self-energy decays to zero within the accessible times in the weak coupling limit.

In contrast, the insulating regime shows a self-energy consistent with a constant in the long-time limit. This behavior is caused by the DMFT mechanism for opening a gap, which requires a pole at zero frequency (at half-filling) in the self-energy, such that $\text{Im}\Sigma(\omega)=\Delta^{2}\delta(\omega)+\text{Im}\Sigma^{\text{reg}}(\omega)$, where $\Delta$ denotes the Mott half-gap size and $\Sigma^{\text{reg}}$ the non-divergent part of the self-energy. In the atomic limit, one would expect $\Delta\sim U/2$. For the lattice problem with non-zero hybridization, the gap sizes are smaller.

Interacting spectral functions in real frequency, $A(\omega)=-\frac{1}{\pi}\text{Im}G^{\mathrm{ret}}(\omega)$, are of principal interest in dynamical mean field calculations because they allow for direct comparison with photoemission experiments. Obtaining $A(\omega)$ in imaginary-time formulations requires an ill-conditioned analytical continuation procedure, such as the maximum entropy algorithm. The real-time formulation avoids this. However, the fact that the real-time Green’s function is only known up to a finite maximum time implies that spectral functions can only be resolved with a frequency resolution proportional to the inverse of that maximum time. Fig. 6 shows spectral functions extracted from the data in Fig. 2 and corresponding to the self-energies in Fig. 5. It demonstrates the metal-insulator crossover as the on-site interaction strength increases. Due to the finite frequency resolution, sharp features are absent and part of the spectral function turns negative, especially in the higher frequency regions.

The linear prediction method described in Sec. II.4 is an interpolation routine designed to replace the sharp cutoff of $G(t)$ at the maximum time $t_{\text{max}}$ with a smoothly decaying function corresponding to a set of poles in the complex plane. As is evident in Fig. 7, data obtained for times up to $tv=2$ and extrapolated up to $tv=4$ using linear prediction approximates well our data directly computed by running the dynamical mean field algorithm up to $tv=4$.

The linear prediction results can then be used to obtain spectral functions that do not suffer from the ‘ringing’ phenomenon. Fig. 8 shows the results of this procedure in practice: while the straightforward continuation of the data up to time $tv=2$ shows only a coarse frequency resolution and has substantial negative contributions between frequencies of $\omega/v=5$ and $10$, the corresponding extrapolated data shows a double-peak feature and is positive. Longer-time data (not based on linear prediction) corroborates the double-peak structure and the slightly larger gap, and similarly does not exhibit a negative contribution to the spectral function.

Linear prediction is an extrapolation procedure that needs to be carefully controlled. Results should not depend on the choice of cutoff time or the number of poles interpolated. These concerns are addressed in Appendix B.

Using linear prediction, we revisit two aspects of single-site dynamical mean field theory: the opening of the Mott gap as interaction strength is increased, and the establishment of a quasiparticle peak as temperature is lowered in the coherent metallic regime.

Fig. 9 shows the data of Fig. 6 obtained within linear prediction. It is evident that the increased frequency resolution leads to additional features in the spectral function. $U/v=2$ is metallic with little change of the band structure due to correlations. $U/v=4$ shows ‘bad metallic’ behavior where the spectral function near zero is suppressed due to the onset of insulating correlations. As the interaction is raised to $U/v=6$, a gap opens and a double-peak feature develops, and by $U/v=8$ a clear insulating structure with a wide Hubbard gap has developed.

As temperature is lowered at $U/v=4$, the ‘bad metal’ dip in the spectral function disappears, and for a temperature of $T/v=0.1$, a clear sign of a ‘quasiparticle peak’ has developed (Fig. 10). Further lowering of the temperature makes the quasiparticle peak more pronounced, while deepening the ‘dip’ between the Hubbard side bands and the peak.

IV Conclusions

We have presented a first application of real-time quantum Monte Carlo methods to real-time dynamical mean field theory. We showed that it is possible to obtain accurate Green’s functions and self-energies directly in real-time. We further showed that if one is interested in frequency-transformed properties, linear prediction methods for extending the range of the available data to longer times work in practice for Monte Carlo data with sufficient accuracy.

We expect that in the near future our algorithms will have applications mainly out of equilibrium, in the study of quenches and driven problems, which the impurity solver algorithm demonstrated here can simulate at no additional cost. In contrast, reaching substantially longer times or substantially lower temperatures than the ones demonstrated here will require additional optimizations and access to supercomputer resources.

To take a long view on the equilibrium problem, it is useful to consider the respective scaling properties of imaginary-time and real-time algorithms with regard to the desired accuracy of real frequency quantities. Whereas imaginary-time algorithms are currently far more efficient, they are limited by the exponential sensitivity to errors in the various analytical continuation methods. Real-time algorithms are substantially more expensive, but their accuracy is limited only by the maximum simulated time. In the inchworm method, extending this time requires a quadratic increase in computational effort. It is therefore entirely conceivable that in the future—as the field grows to be interested in larger, more frustrated problems and in higher frequency resolution—the real-time algorithms will eventually surpass their imaginary-time counterparts in efficiency by the sheer power of scaling.

V Acknowledgments

Funding for this project was provided by DOE ER 46932. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. GC was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 1604/16)

Appendix A Convergence with respect to the inchworm

truncation order

Fig. 11 shows the comparison of the Matsubara Green’s function $G(\tau)$ computed by the inchworm nonequilibrium QMC solver and an equilibrium continuous-time hybridization expansion solverWerner et al. (2006) (we used the solver of the open source TRIQS libraryParcollet et al. (2015); Seth et al. (2016)). The data are obtained for an impurity calculation starting with a semielliptic density of states with a full bandwidth of $4$ at $U/v=4.6$ and half-filling. The upper panel of Fig. 11 demonstrates that at high temperature $T/v=0.2$, a maximum order of $5$ in the inchworm order truncation is sufficient. In contrast, as the temperature is lowered to $T/v=0.04$ (near the onset of the first order coexistence in the dynamical mean field solution of this model), convergence of the inchworm calculations to the equilibrium result is only achieved at orders $5,6$ and $7$. The correspondence between inchworm order truncation and the number of crossings in an N-crossing approximation (such as NCA, OCA, or the two-crossing approximation) implies that a large number of crossings is essential for obtaining correct results in the correlated metallic regime.

Appendix B Robustness of linear prediction

Linear prediction, introduced in Sec. II.4 in the main text, is a method used to extrapolate data known up to a finite cutoff time to much longer times, such that smooth spectral functions can be extracted. The method can only succeed if the underlying data contains enough information to represent the long-time behavior accurately, and is expected to fail if this is not the case. In practice, two control parameters are available. First, the number of components (poles) $p$ that is being fitted. Second, the time interval over which the data is extrapolated.

Fig. 7 and Fig. 8 in the main text show the behavior of the spectral function and the real-time Green’s function on the final time. Fig. 12 and Fig. 13 show the dependence of the converged DMFT spectral functions on the number of poles $p$ and the maximum fitting time. Plotted are the converged dynamical mean field spectral functions obtained by extrapolating the real-time Green’s function from $tv=2$ to $tv=10$.

Clear differences in the data are visible in the height of the quasi-particle peak, the height of the four peaks in the insulating regime, and the size of the dip separating the quasiparticle feature from the Hubbard side bands. The ambiguities become less pronounced when the time accessed is increased. In any case, the main features, in particular the existence of a quasiparticle peak in the metal or of a double peak structure above and below the Mott gap are clearly evident in our data.

Appendix C Self-energy extraction: technical issues

Being limited by the complexity of the inchworm algorithm, we have to take a sparse time grid for the real-time Green’s function $G_{\sigma}(t)$ and its retarded counterpart $G^{\mathrm{ret}}_{\sigma}(t)$. The grid step $\Delta tv$ used throughout this work is equal to 0.05, while $t_{\mathrm{max}}v\in\{2.0,4.0\}$.

The self-energy extraction procedure presented in section II.3 involves a rectangular rule discretization of the time integrals, a first order scheme in $\Delta t$. We have found that effect of the discretization error rapidly grows as we propagate to larger times in $\Sigma^{\mathrm{ret}}_{\sigma}(t)$. It is therefore crucial to introduce a finer time grid, and interpolate $G^{\mathrm{ret}}_{\sigma}(t)$/ $G^{\mathrm{ret}}_{0,\sigma}(t)$ between the original nodes before doing matrix inversions. Fig. 14 shows self-energy extraction results for $T/v=0.5,\ U/v=8.0$, computed with $N_{t}=41$ (the original number of points), $101$, $1001$ and $2001$ interpolation slices. One can clearly see a drastic difference in the self-energy curves for $N_{t}=41$ and $2001$. On the other hand, there is no visible difference between $N_{t}=1001$ and $2001$, which means the interpolation has converged.

The dashed orange curves are obtained by back-substitution of $\Sigma^{\mathrm{ret}}(t)$ into a trapezoidal-rule discretization of the Dyson equation (5) (different choice of $w_{ij}$). For $N_{t}=1001,2001$ they lie on top of the input Green’s function, which also signals convergence.

Appendix D Short time limit of $\Sigma^{\mathrm{ret}}(t)$

[TABLE]

Let us introduce a rescaled frequency $z=\omega t$ and consider a short time limit of the self-energy,

[TABLE]

$\Sigma^{\mathrm{ret}}(z/t)$ is an analytic function in the upper half-plane of $z$ for any positive $t$. We now employ a high-frequency expansion of the self-energy in the absence of a Hartree-Fock term,

[TABLE]

[TABLE]

It is easy to see that only the $m=1$ term contributes in the limit of $t\to 0^{+}$.

[TABLE]

For the symmetric single orbital Anderson model $C^{(1)}=U^{2}/4$ (independent of bath parameters). For a derivation see, for instance, Ref. Potthoff et al., 1997. Given this,

[TABLE]