Quasi continuous-time impurity solver for the dynamical mean-field theory with linear scaling in the inverse temperature
aa r X i v : . [ c ond - m a t . s t r- e l ] M a y Quasi-continuous-time impurity solver for the dynamical mean-field theorywith linear scaling in the inverse temperature
D. Rost,
1, 2
F. Assaad, and N. Bl¨umer Institute of Physics, Johannes Gutenberg University, Mainz, Germany Graduate School Materials Science in Mainz, Johannes Gutenberg University, Mainz, Germany Institute of Theoretical Physics and Astrophysics, University of W¨urzburg, W¨urzburg, Germany (Dated: August 13, 2018)We present an algorithm for solving the self-consistency equations of the dynamical mean-fieldtheory (DMFT) with high precision and efficiency at low temperatures. In each DMFT iteration,the impurity problem is mapped to an auxiliary Hamiltonian, for which the Green function iscomputed by combining determinantal quantum Monte Carlo (BSS-QMC) calculations with amultigrid extrapolation procedure. The method is numerically exact, i.e., yields results which arefree of significant Trotter errors, but retains the BSS advantage, compared to direct QMC impuritysolvers, of linear (instead of cubic) scaling with the inverse temperature. The new algorithm isapplied to the half-filled Hubbard model close to the Mott transition; detailed comparisons withexact diagonalization, Hirsch-Fye QMC, and continuous-time QMC are provided.
PACS numbers: 02.70.Ss, 71.10.Fd, 71.27.+a
I. INTRODUCTION
The dynamical mean-field theory (DMFT)[1–4] andits cluster extensions [5, 6] are powerful approaches forthe numerical treatment of correlated electron systems,both in the model context and for materials science, e.g.,embedded in the LDA+DMFT [7] or GW+DMFT [8–11] frameworks which extend density functional theoryto strongly correlated materials [3, 12]. Recently, manyDMFT studies have also appeared in the context ofultracold fermions on optical lattices [13–15]. TheDMFT reduces electronic lattice models to impurityproblems, which have to be solved self-consistently [16–18]. A challenging part of this iterative procedure isthe computation of the interacting Green function fora given impurity configuration (defined by the fixed localinteractions and the self-consistent Weiss field). Thus,the availability of efficient and reliable impurity solvers determines the complexity of models and the parameterspace that can be accessed using the DMFT.Quantum Monte Carlo (QMC) impurity solvers al-low for numerically exact solutions of the DMFT self-consistency equations at finite temperatures. In thecase of the Hirsch-Fye auxiliary field (HF-QMC) method[16, 19, 20], all raw estimates contain systematic errorsdue to the inherent Trotter decomposition and associatedimaginary-time discretization [19, 21]; unbiased resultscan only be obtained after an extrapolation of the dis-cretization interval ∆ τ → β = 1 /k B T ,which limits their access to low-temperature phases.Exact diagonalization (ED) based impurity solvers [29]require a discrete representation of the impurity action interms of an auxiliary Hamiltonian, which is then solved either by full diagonalization (for evaluations at arbitrarytemperature) or using a Lanczos procedure [30] (e.g., at T = 0). As the numerical effort scales exponentially withthe number N b of auxiliary “bath” sites, N b has to bekept quite small, which introduces, again, a bias and is aparticularly severe limitation for multi-orbital or clusterDMFT studies at finite temperatures.Recently, Khatami et al. proposed anotherHamiltonian-based scheme [31], in which the Greenfunction and other relevant properties of the auxiliaryproblem are computed using the determinantal BSS-QMC method developed by Blankenbecler, Scalapino,and Sugar [32]. The advantage of this scheme, comparedto ED, is the possibility of using more bath sites (dueto cubic instead of exponential scaling with N b ); theadvantage over the direct QMC impurity solvers is thelinear, instead of cubic, scaling in β [33]. The authorsestablished the feasibility of the method and provedthat the associated sign problem (arising at generalband filling in cluster extensions of DMFT, in frustratedlattices, and for generic multiband models) converges tothat of HF-QMC for sufficiently fine bath discretization[31]. However, as all BSS-QMC applications to date, theGreen functions and all observable estimates resultingfrom their implementation suffer from systematic Trottererrors.In this work, we construct a similar algorithm wherethe Trotter bias inherent in the BSS Green functions iseliminated using a multigrid procedure before feedingthem back in the self-consistency cycle. As a DMFTbuilding block, the resulting method is an exact quasi-CT QMC impurity solver with linear scaling in theinverse temperature. Its scaling advantage over directQMC impurity solvers should allow access to lowertemperatures, in particular in multi-orbital and clusterDMFT studies.The paper is organized as follows: In Sec. II, webriefly review the DMFT equations and the BSS-QMCalgorithm and fully specify our multigrid DMFT-BSSapproach. As a test case, the new method is applied(in single-site DMFT) to the half-filled Hubbard modelin the vicinity of the Mott transition in Sec. III. Here,we first focus on the Green function at moderatelylow temperature T = t/
25 and then discuss importantobservables, namely the double occupancy and quasipar-ticle weight, also at lower temperatures. The accuracyof our approach is established by comparisons with theresults of the (multigrid) HF-QMC, ED and CT-QMCsolvers, as well as with the previous (finite ∆ τ ) DMFT-BSS implementation. We show that our elimination ofthe Trotter error improves the results dramatically. Wealso discuss the impact of the bath discretization andestablish convergence to the thermodynamic limit. Asummary and outlook conclude the paper in Sec. IV. II. THEORY AND ALGORITHMS
In this section, we lay out the proposed algorithmfor solving the DMFT self-consistency equations withoutsignificant Trotter errors and with a computational effortthat grows only linearly with the inverse temperature.We start out by reviewing the general DMFT frameworkand established methods (ED, HF-QMC, CT-QMC) forits solution in Sec. II A in sufficient detail to exposethe similarities and differences with respect to the newmethod. Here we also discuss some algorithmic choices,in particular regarding the Hamiltonian representationin our ED implementation, that are essential ingredientsalso for the BSS-QMC based approaches. We then turnto the BSS-QMC method and its applicability in theDMFT context in Sec. II B, and specify, finally, ournew numerically exact implementation in Sec. II C. Forsimplicity, and in line with the numerical results to bepresented in Sec. III, we write down the formalism forthe single-band Hubbard model and the original, single-site variant of DMFT. Extensions to cluster DMFT (andto multiband models) should be straightforward, butrequire some generalizations (e.g. for the treatment ofoffdiagonal Green functions) and will be pursued in asubsequent publication.
A. DMFT and established impurity solvers
The Hubbard model on a lattice or graph – We considerthe single-band Hubbard model H = H + H int = X ij,σ t ij c † iσ c jσ + U X i n i ↑ n i ↓ , (1)where c † iσ ( c iσ ) creates (annihilates) an electron withspin σ ∈ {↑ , ↓} on lattice site i ; n iσ = c † iσ c iσ is thecorresponding density, t ij = t ji the hopping amplitudebetween sites i and j (or the local potential for i = j ); U quantifies the on-site interaction. Usually, the hopping tU DMFT G U Aux. V V V V V V ǫ ǫ ǫ ǫ ǫ ǫ U (a) (b) (c) FIG. 1. (Color online) Mapping of the original lattice problem(a), with local interaction U and hopping t , on a single impu-rity (b), embedded in an effective bath G . (c) Discretizationof the bath in terms of an auxiliary Hamiltonian (treatablewith ED or BSS-QMC), here with star topology. is defined to be translationally invariant, e.g., t ij = − t for nearest-neighbor bonds on an infinite mathematicallattice [as illustrated in Fig. 1(a) for a square lattice]or on a finite cluster with periodic boundary conditions.However, neither the DMFT nor direct QMC approachesto the Hubbard model depend crucially on such assump-tions, as will be discussed in Sec. II B. General DMFT self-consistency procedure – If alllattice sites are equivalent and for spatially homogeneousphases, the DMFT maps the original lattice problem (1),illustrated in Fig. 1(a), onto a single-impurity Andersonmodel [Fig. 1(b)], which has to be solved self-consistently.The impurity problem is defined by its action A [ ψ, ψ ∗ , G ] = β Z β Z dτ dτ ′ X σ ψ ∗ σ ( τ ) G − σ ψ σ ( τ ′ ) (2) − U β Z dτ ψ ∗↑ ( τ ) ψ ↑ ( τ ) ψ ∗↓ ( τ ) ψ ↓ ( τ ) ,here in imaginary time τ ∈ [0 , β ] and in terms ofGrassmann fields ψ , ψ ∗ . G is the “bath” Green function,i.e., the noninteracting Green function of the impurity,which is related to the full impurity Green function G , G σ ( τ ) = −hT τ ψ σ ( τ ) ψ ∗ σ (0) i A (3)(with the time ordering operator T τ ), and the self-energyΣ by the (impurity) Dyson equation G − σ ( iω n ) = G − σ ( iω n ) − Σ σ ( iω n ) , (4)here written in terms of fermionic Matsubara frequencies ω n = (2 n + 1) πT at finite temperature T ; here and in thefollowing, we set ~ = k B = 1.The central DMFT assumption is that of a local self-energy on the lattice [1]: Σ ijσ ( iω n ) = δ ij Σ iiσ ( iω n ), whichis identified with the impurity self-energy. Similarly,the impurity Green function is identified with the localcomponent of the lattice Green function: G σ ( iω n ) = G iiσ ( iω n ) = { t + [ iω n + µ − Σ σ ( iω n )] } − ii = Z ∞−∞ dε ρ ( ε ) iω n + µ − Σ σ ( iω n ) − ε , (5) G Σ lattice Dyson eq. (5)impurity problem (3), (2) G = G [ G ; U, T ] G G impurity Dyson eq. (4) τ → iωn HF-QMC iωn → τ ∆ τ ED G →
Vi, εiNb τ → iωn ∆ τ → ... BSS-QMC ∆ τj G →
Vi, εiNb (a) (b) (c) (d)
FIG. 2. (Color online) (a) Scheme of the general DMFT self-consistency cycle, including the “impurity problem” (dashed box).Established impurity solvers include (b) the Hirsch-Fye (HF-QMC) algorithm and (c) exact diagonalization (ED): cf. the maintext. (d) The proposed algorithm approximates the bath Green function G in terms of the parameters V i , ǫ i of an auxiliaryHamiltonian (6) with N b “bath” sites [like ED (c)]. Corresponding Green functions are computed using BSS-QMC for a grid∆ τ min ≤ ∆ τ j ≤ ∆ τ max of Trotter discretizations. The subsequent extrapolation of ∆ τ → τ , which is easily Fourier transformed and fed back into the self-consistency cycle. where the last expression is valid in the homogeneouscase, ρ ( ε ) denotes the noninteracting density of states,and t is the matrix with elements t ij .The general DMFT iteration scheme is illustrated inFig. 2(a): starting, e.g., with an initial guess Σ = Σ of the self-energy, the Green function G is computedusing the lattice Dyson equation (5). In a second step, Σand G yield the bath Green function G via the impurityDyson equation (4), which defines, in combination withthe local interactions, the impurity problem [illustratedin Fig. 1(b)], the solution of which is the nontrivial partof the algorithm. A second application of the impurityDyson equation (4), to the resulting G and to G , yieldsa new estimate of the self-energy Σ, which closes theself-consistency cycle. In the following, we discuss theprimary options for addressing the impurity problem. Direct impurity solvers – One class of methods directlyevaluates the path integral representation of the Greenfunction [Eqs. (3) and (2)] for a continuous bath G ,which corresponds to a DMFT solution of the originallattice problem in the thermodynamic limit after self-consistency. We will refer to such methods as “directimpurity solvers.”For a long time, the Hirsch-Fye QMC (HF-QMC)algorithm has been the method of choice for nonper-turbative DMFT calculations [16]. HF-QMC is basedon a discretization of the imaginary time τ ∈ [0 , β ]into Λ “time slices” of width ∆ τ = β/ Λ, a Trotterdecomposition of the interaction and kinetic terms in Eq.(2), and a Hubbard-Stratonovich transformation, whichreplaces the electron-electron interaction by an auxiliarybinary field on each time slice; the resulting problem isthen solved employing Wick’s theorem and Monte Carloimportance sampling over the field configurations. Asconfigurations can be updated in the case of a singlespin flip (i.e., an auxiliary-field change on a single timeslice) with a matrix-vector operation of cost O (Λ ) andΛ local updates are needed for a global configurationupdate, the numerical cost of the HF-QMC algorithmscales as Λ . All HF-QMC results have statistical errors (which decay as N − / for N “sweeps,” each consisting ofΛ attempted single-spin updates) and systematic errorsresulting from the Trotter decomposition. As ∆ τ has tobe kept constant for roughly constant systematic errorupon variation of T , the numerical effort of HF-QMCscales as the cube of the inverse temperature, β . This isalso true for the numerically exact (unbiased) “multigrid”HF-QMC method [34].The integration of the (conventional) HF-QMCmethod into the DMFT self-consistency cycle is illus-trated in Fig. 2(b) [as a specification of the lower dashedbox in Fig. 2(a)]: a fixed choice of ∆ τ (diamond-shapedselection box) defines the grid τ l = l ∆ τ with 0 ≤ l ≤ Λfor a Fourier transform (square box) of the Matsubarabath Green function G ( iω n ) (with | ω n | ≤ ω max for somecutoff frequency ω max ) to the imaginary-time equivalent {G ( τ l ) } Λ l =0 . After application of the HF-QMC algorithm(rounded box), the result { G ( τ l ) } Λ l =0 is transformed backto Matsubara frequencies (square box); this step requiresspecial care in order to get around the Nyquist theorem,e.g. using analytic weak-coupling results [23, 35, 36].More recently, conceptionally different QMC ap-proaches have been formulated, which are based ondiagrammatic expansions of the action (2) in continuousimaginary time, either in the interaction U (CT-INT [37])or in the bath hybridization (CT-HYB [25, 26]), and ona stochastic sampling of Feynman diagrams; CT-AUX[38] is related to the HF-QMC method [27]. All of thesecontinuous time (CT-QMC) algorithms require Fouriertransforms, before and (with the exception of CT-INT)after the QMC part; the numerical cost is associatedprimarily with matrix updates, similar to those arisingin HF-QMC, with a total scaling of the computationaleffort, again, as β .Thus, all direct QMC-based impurity solvers are verycostly at low T , which limits their access to low-temperature phases of particular physical interest. Auxiliary Hamiltonian and exact diagonalization –Another class of numerical approaches, such as the “exactdiagonalization” methods, cannot directly be applied tothe action-based formulation of the impurity problem,but requires a Hamiltonian representation [39]. Onepossibility is the “star topology” illustrated in Fig. 1(c),where a central “impurity” site (with the same inter-actions as the impurity problem, here U ) is connectedby hopping matrix elements V iσ to a number N b ofnoninteracting “bath sites,” each characterized by a localpotential ǫ iσ . In general, this representation has to bespin-dependent, leading to the Anderson Hamiltonian H And = ǫ X σ n σ + U n ↑ n ↓ + X σ N b X i =1 h ǫ iσ n iσ + V iσ (cid:16) a † iσ c σ + h.c. (cid:17)i , (6)where c † σ ( c σ ) creates (annihilates) an electron with spin σ ∈ {↑ , ↓} on the impurity site and a † iσ ( a iσ ) creates(annihilates) an electron with spin σ on bath site i ; n σ = c † σ c σ , n iσ = a † iσ a iσ are the corresponding numberoperators. In this work, we consider only nonmagneticphases, which implies spin symmetric bath parameters V i ↓ = V i ↑ , ǫ i ↓ = ǫ i ↑ .For a fixed choice of N b , the bath parameters ǫ iσ , V iσ are determined such that the noninteracting impurityGreen function G And ,σ associated with H And , with G − ,σ ( ω ) = ω + µ σ − N b X i =1 V iσ ω − ǫ iσ , (7)is “close” to the target bath Green function G accordingto some metric (see below). Note that the resultingspectrum − π Im G And ,σ ( ω + i + ) is necessarily discrete(for finite N b ), in contrast to piece-wise smooth spectrumof the true bath Green function; in this sense, themapping to a Hamiltonian implies a “bath discretization”in frequency space; this step clearly introduces a biaswhich has to be controlled particularly carefully withiniterative procedures such as the DMFT.The integration of this type of approach in the DMFTcycle is illustrated for the case of ED in Fig. 2(c): for afixed choice of N b (diamond shaped box), the parameters V i , ǫ i (here and in the following we suppress spin indices)are adjusted (rounded box) as to minimize the bath misfit χ [ { V i , ε i } ] = n c X n =0 w n |G And ( iω n ; { V i , ε i } ) − G ( iω n ) | ,(8)with a cutoff Matsubara frequency iω n c and the weight-ing factor w n , which can be used to optimize the bathparametrization [40] and which we set to w n = 1 [41]. Asthis fit is performed directly on the Matsubara axis, noFourier transform is needed for G . Using ED (roundedbox), the Green function G can be evaluated on theMatsubara axis [29]; therefore, the DMFT cycle is closedwithout any Fourier transform.The minimization of χ [as defined in Eq. (8)] isperformed in our ED and BSS-DMFT calculations using the Newton method, based on analytic expressions for thederivative ∇ χ with respect to the bath parameters. Dueto the multidimensional character of the problem, thisdeterministic method is often trapped in local minima;thus, a naive implementation of Newton based methodswill, in general, not find globally optimal parameters,which can induce unphysical fixed points in the DMFTiteration procedure. Therefore, we use not only thesolution { V i , ǫ i } of the previous iteration as initialization,but perform a large number (up to 1000) of independentNewton searches, starting also from random initial pa-rameters. Of the resulting locally optimal solutions, wechoose the one with minimum χ as the final result ofthe minimization procedure; typically, about 1% of theindividual searches come close to this (estimated) globaloptimum.An advantage of ED, compared to QMC algorithms,is that Green functions and spectra can be computeddirectly on the real axis, without analytic continuation;however, numerical broadening of the resulting discretepeaks is required. This discretization problem is partic-ularly severe as the numerical effort of the matrix diago-nalization scales exponentially with the total number ofsites (here N b + 1), which limits the applicability of EDfor cluster extensions of DMFT or multiband models. B. Principles of the BSS-QMC algorithm andapplication as a DMFT impurity solver
In Eq. (6), we have used the conventional notation forthe auxiliary Hamiltonian that emphasizes its interpre-tation as an impurity model, e.g., with different creationoperators for electrons on the central “impurity” site( c † σ ) and on the bath sites ( a † iσ ), respectively. However,with the changes c σ → c σ , n σ → n σ , a † iσ → c † iσ ,and V iσ → t σ i , it essentially reproduces the Hubbardmodel (1) on a graph, just with nonuniform interaction( U acting only on site 0) and, possibly, spin-dependenthopping amplitudes and local energies.As a consequence, the model (6) is not only treatablewith the universal ED approach, but also with morespecific methods developed for Hubbard-type models. Aspointed out recently by Khatami et al. [31], this includesthe determinantal quantum Monte Carlo approach byBlankenbecler, Scalapino, and Sugar [32, 42] (BSS-QMC), which, thereby, becomes applicable as a DMFTimpurity solver. In the following, we will first sketchthe established BSS-QMC approach (for an extendeddiscussion, including issues of parallelization, see Ref. 43)and then discuss its application in the DMFT context.Similarly to the HF-QMC method (cf. Sec. II A),the BSS-QMC approach is based on a Trotter-Suzukidecomposition, here of the partition function Z = Tr (cid:16) e − β ( H K + H V ) (cid:17) (9) ≈ Z ∆ τ = Tr Λ Y l =0 e − ∆ τH K e − ∆ τH V ! , (10)where H V ( H K ) corresponds to the interaction (kineticand local potential) contribution to the Hubbard typemodels (1) or (6) and ∆ τ = β/ Λ. Again, a discreteHubbard-Stratonovich transformation replaces the inter-action term by a binary auxiliary field { h } with h i ( l ) = ± i and time slice l . The trace in Eq. (10)then simplifies to Z ∆ τ = X { h } det h M { h }↑ i det h M { h }↓ i with (11) M { h } σ = + B Λ ,σ (cid:2) { h i (Λ) } Ni =1 (cid:3) . . . B ,σ (cid:2) { h i (1) } Ni =1 (cid:3) ,where B is defined in terms of the hopping matrix K : B l,σ (cid:2) { h i ( l ) } Ni =1 (cid:3) = e σλ diag[ h ( l ) ,...,h N ( l )] e − ∆ τK . (12)The interaction strength is encoded in the parameter λ =cosh − ( e U ∆ τ/ ). The computation of thermal averages ofphysical observables O takes the form: h O i = X { h } h O { h } P { h } ∆ τ i (13) P { h } ∆ τ = 1 Z ∆ τ det h M { h }↑ i det h M { h }↓ i .At particle-hole symmetry, the weights P { h } ∆ τ are alwayspositive; i.e., the sums can be evaluated at arbitraryprecision, without any sign problem. As in HF-QMC, theproblem is solved by Monte Carlo importance samplingof the auxiliary field { h } and evaluation of the Greenfunction at time slice l , with G { h } l,σ = [ + B l − ,σ . . . B ,σ B Λ ,σ . . . B l,σ ] − . (14)As a spin flip in the auxiliary field h i ( l ) at time slice l and site i only affects B l,σ at this site, the ratio ofthe weights, needed for the decision whether a proposedspin flip is accepted, involves only local quantities; a fullrecomputation of the determinants of N × N matricesappearing in Eq. (11) is not needed. The computationaleffort is further reduced by calculating the Green functionat time slice l +1 from the quantities at time slice l , usingso-called “wrapping”: G l +1 ,σ = B − l,σ G l,σ B l,σ . (15)In order to avoid the accumulation of numerical errors inthe matrix multiplications, it is necessary to recalculatethe full Green function at regular intervals. This isparticularly important at low temperatures.All this considered, the numerical cost scales cubicallywith the number of sites and linearly with the number of time slices; at constant ∆ τ , this translates to a totaleffort O ( N β ), where N = N b + 1. Note that a needfor finer bath discretizations at lower temperatures couldpotentially spoil the scaling advantage of the method overdirect impurity solvers; we will show in Sec. III C thatthis is not the case for our test applications.The application in the DMFT context [31] starts withthe computation of the Hamiltonian parameters (forsome choice of N b ), exactly like in the ED approach. Asin the HF-QMC approach, one then chooses some dis-cretization ∆ τ , computes { G ( l ∆ τ ) } Λ l =0 for the impuritysite, and applies a (nontrivial) Fourier transform back toMatsubara frequencies. The result is an impurity solverwith superior scaling (linear in β ) compared to the directimpurity solvers (cubic in β ), however with a bias due tothe Trotter discretization ∆ τ (in addition to a possiblebias due to the bath discretization with N b sites) which,as we will show in Sec. III, can be quite significant. C. Specification of multigrid BSS-QMC algorithm
The central feature of our new algorithm is the elim-ination of this systematic Trotter error, while retainingthe advantage of linear-in- β scaling inherent in the BSS-QMC method. In the following, we will specify themethod and illustrate it using an example (Fig. 3), thatwill be discussed in detail in Sec. III.In contrast to the previous DMFT-BSS implementa-tion with a unique discretization ∆ τ in all BSS com-putations throughout the DMFT self-consistency cycle,the (impurity) Green function of the Hamiltonian H And at hand is computed in M &
20 parallel BSS runs(indexed by 1 ≤ i ≤ M ), each employing a homogeneousimaginary-time grid with a specific discretization (∆ τ ) i ,chosen from a set { (∆ τ ) i | (∆ τ ) min ≤ (∆ τ ) i ≤ (∆ τ ) max } with typically 6 - 9 different elements. Green functionsresulting from BSS-QMC runs with the same discretiza-tion (∆ τ ) i = (∆ τ ) j are averaged over, thereby reducingthe dependencies on initialization conditions and furtherenhancing the parallelism.This leads to a set of Green functions defined, in gen-eral, on incommensurate imaginary-time grids (symbolsin Fig. 3). In order to apply a local ∆ τ → G (∆ τ ) i have to be transformed to a commongrid. This is possible since the true G ( τ ) is a smoothfunction; however, a direct spline interpolation of theraw QMC results, neglecting higher derivatives, wouldnot be accurate [44]. Instead, we consider differencesbetween the raw data { G (∆ τ ) i ( l ∆ τ ) } Λ i l =0 and a referenceGreen function, obtained via Eq. (5) from a model self-energy Σ ref σ ( iω n ) [35, 36], written here for the single-bandcase (for multi-band generalizations, see Ref. [45]):Σ ref σ ( iω n ) = U (cid:18) h n − σ i − (cid:19) + 12 U h n − σ i (1 − h n − σ i ) × (cid:18) iω n + ω + 1 iω n − ω (cid:19) , (16)which recovers the exact high-frequency asymptotics ofΣ( iω n ) and G ( iω n ) for any choice of the free parameter ω and, therefore, approximates the second and higher-order derivatives of G ( τ ) at τ → τ → β ) well. Thismatch can be further improved by adjusting ω . Con-sequently, the differences { G (∆ τ ) i ( l ∆ τ ) − G ref ( l ∆ τ ) } Λ i l =0 have smaller absolute values and much smaller higherderivatives than the original data; in particular, theircurvature vanishes asymptotically at the boundaries [46].Thus, they are well represented by natural cubic splines.Usually, the parameters of the piecewise polynomialsconstituting such a spline f spline ( x ) are determined fromdiscrete data { f meas ( x i ) } Ni =0 such that the discrete dataare reproduced exactly: f spline ( x i ) = f meas ( x i ) for all 0 ≤ i ≤ N . However, in the QMC context, all measurementshave statistical errors, i.e., the discrete data are betterrepresented as { f meas ( x i ) ± ∆ f meas ( x i ) } Ni =0 with standarddeviations ∆ f meas , which are also estimated within theQMC procedure. It is clear that the usual interpolatingsplines, which do not take the uncertainties of the discretedata into account, contain more features than warrantedby the data (in particular at the Nyquist frequency); inthe context of Green functions this includes the possi-bility of acausal behavior. We use, instead, smoothingspline fits [47, 48] which reproduce the discrete dataonly within error bars, which are typically O (10 − ), (andminimize the curvatures under this constraint); these fitscan be computed in a very similar procedure and at thesame cost as interpolating splines.After combining these approximating “difference”splines with exact expressions for G ref ( τ ) resulting fromEqs. (16) and (5), we obtain smooth approximations ofthe Green functions, as seen in Fig. 3(a); the inset alsodemonstrates slight deviations of the continuous splinefits from the discrete data (within error bars), e.g., for thediscretization ∆ τ = 0 . τ ≈ . τ fine = 0 . τ →
0. This is illustrated in Fig.4 for the representative values of τ denoted by arrowsin Fig. 3(a). Even though most of the raw BSS-QMCdata do not include estimates of the Green functions atthese precise values of τ , the transformed data (symbolsin Fig. 4) depend very regularly on ∆ τ , falling on nearlystraight lines as a function of (∆ τ ) . Therefore, theycan accurately and reliably be extrapolated to ∆ τ → τ = 0); an application ofthis procedure at all τ (on the fine grid) leads to quasi-continuous Green functions without significant Trottererrors, shown as solid lines in Fig. 3. These resultscan be Fourier transformed to Matsubara frequenciesin a straightforward manner [cf. Fig. 1(d)]. A similarapproach has also been useful for computing unbiasedspectra from BSS-QMC [49].At first sight, the computational advantage of themultigrid procedure is less obvious in the BSS-QMC G ( τ ) T = 0.04, U = 4.40 G ( τ ) τ G ( τ ) τ T = 0.04, U = 5.10 ∆τ=0.4∆τ=0.7 ∆τ=0.9∆τ → 0 ED G ( τ ) τ (a)(b) FIG. 3. (Color online) BSS-QMC impurity Green functions at T = 0 .
04 (symbols) using a bath representation with N b = 4sites (with parameters of converged DMFT-ED solution, long-dashed lines) and results of multigrid extrapolation to ∆ τ = 0(solid lines). Upper panel: metallic phase ( U = 4 . U = 5 . τ valuesfor which the extrapolation is shown in Fig. 4. context than for HF-QMC [34, 50], since the numericaleffort for direct computations at small ∆ τ grows onlylinearly, not cubically, with (∆ τ ) − in the BSS case(while the systematic errors decay generically as (∆ τ ) for a given impurity problem). However, even for afixed Hamiltonian, so much accuracy can be gained byextrapolation that it more than offsets the cost of theadditional grid points. This is true, in particular, sincestable results are best obtained by averaging over inde-pendent BSS-QMC runs; performing these on variablegrids then allows for extrapolation without additionalcost. Furthermore, the individual runs thermalize fasterin the multigrid variant, due to the smaller numberof time slices (and proportionally shorter run time persweep), which enhances the parallelism. Most impor-tantly, as we will see below, the DMFT self-consistencycan magnify any bias of the employed impurity solvers incomplicated ways (in the vicinity of phase transitions), sothat controlled results are really dependent on unbiasedmethods, such as our multigrid approach. G ( τ ) ( ∆τ ) U = 4.40, T = 0.04 τ = 1.0τ = 2.0 τ = 5.0τ = 12.5
FIG. 4. (Color online) BSS-QMC estimates of imaginary-time Green functions G ∆ τ ( τ ) at T = 0 . U = 4 . τ (symbols) and extrapolation to∆ τ = 0 using least-squares fits (lines). III. RESULTS
In this section, we compare results of the new nu-merically exact “multigrid” BSS-QMC method with rawBSS-QMC results (at finite Trotter discretization), withreference ED results (which are exact at the level ofthe auxiliary Hamiltonian), and with the predictions ofestablished impurity solvers (multigrid HF-QMC [34] andCT-HYB [25, 26, 51]). These comparisons are performedin three stages: In Sec. III A, we keep the bath G andits approximation by an auxiliary Hamiltonian fixed anddiscuss the impact of the Trotter error and its eliminationwithout the complications of the DMFT self-consistency.In Sec. III B, we compare full DMFT solutions obtainedusing the various algorithms at moderate temperature( T = 0 . T ≥ .
01 (with DMFT self-consistency), where the impact of the bath discretizationbecomes particularly relevant, in Sec. III C.Following the established practice for the evaluation ofDMFT impurity solvers [23, 28], all of these comparisonsare performed for the half-filled Hubbard model withsemi-elliptic “Bethe” density of states [52] (full bandwidth W = 4) within the paramagnetic phase. Specifi-cally, we choose temperatures T ≤ .
04, which are belowthe critical temperature T ∗ ≈ .
055 [36, 50] of the first-order metal-insulator transition, and interactions close toor within the coexistence region of metallic and insulatingsolutions, which arises from the mean-field character ofthe DMFT.
A. Green function extrapolation at fixed bathHamiltonian parameters
In general, a bias present in an impurity solver has atwo-fold impact: On the one hand, it affects estimates of Green functions and all other properties for a givenimpurity problem, defined by its bath Green function G .On the other hand, it shifts the fixed point of the DMFTself-consistency cycle, i.e., it also modifies the convergedbath Green function, which, in turn, also affects themeasured Green functions and all other properties. Inthis subsection, we study the first effect in isolationby fixing the bath Green function to the convergedsolution of the ED procedure for N b = 4 bath sites (withHamilton parameters { ǫ i , V i } i =1 ). As the same auxiliaryHamiltonian is used also in the BSS-QMC algorithm, theED estimates of the Green function are exact for thepurpose of the current comparison; the impact of thebath discretization (which corresponds to a bias on theDMFT level) will be discussed in Sec. III C.Local imaginary-time Green functions G ( τ ) are shownin Fig. 3(a) for the metallic phase, at U = 4 .
4, and inFig. 3(b) for the insulating phase, at U = 5 .
1. Here andin the following, we restrict the imaginary-time rangeto 0 ≤ τ ≤ β/
2; data for τ > β/ G ( β − τ ) = G ( τ ). Symbols (inthe magnified insets) represent raw BSS-QMC results(with discretizations ∆ τ = 0 .
4, ∆ τ = 0 .
7, and ∆ τ =0 . τ ≈
2. In contrast, multigrid BSS-QMCGreen functions (solid lines) are indistinguishable fromthe ED data at U = 5 . U = 4 .
4, with deviations of the order of statistical errors.Thus, our method yields, indeed, quasi-continuous Greenfunctions without significant Trotter errors in both testcases, although the discretizations of the underlying rawBSS-QMC computations (with 0 . ≤ ∆ τ ≤ .
0) would beconsidered much too coarse in conventional applications.A very similar picture emerges in an analogous compar-ison for the two coexisting solutions at U = 4 .
74, shownin Fig. 5. Again, the raw BSS-QMC results (symbolsand colored broken lines) show a strong systematic bias,towards more metallic Green functions and of differentmagnitude in the different phases, while the extrapolatedGreen functions agree nearly perfectly with the EDreferences. In fact, some of the BSS Green functionscalculated for an insulating bath (lower set of symbolsand broken lines) show such large discretization errorsat small τ .
2, that they approach the exact Greenfunction of the metallic DMFT solution (upper solidand long-dashed lines). One may suspect from thisobservation that these biased “insulating” solutions willnot be associated with stable DMFT fixed points if theyare fed back in the self-consistency cycle; such shifts ofstability regions induced by the Trotter bias at ∆ τ > { (∆ τ ) i } of discretizations G ( τ ) τ metallicinsulating ∆τ=0.4∆τ=0.7 ∆τ=0.9∆τ → 0 ED G ( τ ) τ T = 0.04, U = 4.74 metallicinsulating
FIG. 5. (Color online) BSS-QMC impurity Green functionsat T = 0 .
04 and U = 4 .
74 (symbols and colored brokenlines) using a bath representation with N b = 4 sites (withparameters fixed by converged DMFT-ED solution, long-dashed lines) and extrapolation to ∆ τ = 0 (solid lines).Upper (lower) set of curves: metallic (insulating) bath. G ( τ ) τ T = 0.04, U = 4.70insulating G ( τ ) τ FIG. 6. (Color online) Green functions in the insulatingphase at T = 0 . U = 4 .
7, extrapolated from BSS-QMCresults using different imaginary-time grids [53] (at fixed bathrepresentation). The excellent agreement shows that themultigrid procedure is stable with respect to its technicalparameters. in the underlying BSS-QMC runs, i.e., if no sensiblechoice leads to a significant bias. This is demonstratedin Fig. 6 for the insulating phase at U = 4 .
7: the Greenfunctions for the same auxiliary problem obtained frommultigrid extrapolations with three different ∆ τ grids[53] agree perfectly within the precision of the method.The latter is primarily determined by the statistical errors, i.e., by the number of sweeps and, possibly,by the numerical precision in the matrix operations.Only if raw BSS-QMC data of much higher precisionwas available (with many millions of sweeps per run),additional accuracy could be gained by choosing smallerdiscretizations (e.g., 0 . ≤ ∆ τ ≤ . τ that are 3 to 10 times as large as thediscretization that one would choose in a conventionalBSS-QMC procedure. B. Comparisons of impurity solvers at full DMFTself-consistency: impact of Trotter errors
So far, we have compared different algorithms justat the impurity level, i.e., for a fixed bath Greenfunction (determined from a self-consistent DMFT-EDcalculation). In contrast, we will now discuss results ofcompletely independent DMFT solutions, each of whichcorresponds to full self-consistency for a given impuritysolver (cf. Fig. 2). For all Hamiltonian based methods(ED, BSS-QMC, and multigrid BSS-QMC), the numberof bath sites is restricted to N b = 4 (as above); the impactof this parameter will be studied in Sec. III C.Specifically, we discuss static observables that areparticularly useful for discriminating between metallicand (possibly coexisting) insulating DMFT solutions,namely the double occupancy D = h n ↑ n ↓ i , (17)which is proportional to the interaction energy E int = U D , and the quasiparticle weight Z = (cid:20) − ∂ Re Σ( ω ) ∂ω (cid:12)(cid:12)(cid:12) ω =0 (cid:21) − ≈ (cid:20) iω ) πT (cid:21) − . (18)Open symbols in Fig. 7 denote estimates resulting fromself-consistent DMFT solutions using the conventionalBSS-QMC impurity solver at finite discretization 0 . ≤ ∆ τ ≤ .
5, i.e., using the scheme established in Ref. [31].The estimated values of Z , shown in Fig. 7(b), have anearly uniform offset in the metallic phase at U . . U ≈ . τ = 0 .
4) than the ED reference solution.This is also seen in corresponding estimates of thedouble occupancy [Fig. 7(a)]; however, for these ob-servables the Trotter bias is highly nonuniform (in themetallic solution): at U = 4 . U & . U . .
5. At the same time, nearly all D T = 0.04, N b = 4 Z U (a)(b) multigrid HFBSS ( ∆τ=0.5 )BSS ( ∆τ=0.4 )BSS ( ∆τ=0.3 )BSS ( ∆τ→ 0 )ED FIG. 7. (Color online) Estimates of double occupancy D ( U )and quasiparticle weight Z ( U ) obtained in independent self-consistent DMFT calculations using various impurity solvers:multigrid HF-QMC (crosses), conventional BSS-QMC (opensymbols), multigrid BSS-QMC (circles), and ED (diamonds).In each panel, the upper (lower) sets of curves correspondto metallic (insulating) solutions. Lines are guides to theeye only. Arrows in (a) indicate parameters for which thediscretization dependence is studied in Fig. 8. of these data deviate significantly (and without obvioussystematics) from the reference ED result (diamonds), sothat an a posteriori elimination of the Trotter bias seemsimpossible.In contrast, the new multigrid BSS-QMC procedure, asdiscussed in Sec. II C and illustrated in Fig. 2(d), leads toestimates of both D and Z (filled circles) which perfectlyrecover the ED solutions, even though they are based onBSS-QMC runs with ∆ τ ≥ . D andvery nonuniform for Z . Again, the multigrid BSS-QMCresults agree perfectly with the ED reference data.For comparison, crosses and black solid lines in Fig. 7denote estimates of an unbiased direct impurity solver,namely the multigrid HF-QMC method [34]; these showgood overall agreement with both the ED and themultigrid BSS-QMC data. A slight negative deviationin the estimates of D of the latter, Hamiltonian based,methods can be traced back to the relatively poor bathdiscretization with N b = 4 auxiliary sites (cf. Sec. III C).Since the double occupancy D is best computed di-rectly on the impurity level (in QMC based approaches), D ( ∆τ ) T = 0.04, N b = 4 U = 4.7U = 5.0 FIG. 8. (Color online) Discretization dependence of thedouble occupancy D as estimated from BSS-QMC, usingeither the multigrid scheme (filled symbols) or self-consistentBSS-QMC solutions at finite ∆ τ (open symbols), within themetallic phase at U = 4 . U = 5 . τ ≥ .
3; dotted lines denote fits that include alsodata at ∆ τ = 0 . τ = 0 . its physical value has to be extrapolated from rawestimates D ∆ τ , with discretizations corresponding to thedifferent grid points used within the multigrid procedure(in contrast to the quasiparticle weight Z , which followsfrom the self-energy Σ, which, in turn, is determined fromunbiased Green functions).As seen in Fig. 8, the Trotter bias inherent in these rawestimates (filled symbols) is perfectly regular [54] even atlarge ∆ τ , so that reliable extrapolations ∆ τ → U = 4 . U = 5 . D resulting from conventionalBSS-QMC calculations in the same range of discretiza-tions ∆ τ ≥ . τ that quadratic least-square fits (solid lines) lead to extrapolations ∆ τ → within the self-consistency cycle, as introducedby our multigrid approach, can efficiently generate high-precision results. C. Comparisons of impurity solvers at full DMFTself-consistency: impact of bath discretization
So far, we have restricted the bath representation inall Hamiltonian based impurity solvers (ED and bothvariants of BSS-QMC) to only N b = 4 bath sites andfocused on the impact of the Trotter errors and theirelimination. From the mutual agreement with multigridHF-QMC, an impurity solver which treats the bathdirectly on the action level, we can conclude that this0 D T = 0.04 multigrid HFmultigrid BSS N b =6ED N b =5ED N b =4ED N b =30.00.1 4.5 4.6 4.7 4.8 4.9 5.0 5.1 Z U (a)(b) FIG. 9. (Color online) Estimates of double occupancy D ( U )and quasiparticle weight Z ( U ) at T = 0 .
04, obtained inselfconsistent DMFT calculations using Hamiltonian basedimpurity solvers with 3 ≤ N b ≤ N b → ∞ . In eachpanel, the upper (lower) sets of curves correspond to metallic(insulating) solutions. Lines are guides to the eye only. coarse bath discretization allows for reasonably accurateestimates of D and, in particular, Z at the moderatelylow temperature T = 0 .
04. However, the ED andmultigrid BSS-QMC estimates of D were found in Fig.7 to lie a bit below the multigrid HF-QMC data; thisdeviation must be an artifact of the bath discretizationif the multigrid HF-QMC reference data are correct.Moreover, we must suspect that the bath discretizationbias gets worse (at constant N b ) at lower temperatures.Figure 9 shows estimates of D ( U ) and Z ( U ) at T =0 .
04, similarly to Fig. 7 and with the same multigridHF-QMC reference data (crosses), but now using Hamil-tonian based impurity solvers with 3 ≤ N b ≤ N b = 6). Smaller bathsizes ( N b = 3, N b = 4, and N b = 5) are representedonly by the ED solution, which is cheaper and free ofstatistical noise. At the resolution of Fig. 9, the estimatesassociated with the finer bath discretizations N b = 5(triangles) and N b = 6 (circles) agree with each other.Therefore and since they are also consistent with theunbiased multigrid HF-QMC data (crosses), we concludethat convergence with respect to the bath discretizationis reached already at N b = 5 at T = 0 .
04. In contrast,the ED estimates of D are apparently slightly too small D T = 0.02 multigrid HFmultigrid BSS N b =6ED N b =5ED N b =40.00.1 4.5 4.6 4.7 4.8 4.9 5.0 5.1 Z U (a)(b) FIG. 10. (Color online) Estimates of double occupancy D ( U )and quasiparticle weight Z ( U ) at T = 0 .
02, using bathdiscretizations with 4 ≤ N b ≤ at N b = 4 (pentagons); corresponding results at N b = 3(squares) are far off both for D and Z .Note that consistent convergence of observable esti-mates with N b , as demonstrated in Fig. 9 (as well asFig. 10 and Fig. 11) can only be observed when optimalHamiltonian parameters are determined with great care,as described in Sec. II A, within each self-consistencycycle; otherwise some bath sites may remain ineffective orthe estimates can even get worse upon increasing N b . Inaddition (as always in the DMFT context), it is essentialthat enough DMFT iterations are performed at eachphase point in order to ensure convergency with respectto the self-consistency cycle (cf. Fig. 2).Halving the temperature amplifies the bath discretiza-tion effects, as seen in Fig. 10: At T = 0 .
02, onlythe best Hamiltonian representation ( N b = 6, evaluatedwith multigrid BSS-QMC, circles) recovers all referencemultigrid HF-QMC results (crosses) within their accu-racy. At N b = 5, the estimates of D in the insulatingphase are already slightly too small; at N b = 4, strongnegative deviations in D ( U ) are apparent also for themetallic solution. The impact of the bath discretizationbecomes even much stronger at T = 0 .
01, as shown inFig. 11 for the metallic phase, which is interesting as astrongly renormalized Fermi liquid (while the propertiesof the insulating phase are asymptotically independent oftemperature). We find that even the results for N b = 5and N b = 6 deviate significantly in Fig. 11(a) from eachother and from the multigrid HF-QMC reference result,especially at U = 5 .
3, near the edge of the stabilityregion of the metallic phase. Only the multigrid BSS-QMC results using N b = 7 bath sites (circles) agree1 D T = 0.01
CT-HYBmultigrid HFmultigrid BSS N b =7ED N b =6ED N b =5ED N b =40.000.050.105.0 5.1 5.2 5.3 Z U (a)(b) FIG. 11. (Color online) Estimates of double occupancy D ( U )and quasiparticle weight Z ( U ) at T = 0 .
01, using bathdiscretizations with 4 ≤ N b ≤ N b → ∞ . with the reference data within their precision; even betteragreement is observed with data obtained using the CT-HYB impurity solver (diamonds). Note that BSS-QMCis much more efficient than ED at N b = 7; in particular,the latter would need orders of magnitude more mainmemory than the former.As noted in Sec. II B, a strong increase with inversetemperature of the number N b of bath sites neededfor a given accuracy could, in principle, eliminate thescaling advantage of the DMFT-BSS approach, as itscomputational cost does not only include the directfactor β , but also a factor N b ≡ N b ( β ). However,our results indicate that this effect is minor: In ourtest case, we needed to add one bath site upon halvingthe temperature for roughly constant accuracy; this isconsistent with a scaling N b ∝ ln( β ), i.e., an overallcomputational cost proportional to β [ln( β )] which isstill linear up to logarithmic corrections. IV. CONCLUSIONS
The DMFT and its extensions are invaluable tools forthe study of phenomena associated with strong electroniccorrelations and for quantitative predictions of propertiesof correlated materials. However, the numerical solutionof the DMFT self-consistency equations remains a greatchallenge: the established, direct, QMC impurity solversyield unbiased results, but provide only limited accessto the low- T phase regions of interest, due to the cubicscaling of their computational cost with the inversetemperature β . Exact diagonalization (ED) approaches,on the other hand, are limited by their exponentialscaling with the number of sites N of the auxiliaryHamiltonian.The multigrid BSS-QMC algorithm presented in thiswork allows for solving the DMFT self-consistency equa-tions with an effort that grows only linearly with β ;in contrast to an earlier BSS-QMC based method [31],it is free of significant Trotter errors, i.e., numericallyexact at the level of the auxiliary Hamiltonian. Since thecomputational cost grows only cubically with N , muchbetter representations of the bath are possible than forED. As demonstrated by applications to the half-filledHubbard model in and near the coexistence region ofmetallic and insulating solutions and by comparisonswith direct QMC impurity solvers, the new method yieldsunbiased results (for sufficiently fine bath discretization),in spite of using quite coarse Trotter discretizations in theunderlying BSS-QMC evaluations.The new unbiased quasi-CT impurity solver shouldshow its full potential in multi-band cases and in clusterextensions of DMFT, where the prefactor N of the BSS-QMC scheme (compared to a factor of 1 in HF-QMCcalculations in single-site DMFT) is leveled off by theincreased complexity of the original DMFT problem. Ourapproach can also be extended beyond Hubbard models;it could be particularly valuable for the cellular DMFTtreatment of the Kondo lattice model, where interestingtemperature regimes are out of reach of the existingimpurity solvers [26, 55].We thank E. Gorelik, E. Khatami, R. Scalettar, andP. Werner for valuable discussions. Financial support bythe Deutsche Forschungsgemeinschaft through FOR 1346is gratefully acknowledged. [1] W. Metzner and D. Vollhardt, Phys. Rev. Lett. , 324(1989).[2] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozen-berg, Rev. Mod. Phys. , 13 (1996).[3] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. ,865 (2006).[4] D. Vollhardt, Ann. Phys. , 1 (2012). [5] G. Kotliar, S. Y. Savrasov, G. P´alsson, and G. Biroli,Phys. Rev. Lett. , 186401 (2001).[6] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler,Rev. Mod. Phys. , 1027 (2005).[7] K. Held, I. A. Nekrasov, G. Keller, V. Eyert, N. Bl¨umer,A. K. McMahan, R. T. Scalettar, T. Pruschke, V. I.Anisimov, and D. Vollhardt, Phys. Stat. Sol. B ,2599 (2006). [8] S. Biermann, F. Aryasetiawan, and A. Georges, Phys.Rev. Lett. , 086402 (2003).[9] K. Held, C. Taranto, G. Rohringer, and A. Toschi,in The LDA+DMFT approach to strongly correlatedmaterials , edited by E. Pavarini, E. Koch, D. Vollhardt,and A. Lichtenstein (Forschungszentrum J¨ulich GmbH,Institute for Advanced Simulations, 2011).[10] C. Taranto, M. Kaltak, N. Parragh, G. Sangiovanni,G. Kresse, A. Toschi, and K. Held, arXiv:1211.1324.[11] J. M. Tomczak, M. Casula, T. Miyake, F. Aryasetiawan,and S. Biermann, Europhys. Lett. , 67001 (2012).[12] G. Kotliar and D. Vollhardt, Phys. Today , 53 (2004).[13] R. J¨ordens, L. Tarruell, D. Greif, T. Uehlinger,N. Strohmaier, H. Moritz, T. Esslinger, L. De Leo,C. Kollath, A. Georges, V. Scarola, L. Pollet,E. Burovski, E. Kozik, and M. Troyer, Phys. Rev. Lett. , 180401 (2010).[14] E. V. Gorelik, I. Titvinidze, W. Hofstetter, M. Snoek,and N. Bl¨umer, Phys. Rev. Lett. , 065301 (2010).[15] E. V. Gorelik, D. Rost, T. Paiva, R. Scalettar,A. Kl¨umper, and N. Bl¨umer, Phys. Rev. A ,061602(R) (2012).[16] M. Jarrell, Phys. Rev. Lett. , 168 (1992).[17] A. Georges and G. Kotliar, Phys. Rev. B , 6479 (1992).[18] V. Janiˇs and D. Vollhardt, Int. J. Mod. Phys. B , 731(1992).[19] J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. , 2521(1986).[20] N. Bl¨umer, in The LDA+DMFT approach to stronglycorrelated materials , edited by E. Pavarini, E. Koch,D. Vollhardt, and A. Lichtenstein (ForschungszentrumJ¨ulich GmbH, Institute for Advanced Simulations, 2011).[21] R. M. Fye and R. T. Scalettar, Phys. Rev. B , 3833(1987).[22] N. Bl¨umer and E. Kalinowski, Phys. Rev. B , 195102(2005); Physica B , 648 (2005).[23] N. Bl¨umer, Phys. Rev. B , 205120 (2007).[24] A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein,Phys. Rev. B , 035122 (2005).[25] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, andA. J. Millis, Phys. Rev. Lett. , 076405 (2006).[26] P. Werner and A. J. Millis, Phys. Rev. B , 155107(2006).[27] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,M. Troyer, and P. Werner, Rev. Mod. Phys. , 349(2011).[28] E. Gull, P. Werner, A. Millis, and M. Troyer, Phys. Rev.B , 235123 (2007).[29] M. Caffarel and W. Krauth, Phys. Rev. Lett. , 1545(1994).[30] E. Koch, in The LDA+DMFT approach to stronglycorrelated materials , edited by E. Pavarini, E. Koch,D. Vollhardt, and A. Lichtenstein (ForschungszentrumJ¨ulich GmbH, Institute for Advanced Simulations, 2011).[31] E. Khatami, C. R. Lee, Z. J. Bai, R. T. Scalettar, andM. Jarrell, Phys. Rev. E , 056703 (2010).[32] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,Phys. Rev. D , 2278 (1981).[33] E. Khatami, K. Mikelsons, D. Galanakis, A. Macridin,J. Moreno, R. T. Scalettar, and M. Jarrell, Phys. Rev.B , 201101 (2010).[34] N. Bl¨umer, arXiv:0801.1222; E. V. Gorelik andN. Bl¨umer, Phys. Rev. A , 051602 (2009). [35] C. Knecht, Numerische Analyse des Hubbard-Modellsim Rahmen der Dynamischen Molekularfeld-Theorie ,Master’s thesis, Johannes Gutenberg-Universit¨at Mainz(2003).[36] N. Bl¨umer,
Mott-Hubbard Metal-Insulator Transition andOptical Conductivity in High Dimensions , Ph.D. thesis,University of Augsburg (2002).[37] A. Rubtsov and A. Lichtenstein, Sov. Phys. JETP , 61(2004).[38] E. Gull, P. Werner, O. Parcollet, and M. Troyer,Europhys. Lett. , 57003 (2008).[39] E. Koch, G. Sangiovanni, and O. Gunnarsson, Phys.Rev. B , 115102 (2008).[40] D. S´en´echal, Phys. Rev. B , 235125 (2010).[41] In our experience, this choice for w n is not crucial: otherchoices lead to very similar results.[42] F. F. Assaad and H. G. Evertz, in Computational ManyParticle Physics , Lecture Notes in Physics, Vol. 739,edited by H. Fehske, R. Schneider, and A. Weiße(Springer Verlag, Berlin, 2008) p. 277.[43] C.-R. Lee, I.-H. Chung, and Z. Bai, in
Parallel Dis-tributed Processing (IPDPS) (2010) pp. 1 –9.[44] J. Joo and V. Oudovenko, Phys. Rev. B , 193102(2001).[45] C. Knecht, N. Bl¨umer, and P. G. J. van Dongen, Phys.Rev. B , 081103 (2005).[46] A reasonably good approximation of the curvature of G ( τ ) at the boundaries by the reference Green functionis the only crucial requirement for this procedure; it canalso be achieved using other model self-energies or, alter-natively, via maximum entropy analytic continuation ofthe “best” measured G (∆ τ ) i .[47] P. Craven and G. Wahba, Numer. Math. , 377 (1978).[48] I. G. Enting, C. M. Trudinger, and D. M. Etheridge,Tellus B , 305 (2006).[49] D. Rost, E. V. Gorelik, F. Assaad, and N. Bl¨umer, Phys.Rev. B , 155109 (2012).[50] N. Bl¨umer and E. V. Gorelik, Phys. Rev. B , 085115(2013).[51] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin,J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull,S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov,D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet,G. Paw lowski, J. D. Picon, L. Pollet, E. Santos, V. W.Scarola, U. Schollw¨ock, C. Silva, B. Surer, S. Todo,S. Trebst, M. Troyer, M. L. Wall, P. Werner, andS. Wessel, J. Stat. Mech. Theor. Exp. , P05001(2011).[52] M. Kollar, M. Eckstein, K. Byczuk, N. Bl¨umer, P. vanDongen, M. Radke De Cuba, W. Metzner, D. Tanaskovi´c,V. Dobrosavljevi´c, G. Kotliar, and D. Vollhardt, An-nalen der Physik , 642 (2005).[53] The results shown in Fig. 6 correspond to the followinggrids of the discretization ∆ τ : [0.29, 0.40, 0.50, 0.60,0.69, 0.78, 0.89, 0.96] (yellow line), [0.29 0.34 0.40 0.440.50 0.54 0.60] (light blue line), and [0.19 0.25 0.29 0.340.40 0.44 0.50 0.54 0.60] (dark blue line).[54] Specifically, we have used fit laws of the formln[ D (∆ τ )] = D (∆ τ = 0) + A (∆ τ ) .[55] L. C. Martin and F. F. Assaad, Phys. Rev. Lett. ,066404 (2008); L. C. Martin, M. Bercx, and F. F.Assaad, Phys. Rev. B82