Chebyshev Matrix Product State Impurity Solver for the Dynamical Mean-Field Theory
F. Alexander Wolf, Ian P. McCulloch, Olivier Parcollet, Ulrich Schollwöck
CChebyshev Matrix Product State Impurity Solverfor the Dynamical Mean-Field Theory
F. Alexander Wolf, Ian P. McCulloch, Olivier Parcollet, and Ulrich Schollw¨ock Theoretical Nanophysics, Arnold Sommerfeld Center for Theoretical Physics,LMU M¨unchen, Theresienstrasse 37, D-80333 M¨unchen, Germany Centre for Engineered Quantum Systems, School of Physical Sciences,The University of Queensland, Brisbane, Queensland 4072, Australia Institut de Physique Th´eorique, CEA, IPhT, CNRS, URA 2306, F-91191 Gif-sur-Yvette, France (Dated: July 8, 2014)We compute the spectral functions for the two-site dynamical cluster theory and for the two-orbitaldynamical mean-field theory in the density-matrix renormalization group (DMRG) framework us-ing Chebyshev expansions represented with matrix product states (MPS). We obtain quantitativelyprecise results at modest computational effort through technical improvements regarding the trun-cation scheme and the Chebyshev rescaling procedure. We furthermore establish the relation of theChebyshev iteration to real-time evolution, and discuss technical aspects as computation time andimplementation in detail.
I. INTRODUCTION
The dynamical mean-field theory (DMFT) and itscluster extensions are among the most successful meth-ods to study strongly correlated electron systems in di-mensions higher than one. The impurity problem withinDMFT is usually solved with continuous-time quantumMonte Carlo (CTQMC) algorithms, the numericalrenormalization group (NRG) or exact diagonalization(ED). While CTQMC is computationally feasibleeven for problems with many bands or a high numberof cluster sites, it provides numerically exact results onlyon the imaginary frequency axis. Many experimentallyrelevant frequency-dependent quantities like e.g. the con-ductivity therefore can only be obtained via the numer-ically ill-conditioned analytical continuation. NRG, bycontrast, solves the problem on the real frequency axis.But it badly resolves spectral functions at high energiesand cannot treat DMFT calculations with more than e.g.two bands. The limiting factor for this is the exponen-tial growth of the local
Hilbert space with the numberof bands. Only recently, a reformulation of the mappingproblem could avoid this exponential growth, but itis still unclear whether this can be efficiently exploitedin the context of DMFT. ED faces the problem of a lim-ited spectral resolution due to the limited number of bathsites it can treat, although recent publications could sub-stantially improve that. As the impurity problem of DMFT is one-dimensional,there has been a long-time interest to solve it using den-sity matrix renormalization group (DMRG), whichoperates on the class of matrix product states (MPS).DMRG features an unbiased energy resolution and showsno exponential growth of the local Hilbert space withrespect to the number of baths. It also works directlyon the real-frequency axis, avoiding analytic continua-tion. The earliest DMRG approach to spectral func-tions, the Lanczos algorithm approach, is computa-tionally cheap, but does not yield high-quality DMFT results due to its intrinsic numerical instability. Recentimprovements using a fully MPS-based representation ofthis algorithm are not sufficient to resolve this issue. The dynamical DMRG (DDMRG) approach yieldsvery precise results for single-site DMFT on the real fre-quency axis, but is computationally extremely costlyand therefore not competitive with other impurity solversfor DMFT.Recently, a new approach to spectral functions basedon expansions in Chebyshev polynomials representedwith matrix product states (CheMPS) was intro-duced by two of us in Ref. 28, which gave essen-tially the accuracy as the DDMRG approach at a frac-tion of the computational cost. At the same time,the availability of real-time evolution within time-dependent DMRG (tDMRG) and closely related meth-ods generally also permits access to spectral functions bya Fourier transformation. Both Chebyshev expansions(CheMPS) and tDMRG were recently seen to be ap-plicable to the solution of the DMFT. Both approachesare computationally cheaper than DDMRG and numer-ically stable. For the single-impurity single-band case,results on the real-frequency axis are excellent, but formore typical present-day DMFT setups involving clus-ters or multiple bands, results are not available in thecase of Chebyshev expansions or do so far not reach thequality of the competing QMC and NRG methods in thecase of real-time evolution.In this paper, we push the application of CheMPS toDMFT further: (i) We solve the dynamical cluster ap-proximation (DCA) for a two-site cluster and the DMFTfor a two-band Hubbard model. The accuracy of the re-sults for the latter case is better than those shown in Ref.35, where the problem has been solved using tDMRG.(ii) We consider the experimentally relevant case of finitedoping, which is significantly more complicated than thehalf-filled cases treated so far. (iii) We suggest a newtruncation scheme for CheMPS, which allows to main-tain the same error level at strongly reduced computa- a r X i v : . [ c ond - m a t . s t r- e l ] J u l tional cost. (iv) We establish that the Chebyshev recur-rence iteration can be interpreted as a discrete real-timeevolution. (v) By comparing different methods to setup CheMPS, we obtain another substantial increase incomputation speed. (vi) We discuss limitations of post-processing methods, which have been crucial to the suc-cess of DMRG as an DMFT impurity solver.With these improvements, CheMPS immediately pro-vides an efficient, precise and controlled way to solveDMFT problems with two baths (two-site clusters) on thereal-frequency axis with feasible extensions to problemswith more bands. The presentation proceeds as follows.After a general introduction to Chebyshev expansions ofspectral functions in Sec. II, we move on to discuss itsimplementation in the approximate framework of MPS:in Sec. III, we present a new truncation scheme, and inSec. IV, we discuss the mapping of the Hamiltonian tothe [ − , +1] convergence interval of Chebyshev polyno-mials, because this interacts non-trivially with efficientMPS calculations. Sec. V treats the post-processing ofChebyshev moments obtained in the expansion. Theseimprovements are then applied to various DMFT prob-lems. As the case of the single-impurity single-bandDMFT has been treated extensively in the literature andjust serves as an initial benchmark, we move those re-sults to the Appendix. In the main text, we give exam-ples for the relevance of our improvements to CheMPSby solving a two-site DCA in Sec. VI A and a single-sitetwo-orbital DMFT in Sec. VI B. Technical details of thesecalculations are again found in the Appendix. Sec. VIIconcludes the paper. II. CHEBYSHEV EXPANSION OF SPECTRALFUNCTIONS
In this Section, we establish notation and explain thegeneral ideas behind Chebyshev expansions of spectralfunctions. The zero-temperature single-particle Green’sfunction associated with a many-body hamiltonian H is G ( ω ) = (cid:104) E | c ω + i + − ( H − E ) c † | E (cid:105) , (1)where c † creates a particle in a particular quantum stateand | E (cid:105) is the ground state with energy E . The spectralfunction A ( ω ) = − π Im G ( ω ) reads A ( ω ) = (cid:104) E | c δ ( ω − ( H − E )) c † | E (cid:105) = (cid:88) n W n δ ( ω − ( E n − E )) , (2)with weights W n = |(cid:104) E n | c † | E (cid:105)| . If evaluated exactly ina finite system, A ( ω ) is a comb of delta peaks, which onlyin the thermodynamic limit becomes a smooth function A lim ( ω ). If evaluated in an approximate way that aver-ages over the finite-size structure of A ( ω ), it is possibleto extract A lim ( ω ) also from a sufficiently big finite-size system. Among various techniques that provide such anapproximation, the most popular one is the definitionof a broadened representation of A ( ω ) A η ( ω ) = (cid:88) n W n h η ( ω − E n ) (3)where the broadening function h η ( ω − E n ) is given by theGaussian kernel h η ( x ) = 1 √ πη e − x η . (4)Besides the Gaussian kernel, a Lorentzian kernel h η ( x ) = ηπ x + η (5)is often implicitly used as it emerges automatically whencomputing the spectral function A η = − π Im G ( ω + iη )from the shifted Green’s function G ( ω + iη ). In general, A η ( ω ) is indistinguishable from A lim ( ω ) if the latter hasno structure on a scale smaller than η .An efficient way to generate the broadened version A η ( ω ) of A ( ω ) is via iterative expansions in orthogonalpolynomials. Historically most frequently used in thiscontext is the Lanczos algorithm, which is intrinsicallynumerically unstable, though. By contrast, expansionsin Chebyshev polynomials can be generated in a numer-ically stable way. As they haven’t been used much in ei-ther the DMRG or DMFT community so far, we brieflyintroduce them based on Ref. 27. A. General implementation
The Chebyshev polynomials of the first kind T n ( x ) canbe represented explicitly by T n ( x ) = cos ( n arccos( x )) (6)or generated with the recursion T n ( x ) = 2 xT n − ( x ) − T n − ( x ) , T = 1 , T = x, (7)which is numerically stable if | x | ≤
1. Chebyshev polyno-mials are orthonormal with respect to the weighted scalarproduct (cid:90) − dx w n ( x ) T m ( x ) T n ( x ) = δ nm , (8a) w n ( x ) = 2 − δ n π √ − x . (8b)Any sufficiently well-behaved function f ( x ) | x ∈ [ − , canbe expanded in Chebyshev polynomials f ( x ) = ∞ (cid:88) n =0 w n ( x ) µ n T n ( x ) , (9a) µ n = (cid:90) − dxf ( x ) T n ( x ) , (9b)where the definition of the so-called Chebyshev moments µ n via the non-weighted scalar product follows when ap-plying (cid:82) − dx T m ( x ) . . . to both sides of (9a).If f ( n ) is smooth, the envelope of µ n decreases atleast exponentially to zero with respect to n ; if f ( n ) isthe step function, the envelope decreases algebraically;and if f ( n ) is the delta function, the envelope remainsconstant. For a smooth function, the truncated ex-pansion f N ( x ) = (cid:80) Nn =0 w n ( x ) µ n T n ( x ) therefore approx-imates f ( x ) very well if N is chosen high enough. Butfor the delta function, any truncated expansion yieldsan approximation with spurious (Gibbs) oscillations. Acontrolled damping scheme for the oscillations, the so-called kernel polynomial approximation (KPM), can beobtained with a simple modification of the Chebyshevexpansion, f kernel N = N (cid:88) n =0 w n ( x ) g n µ n T n ( x ) , (10a) g n = ( N − n + 1) cos πnN +1 + sin πnN +1 cot πN +1 N + 1 , (10b)where g n is the so-called Jackson kernel that leads toa very good Gaussian approximation h η ( x ) ( x ) with x -dependent width η ( x ) = √ − x π/N of the delta func-tion, and hence directly leads to (4).In the case of the spectral function (2), one aims atan expansion of a superposition of delta functions. Thiscan in practice often be done without damping: Whenexpanding (2) in Chebyshev polynomials, the integrationin (9b) averages over the delta-peak as well as over thefinite-size peak structure of A ( ω ). If the weights W n varyslowly on the scale of the spacing of finite-size peaks, thesequence µ n approaches zero as soon as the characteristicform of this slow variation is resolved. The value of n at which this pseudo -convergence occurs is the one thatresolves the spectral function in the thermodynamic limit A lim ( ω ), provided that A lim ( ω ) has no structure on asmaller scale than the spacing of finite-size peaks. Onlyfor much higher values of n , the Chebyshev momentsstart deviating from zero again to then oscillate forever,resolving first the finite-size structure of A ( ω ) and finallythe delta-peak structure. Therefore, if one can generatethe sequence up to pseudo -convergence, then there is noneed for Jackson damping. B. Operator valued Chebyshev expansion
In order to expand the spectral function (2), one usu-ally introduces a rescaled and shifted version of H in or-der to map its spectrum into the interval [ − , H (cid:48) = H − E + ba , ω (cid:48) = ω + ba . (11) Obviously, there is a lot of leeway in the choice of a and b , which will be found to have large implicationsfor CheMPS (Sec. IV). Generally, A ( ω ) = 1 a A (cid:48) (cid:0) ω + ba (cid:1) , where A (cid:48) ( ω (cid:48) ) = (cid:104) t | δ ( ω (cid:48) − H (cid:48) ) | t (cid:105) , | t (cid:105) = c † | E (cid:105) . (12)Expanding A (cid:48) ( ω (cid:48) ) in Chebyshev polynomials yields themoments µ n = (cid:90) − dω (cid:48) (cid:104) t | δ ( ω (cid:48) − H (cid:48) ) | t (cid:105) T n ( ω (cid:48) )= (cid:88) i (cid:90) − dω (cid:48) (cid:104) t | δ ( ω (cid:48) − E (cid:48) i ) T n ( ω (cid:48) ) | E i (cid:105)(cid:104) E i | t (cid:105) = (cid:104) t | t n (cid:105) , | t n (cid:105) = T n ( H (cid:48) ) | t (cid:105) , (13)Inserting the recursive definition (7) of T n ( H (cid:48) ) in the def-inition of | t n (cid:105) one obtains a practical calculation schemefor the power series expansion of T n ( H (cid:48) ) | t n (cid:105) = 2 H (cid:48) | t n − (cid:105) − | t n − (cid:105) , (14a) | t (cid:105) = c † | E (cid:105) , | t (cid:105) = H (cid:48) | t (cid:105) . (14b)One can double the expansion order with the followingrelation µ n − = 2 (cid:104) t n | t n − (cid:105) − µ , (15a) µ n = 2 (cid:104) t n | t n (cid:105) − µ , (15b)but has to be aware of the fact that moments computedthis way are more prone to numerical errors. C. Retarded fermionic Green’s function
In the case of fermionic problems, as encountered inDMFT, an additional technical complication comes up.The spectral representation of the fermionic retardedGreen’s function is the sum of its particle and hole parts A ( ω ) = A > ( ω ) + A < ( − ω ) ,A > ( ω ) = (cid:104) E | c δ ( ω − ( H − E )) c † | E (cid:105) ,A < ( ω ) = (cid:104) E | c † δ ( ω − ( H − E )) c | E (cid:105) . (16)As A ≶ ( ω ) have steps at ω = 0, their representation interms of smooth polynomials is notoriously ill-conditio-ned. One should therefore try to represent the smoothfunction A ( ω ) by a single Chebyshev expansion: Allow-ing for two different rescaling prescriptions, one has A > ( ω ) = 1 a (cid:88) n w n ( ω (cid:48) ( ω )) µ >n T n ( ω (cid:48) ( ω )) (17a) A < ( − ω ) = 1 a (cid:88) n w n ( ω (cid:48) ( − ω )) µ So far everything has been general, or it was somehowassumed that all calculations can be carried out exactly,which meets severe limitations in computational practice.Representing Chebyshev states | t n (cid:105) with matrix prod-uct states (MPS) enables more efficient computationsthan in an exact representation, as the size of the effec-tive Hilbert space can be tremendously reduced. As anMPS is usually only an approximate representation of astrongly correlated quantum state, the issue of optimalcompression, i.e. the representation of a quantum stateas an MPS using finite-dimensional matrices with a min-imal loss of accuracy (information), is crucial. Here, weargue in the following that instead of controlling the max-imal matrix dimension, one should rather controlthe cumulated truncated weight (a proxy measure of theloss of accuracy), allowing for more efficient and morecontrolled calculations of Chebyshev moments. A. Adaptive matrix dimension If one follows through the recursive scheme for Cheby-shev vectors, one starts out from a ground state, whichwe may assume has been obtained by a standard DMRG(MPS) calculation to extremely high precision, thismeans that an optimally compressed starting MPS isavailable where matrices have some computationally fea-sible dimension at very small loss of accuracy comparedto the exact starting state. This, in turn, yields an ex-tremely precise starting Chebyshev state | t (cid:105) . Now, ineach step of the recursion (14a), one applies H (cid:48) and sub-tracts a preceding Chebyshev state. As is well-knownfor MPS, the application of H (cid:48) (and to a lesser extentthe subtraction) lead to a drastic increase in matrix di-mension, which necessitates a state compression (Sec. 4.5of Ref. 17) of the new Chebyshev state | (cid:101) t n (cid:105) to a com-putationally manageable state | t n (cid:105) with smaller matrix dimension m , which generates the error δµ n = (cid:104) t | (cid:101) t n (cid:105) = (cid:104) t | t n (cid:105) ± δ, (19) | (cid:101) t n (cid:105) = 2 H (cid:48) | t n − (cid:105) − | t n − (cid:105) ,δ = (cid:12)(cid:12) (cid:104) t | ( | (cid:101) t n (cid:105) − | t n (cid:105) ) (cid:12)(cid:12) < || t (cid:105)| (cid:12)(cid:12) | (cid:101) t n (cid:105) − | t n (cid:105) (cid:12)(cid:12) < || t (cid:105)| ε compr ( m ) . Here, we used the upper error bound provided by thecumulated truncated weight ε compr ( m ) (cid:12)(cid:12) | (cid:101) t n (cid:105) − | t n (cid:105) (cid:12)(cid:12) ≤ ε compr ( m ) = L − (cid:88) i =1 (cid:15) i ( m ) , (20)where (cid:15) i ( m ) is the sum over the discarded reduceddensity-matrix eigenvalues per bond and the sum over i is over all bonds. This error bound for a single step ofthe recursion unfortunately does not provide a statementabout the total error that accumulates over all compres-sion steps in preceding Chebyshev recursion steps. Still,we experienced that the numerical stability of the Cheby-shev recursion rather leads to a helpful compensation oferrors of single recursion steps. Fig. 1 shows that the total error stays at the order of the error of a single step || t (cid:105)| ε ( m ) also for high iteration numbers n . In the casein which one fixes the matrix dimension m , Fig. 1 showsa steady, uncontrolled increase of the total error. Thisis particularly undesirable in view of the desired post-processing of Chebyshev moments (Sec. V).Another possibility would be to fix the local discardedweight (cid:15) i ( m ) as defined in (20). But this does in general not lead to a viable computation scheme for impuritymodels: In the simplest and most-employed chain repre-sentation of impurity models, the impurity site is locatedat an edge of the chain. Fixing the same value for (cid:15) i ( m )for all bonds then leads to extremely high matrix dimen-sions in the center of the chain, i.e. in the center of thebath, where entanglement for systems with open bound-ary conditions is maximal. The relevant entanglement,by contrast, is the one between the impurity site andthe bath. This becomes clear when noticing that uponprojecting the Chebyshev state | t n (cid:105) on | t (cid:105) to compute µ n , only correlations with respect to the local excitation c † | E (cid:105) are measured. The high computational effort ofhigh matrix dimensions that follows when faithfully rep-resenting entanglement within the bath, is therefore invain. For geometries with the impurity at the center, likethe two-chain geometry used for the two-bath problemsin this paper, the preceding argument is not valid. An in-homogeneous distribution of matrix dimensions with highvalues at the center and low values at the boundaries is a priori consistent with open boundary conditions. Thisdistribution can therefore be achieved by fixing a con-stant value for (cid:15) i ( m ) for each bond. Another possibletruncation scheme could be obtained by using an esti-mator for the correlations of the impurity with the bath,which then fixes the matrix dimensions as a function ofbonds m ( i ) (distance to the impurity). Both approaches > n × - compr =10 -3 m=70 FIG. 1. (Color online) Error of Chebyshev moments µ >n (asthey appear in (17a)), computed as ∆µ >n = (cid:12)(cid:12) µ >n − ˜ µ >n (cid:12)(cid:12) , where˜ µ >n is obtained with a quasi-exact calculation with high matrixdimension m = 200. If one fixes the matrix dimension m , theerror steadily increases. If, instead, one fixes the cumulatedtruncated weight ε compr , the error remains approximately con-stant and does not accumulate. This is the procedure followedin this paper. As here, || t (cid:105)| = 1, ε compr equals the upper er-ror bound of a single compression step. Results shown are forthe spectral function of the half-filled single-impurity Ander-son model (SIAM) (Appendix C 1) with semi-elliptic densityof states of half-bandwidth D , interaction U = 2 D , repre-sented on a chain with L = 40 lattice sites. This is equivalentto considering the local density of states at the first site of afermionic chain with constant hopping t = D/ U = 4 t that acts solely at the first site. constitute possible future refinements. For simplicity, inthis paper, we consider the truncation scheme that fixesa constant value of m based on the cumulative truncatedweight. B. State compression During the repeated solution of (14a) we monitor thetruncated weight ε compr . If ε compr exceeds a certainthreshold of the order of 10 − to 10 − , we slightly in-crease the matrix dimension m , and repeat the compres-sion. For the first compression step we take as an initialguess the previous Chebyshev state | t n − (cid:105) . For repeatedcompression steps we take as an initial guess the stateof the previous compression step. It turns out that inpractice one almost never faces repeated compressions,which gains one approximately a factor 2 in computationspeed compared to the error monitoring of Ref. 28: inRef. 28, the authors keep the matrix dimension fixed andvariationally compress an exact representation of theright hand side of (14a) for fixed m by repeated itera- tions (“sweeps”) until the error (cid:12)(cid:12)(cid:12)(cid:12) − (cid:104) t (cid:48) n | t n (cid:105)|| | t (cid:48) n (cid:105) || || | t n (cid:105) || (cid:12)(cid:12)(cid:12)(cid:12) , (21)drops below a certain threshold. Here, | t (cid:48) n (cid:105) denotes thestate before a sweep, and | t n (cid:105) the state after a sweep.This error measure is not related to the factual error ofChebyshev moments, for any but the first sweep. Itsmonitoring is costly to compute and leads to at least twocompression sweeps. IV. OPTIMAL CHEBYSHEV SETUP One can generally state that the effectiveness of theMPS evaluation of the Chebyshev recursion (14a) fora certain system is unknown a priori but must be ex-perienced by observing how strong entanglement in theChebyshev vectors, and therefore matrix dimension m needed for a faithful representation grows as comparedto the speed of convergence of µ n . For very high iter-ation numbers one will always reach a regime in whichmatrix dimensions have grown so much that further cal-culations become too expensive computationally. This isknown from tDMRG as hitting an exponential wall anddefines an accessible time scale , or in our case, an acces-sible expansion order. In the case of the computation ofChebyshev moments, the accessible time scale stronglydepends on the choice of the shifting parameter b , whichleads us to consider the two cases b = 0 and b = − a .Comparing these cases, one finds a much slower speedof convergence of the Chebyshev moments in the case b = 0 than in the case b = − a . Putting that differently:per fixed amount of entanglement growth (application of H in one step of (14a)), much less information aboutthe spectral function is extracted in case b = 0 than incase b = − a . Independent of that, one finds that theadvantage of the choice b = 0 to provide one with an an-alytic expression for A ( ω ) in terms of a single Chebyshevexpansion (Sec. II C) can be detrimental. We thereforeneed to study both cases in more detail. A. No shift: b = 0 If choosing b = 0, one can derive a scaling property ofChebyshev moments that simplifies extracting the ther-modynamic limit as well as the examination of computa-tional performance.The spectral function of a one-particle operator A ( ω )is non-zero only in the vicinity of the groundstate en-ergy ω = 0, up to a distance of the order of the single-particle bandwidth W single . The rescaled spectral func-tion A (cid:48) ( ω (cid:48) ) is non-zero up to a distance of W single /a from ω (cid:48) = 0. For all rescaling parameters a that have beenproposed up to now, one has W single /a < . Usu-ally W single /a is much smaller than the upper bound .As arccos( x ) = π/ − x − x / . . . is well approxi-mated by its linear term already for | x | < . 5, Chebyshevpolynomials (6) behave like a shifted cosine function inthe region where A (cid:48) ( ω (cid:48) ) is non-zero. The expansion of A (cid:48) ( ω (cid:48) ) in Chebyshev polynomials is therefore essentiallyequivalent to a Fourier expansion. This means that theiteration number n of the Chebyshev expansion has thesame meaning as a discrete propagation time, the evo-lution of which is mediated by simple applications of H instead of the ordinary continuous time propagation e − iHt . To answer the question of whether an ordinarytime evolution is more effective in generating informa-tion about the spectral function, one has to study the en-tanglement entropy production of repeated applicationsof H compared to the one of e − iHt . The following resultsare first steps in this direction.In discrete time evolution, the rescaling of the fre-quency directly translates to an inverse scaling of time.Considering two calculations of Chebyshev moments, onefor µ (1) n performed with H (cid:48) and another for µ ( a ) n per-formed with H (cid:48) /a , one therefore has the simple approxi-mate relation µ (1) n ∼ (cid:104) t | cos( nH (cid:48) ) | t (cid:105) = (cid:104) t | cos( anH (cid:48) /a ) | t (cid:105) ∼ µ ( a ) na . (22)This means that if rescaling with a , one has to compute a times more Chebyshev moments than in the case withoutrescaling. An exact version of statement (22) is givenin (A2) in Appendix A. Fig. 2(a) illustrates the scalingproperty (22) for a system of fixed size. a. Extracting the thermodynamic limit. One directapplication of the scaling property (22), lies in the studyof the thermodynamic limit by comparing systems of in-creasing size L . For low values of n , even small systemshave the same Chebyshev moments as in the thermody-namic limit. Finite-size features are averaged out in theintegral (9b) as long as T n ( x ) oscillates slowly enough. T n ( x ) oscillates n times on [ − , N th order Cheby-shev expansion therefore resolves features on the scale2 /N , which on the original energy scale is 2 a/N . Finite-size oscillations appear at a spacing of W single /L , where W single is the single-particle bandwidth. Equating reso-lution with the spacing of finite-size oscillations2 a/N finsize = W single /L, (23)gives the expansion order N finsize at which finite-size fea-tures are first resolved. Fig. 2(b) illustrates these state-ments by comparing Chebyshev moments computed fordifferent system sizes. b. Optimizing computation time. Fig. 3 shows howcomputation time depends on the rescaling constant a for the example of the moments shown in Fig. 2(a). Asalready qualitatively stated previously , one observesthat upon using a lower value of a computation time isreduced. In all cases, computation time diverges expo-nentially (Fig. 2(b)). Note that rescaling with a higher + n (a) a =15Da =20Da =40Da =80D W + n (b) L =10L =20L =40 FIG. 2. (Color online) Panel (a): Chebyshev moments µ >n vs n/a for fixed system size and different values of a and b = 0. Except for a different total number of points, therescaled moments all lie on the line obtained when a → ∞ and n/a becomes continuous. Here, we study the half-filledSIAM (Appendix C 1) with semi-elliptic density of states ofhalf-bandwidth D and U = 2 D , represented on a chain withlength L = 80. The full many-body bandwidth is W (cid:39) D .Panel (b): Chebyshev moments for different system sizes L .Except for the system size and the scaling parameter, param-eters are the same as in panel (a). Here all calculations weredone with a rescaling constant of a = 20 D . For low values of n , the results for different system sizes are virtually indistin-guishable. For higher values of n , moments start to disagreeas finite-size features start to be resolved. The L = 80 andthe L = 40 results would be indistinguishable in this plot. value of a allows to compute at smaller matrix dimen-sions. Note further that if choosing a too small, numer-ical errors can render the recursion (14a) unstable. Incontrast to common belief, it is possible to use muchsmaller values of a than the full many-body bandwidth.Achieving even smaller values of a can be done with theso-called energy truncation , but after several tests, wedid not find this to lead to an effective speed-up of cal-culations. We therefore discard it in our calculations as m a t r i x d i m e n s i o n m (a) a =80D Wa =40Da =20Da =15D c p u t i m e ( m i n ) (b) a =80D Wa =40Da =20Da =15D FIG. 3. (Color online) Performance of the adaptive matrixdimension algorithm (Sec. III A) for the example described inthe caption of Fig. 2. Panel (a): Adaption of matrix dimen-sions for different rescaling factors, fixing a truncation errorof ε compr = 10 − . Panel (b): Computer time needed to gen-erate the same amount of information for different scalingsrunning on a single-core 2.0 GHz workstation. Solid lines:fixing a truncation error of ε compr = 10 − . Dashed lines: ε compr = 5 × − . The iteration number where the irregularbehavior of the dashed line for a = 15 D starts correspondsto the point where numerical errors render the Chebyshev re-cursion unstable. Note that while small a leads to the largestmatrix sizes, which is costly in MPS, the overall cost of CPUtime nevertheless is lowest, as a smaller expansion order isneeded. a source of additional tuning parameters. We have alsotested the idea of Ganahl et al. to map the spectrumof H into [ − , 1] via 1 − exp( βH ). The idea might beworth to study in more detail, but again, we could notgain any performance improvement over a simple rescal-ing procedure. B. Shifting by b = − a . The choice b = − a in (11) makes an analytic expres-sion of the complete spectral function A ( ω ) = A + ( ω ) + A − ( − ω ) in terms of a single Chebyshev expansion impos-sible, but has beneficial effects on the computation time.This is to be understood in the following sense: Due tothe increased oscillation frequency of T n ( x ) close to theinterval boundaries of [ − , − , 1] (see the discussion below(10b)). It is therefore desirable to shift the relevant partof the spectral function, the part slightly above the Fermiedge, to match the left boundary − 1. This is achieved bythe choice b = − a . In practice, one adds a small correc-tion a(cid:15) , (cid:15) ∼ − , to avoid problems with the divergingweight function w n ( x ) in (8b).Another advantage of the b = − a setup is that onecan use a smaller scaling constant a than in the b = 0setup. The Chebyshev iteration becomes unstable whenthe iteration number n becomes so high that | t n (cid:105) hasaccumulated erroneous contributions from eigen stateswith eigen energies E (cid:48) n = ( E n − E + b ) /a > 1. Forfixed a , the additional subtraction in the b = − a setupensures that the instability appears for a higher iterationnumber than in the b = 0 setup. Therefore, the b = − a setup allows smaller values of a . We finally note that thechoice b = − a is equivalent to the choice suggested byWeiße et al. , if one rescales with the full many-bodybandwidth a = W . In this case, the computation canbe carried out to arbitrarily high order and will neverbecome unstable. In the b = 0 setup, one would haveto choose a = 2 W to reach arbitrarily high expansionorders.In Fig. 4(a), we plot Chebyshev moments for bothtypes of shifts b = 0 and b = − a . The moments ob-tained for b = 0 show a slow structureless oscillationwhereas the moments obtained for b = − a show a muchfaster oscillation. Fig. 4(b) shows that upon using thesame rescaling constant a and the same expansion or-der N = 100, which leads to very similar entanglementgrowth, both shift types differ strongly in the achievedresolution. To resolve at least the right Hubbard peakwith a b = 0 calculation at the resolution of b = − a cal-culation, one needs N = 300 moments. As computationtime increases exponentially (Fig. 3(b)) with respect toexpansion order N in both cases, this difference is highlyrelevant.We apply both setups, b = 0 and b = − a , to the bench-mark test of the DCA in Sec. VI A, and find a signif-icant speed-up for b = − a at a small loss in accuracy.Previously, only b = 0 has been considered for the so-lution of the DMFT. > n (a) b =(-1 + )ab =0 A > () D (b) b =(-1 + )a N=100b =0 N=100b =0 N=300 FIG. 4. (Color online) Local particle density of states of thehalf-filled SIAM (Appendix C 1) with semi-elliptic density ofstates of half-bandwidth D . L = 40, U = 2 D and a = 30 D in all cases. Panel (a): Chebyshev moments. Lines connectevery 4th moment and by that reveal the relevant slow oscilla-tion. They are a guide to the eye. Panel (b): Correspondingspectral functions evaluated using Jackson damping (10b).The b = 0 calculation requires three times more iterationsthan the b = − a calculation to resolve the right Hubbardpeak with the same resolution. In this case the central peakis still much better resolved for b = − a . V. POST-PROCESSING MOMENTS Whereas Jackson damping (10b) can be seen as onepossibility to post-process Chebyshev moments in orderto achieve uniform convergence even for the truncatedChebyshev expansion of a delta function, there is an-other, fundamentally different approach.The computation of the Chebyshev moments becomesvery costly for high iteration numbers. In the case inwhich Chebyshev moments start to follow a regular pat-tern when n exceeds a certain threshold, it is possible tocontinue this pattern to infinity, and one can avoid thecostly computation of moments. Consider the typical example in which the spectral function is a superposi-tion of Lorentzians (quasiparticle peaks) and of a slowlyvarying background density. As for low values of n , T n ( x )extracts information via (9b) only about the slowly vary-ing background density, while for high values of n , T n ( x )extracts information only about the sharp and regularLorentzian structures, µ n starts to follow a regular pat-tern for high numbers of n . For a sum of Lorentzians,with weights α i , widths η i , and positions ω i , this patterncan be obtained analytically: A Lor ( ω ) = (cid:88) i α i η i π ω − ω i ) + η i , ⇒ µ n (cid:39) (cid:88) i α i cos( n ( ω i − π e − nη i , (24)as shown in Appendix B. If one recalls (Sec. IV) thatthe Chebyshev recursion corresponds to a discrete timeevolution if choosing b = 0, the result of (24) could havebeen anticipated.Fig. 5(a) shows the spectral density for a SIAM to-gether with a fitted superposition of three Lorentzians.Their difference corresponds to a background densitythat is composed of either slowly varying features or fea-tures with negligible weight. Fig. 5(b) shows the corre-sponding Chebyshev moments. The slowly varying back-ground density only contributes for the first 200 mo-ments. After that, the Chebyshev moments for the su-perposition of Lorentzians starts to be a very good ap-proximation to the original moments, and it seems un-necessary to compute more than about 400 moments.For 200 < n < linear prediction . The linear problem can beanalytically reformulated as a matrix inversion problem.Its solution is faster and more stable than that of theoriginal non-linear problem. This allows in principle tooptimize a superposition of many more Lorentzians thanin the non-linear case. 1. Linear prediction In the context of time evolution linear prediction hasbeen long established in the DMRG community, butit has only recently been applied to the computation ofChebyshev moments. The optimization problem for thesequence µ n becomes linear, if the sequence can be de-fined recursively ˜ µ n = − p (cid:88) i =1 a i µ n − i , (25) A () W (a) A SIAM A Lor A SIAM -A Lor n (b) FIG. 5. (Color online) Panel (a): A SIAM ( ω ) for a semi-ellipticdensity of states, half filling and U = 2 D (Appendix C 1).Quantities are shown in units of the full many-body band-width W . The superposition of three Lorentz peaks A Lor ( ω )has been fitted to A SIAM ( ω ). Panel (b): CorrespondingChebyshev moments. The result presented here was obtainedwith a L = 40 fermionic chain and CheMPS. It agrees withthe result of Raas et al. , see Appendix C 1. The legend inpanel (a) is valid also for panel (b). which is easily found to be equivalent to (24) . Thestrategy is then as follows. Compute n = N c Chebyshevmoments, and predict moments for higher values of n using (25). The coefficients a i are optimized by minimiz-ing the least-square error (cid:80) n ∈N fit | ˜ µ n − µ n | for a subset N fit = { N c − n fit , . . . , N c − , N c } of the computed data.We confirmed n fit = N c / small enough to go beyond spurious short-time behaviorand large enough to have a good statistics for the fit.Minimization yields R a = − r , a = − R − r , (26) R ji = (cid:88) n ∈N fit µ ∗ n − j µ n − i , r j = (cid:88) n ∈N fit µ ∗ n − j µ n . We found that linear prediction loses its favorable fil-ter properties if choosing p to be very high. There-fore one should restrict the number of Lorentzians to p = min( n fit / , δ = 10 − to the diagonal of R in order to enable the inversion of the singular matrix R . Defining M = − a − a − a . . . − a p . . . 00 1 0 . . . . . . , one obtains the predicted moments ˜ µ N c + n = ( M n µ N c ),where µ N c = ( µ N c − µ N c − . . . µ N c − p ) T . The matrix M usually has eigenvalues with absolute value larger than1, either due to numerical inaccuracies or due to the factthat linear prediction cannot be applied as µ n rather in-creases than decreases on the training subset N fit . In or-der to obtain a convergent prediction, we set the weightsthat correspond to these eigenvalues to zero measuringthe ratio of the associated discarded weight compared tothe total weight. If this ratio is higher than a few per-cent, we conclude that linear prediction cannot yet beapplied and restart the Chebyshev calculation to increasethe number of computed moments N c . 2. Failure of linear prediction It is not a priori clear that the spectral function canbe well approximated by a superposition of Lorentzians,although this is true for the SIAM as shown in Fig. 5.Other types of smooth functions lead to a different func-tional dependence of the moments on n than the expo-nentially damped behavior. Close to phase transitions,e.g. one might find an algebraic decay in the time evolu-tion, corresponding to an algebraic decay in the Cheby-shev moments. If the spectral function has rather Gaus-sian shaped peaks, the decrease of Chebyshev momentsis ∝ e − ( σn ) (Appendix B). For both scenarios, linearprediction is a non-controlled extrapolation scheme. Itstill extracts oscillation frequencies (peak positions) withhigh reliability, but predicts a wrong decrease of the en-velope, which often leads to an overestimation of peakweights.In practice it turns out that a combination of damping with a Jackson kernel (Kernel Polynomial Method) and linear prediction is a powerful way to get controlled esti-mates for the spectral function. While damping alwaysunderestimates peak heights, linear prediction typicallyoverestimates peak heights. Both methods trivially con-verge to the exact result, when N c → ∞ . One thereforeobtains upper and lower bounds for the spectral function.This is particularly valuable in the DMFT as overesti-mated (diverging) peak heights can spoil convergence ofthe DMFT loop.A historically much used alternative to linear predic-tion, suitable for arbitrary forms of the spectral function,is an extrapolation of Chebyshev moments using maxi-mum entropy methods . These suffer from severe nu-merical instabilities, though. Of course, one might also0think of fitting another ansatz than the one of the expo-nential decrease. As it is a priori not clear which ansatzshould be better, it is meaningful to stick to the easily im-plemented linear prediction that is moreover known to beapplicable for the description of quasi-particle features. VI. RESULTS FOR DMFT CALCULATIONSWITH TWO BATHSA. Results for two-site DCA (VBDMFT) In order to benchmark the Chebyshev technique for atwo-bath situation, which goes beyond previous work (see Appendix C), we study the Hubbard model on thetwo dimensional square lattice H Hub = (cid:88) k σ ε k c † k ,σ c k ,σ + U (cid:88) i n i ↑ n i ↓ , (27) ε k = − t (cos( k x ) + cos( k y )) − t (cid:48) cos( k x ) cos( k y ) . in a two-site dynamical cluster approximation (DCA)developed by Ferrero et al. . This so-called valencebond DMFT (VBDMFT) is a minimal description of thenormal phase of the high-temperature superconductors,using a minimal two patches DCA cluster. It leads to asimple physical picture of the pseudogap phase in termsof a selective Mott transition in the momentum space.We choose this model here as a benchmark since its solu-tion contains low energy features in the spectral functions(pseudogap), which have required high-precision QMCcomputations followed by a careful Pad´e analytic con-tinuation. Moreover, real-frequency computations arevery important for the comparison with experiments thatmeasure e.g. the optical conductivity along c-axis. Itis therefore a non-trivial case where DMRG impuritysolvers would bring significant improvements over theQMC in practice.To set up the VBDMFT, one splits the Brioullin zoneinto a central patch P + = { k (cid:12)(cid:12) | k x | < k ∧ | k y | < k } ,where k = π (1 − / √ border patch P − = { k (cid:12)(cid:12) k / ∈ P + } . In the DCA, the k -dependence of theself-energy Σ κ ( ω ) within each patch is neglected and onecomputes a Green’s function for a patch by averagingover all k vectors in the patch G κ ( ω ) = 1 | P κ | (cid:88) k ∈ P κ ω + µ − ε k − Σ κ ( ω ) , (28a) Σ κ ( ω ) = G κ ( ω ) − − G κ ( ω ) − . (28b)Representing the non-interacting baths in a chain-geometry, and taking the two impurities to be the firstof two chains c κσ ≡ c κσ , the model Hamiltonian that needs to be solved is H = H d + H b, + + H b, − H d = (cid:88) κ = ± σ = ↑ , ↓ (cid:0) t κ + ε (cid:1) n κσ + U (cid:88) κ = ± κ = − κ (cid:0) n κ ↑ n κ ↓ + n κ ↑ n κ ↓ + c † κ ↑ c † κ ↓ c κ ↓ c κ ↑ + c † κ ↑ c † κ ↓ c κ ↓ c κ ↑ (cid:1) , (29) H b,κ = L κ (cid:88) i =0 ,σ t iκ ( c † iκσ c i +1 ,κσ + h.c.) + L κ (cid:88) i =1 ,σ ε iκ n iκσ , where ε = − µ and the term t κ = | P κ | (cid:80) k ∈ P κ ε k ac-counts for high-frequency contributions of the hybridiza-tion function (see Appendix D 4).The κ -space interaction term in (29) arises when diag-onalizing the hybridization function of a real-space two-site cluster c ± σ = √ ( c σ ± c σ ), where c σ , c σ are an-nihilation operators for the cluster sites in real-space,and c ± σ for the cluster sites in κ space. In real-space,the interaction is a simple Hubbard expression, but thenthe hybridization function is non-diagonal. A diagonalhybridization function, which leads to two uncoupledbaths for the patches and by that allows a simple chain-geometry for the whole system, is therefore only possiblein κ -space. The more complex form of the interaction in κ -space does not affect the efficiency of DMRG.We iteratively solve the self-consistency equation ob-tained by inserting the self-energy estimates of the impu-rity model (29) into the lattice Green functions (28a). Wedo that on the real-energy axis with an unbiased energyresolution. The details of this calculation are describedin Appendix D.In Fig. 6(a) and (b), we compare our CheMPS resultsfor the spectral densities of the two momentum patcheswith those of Ferrero et al. obtained using CTQMCand analytical continuation. We observe a good over-all agreement between the two methods, in particularat low frequencies. Low energy features (pseudogap), inparticular in A − ( ω ), are well reproduced by both meth-ods. At high energy (Hubbard bands) however, thereare some differences between QMC and CheMPS (andalso between the two variants of CheMPS). This is tobe expected since the Pad´e analytic continuation tech-nique used on the QMC data in Ref. 44 is not a precisionmethod at high energy.In Fig. 6(c) and (d), we do the analogous compari-son on the imaginary axis, and find much better agree-ment. On the imaginary axis, the QMC results can beconsidered numerically exact. The very low tempera-ture ( βD = 200) used for QMC should yield results thatare indistinguishable from a zero-temperature calcula-tion. The slight disagreement of our data and the QMCdata on the Matsubara axis could probably be removedif we were able to reach higher expansion orders. OneDMFT iteration for the presented b = 0 calculation tookaround 5 h running on four cores with 2 . A + () = - I m G + () (a)2 1 0 1 2 3 4/D0.00.20.40.60.81.0 A - () = - I m G - () (b) this workFerrero (2009) G + ( i n ) (c)0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4i n /D0.80.60.40.20.0 G - ( i n ) (d) ReG(i n )ImG(i n ) FIG. 6. (Color online) Spectral functions (a,b) and Green’sfunctions on the imaginary axis (c,d) within VBDMFT for U = 2 . D and n = 0 . 96. We compare our zero-temperatureCheMPS results (solid lines) with CTQMC data for T =1 / 200 (dashed lines) from Ferrero et al. . For this compu-tation, we used the b = 0 setup, a chain length of L = 30 perbath, a truncation error of ε compr = 10 − , and N/a = 60 /D , a = 40 D . A + () = - I m G + () (a)2 1 0 1 2 3 4/D0.00.20.40.60.81.0 A - () = - I m G - () (b) this workFerrero (2009) G + ( i n ) (c)0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4i n /D0.80.60.40.20.0 G - ( i n ) (d) ReG(i n )ImG(i n ) FIG. 7. (Color online) The same comparison as in Fig. 6. Forthis computation, we used the b = − a setup, a chain lengthof L = 40 per bath, a truncation error of ε compr = 10 − , and N = 450 and a = 15 D . For the b = − a setup, one can usea smaller value of a as in the b = 0 setup, as discussed inSec. IV. L = 30 lattice siteseach. We did not observe changes for higher bath sizesup to L = 40, but could not reach high enough expan-sion orders for chains longer than L = 40. We computed N = 2500 moments using a scaling constant a = 40 D ,which corresponds to the full bandwidth.The calculation can be accelerated significantly by us-ing the b = − a setup of Sec. IV B and avoiding linearprediction . This leads to the same quality of agreementwith QMC on the Matsubara axis, but on the real axis,peaks are a bit less pronounced while the pseudogap isstill well resolved (Fig. 7). While the study of systemswith higher bath sizes increases the computational costtremendously in the b = 0 setup, we could easily go to L = 50 within the b = − a setup. This did not changethe results. Computation times varied from 1 . L = 30, over 3 h for L = 40 to around 10 h forthe L = 50 calculation. We computed N = 450 momentsusing a scaling of a = 15 D in all cases. B. Single-site two-orbital DMFT In the following, we apply CheMPS to the DMFTtreatment of the two-orbital Hubbard model H = (cid:88) k νσ ε k ν n k νσ + U (cid:88) iν n iν ↑ n iν ↓ + (cid:88) iσσ (cid:48) ( U − δ σσ (cid:48) J ) n i σ n i σ (cid:48) + J (cid:88) iνσ c † iνσ ( c † iν σ c iνσ + c † iνσ c iν σ ) c † iνσ (30)on the Bethe lattice. We study a parameter regime closeto the Metal-Insulator phase transition. This regime iscomputationally particularly expensive and we had to usea logarithmic discretization to reach Chebyshev expan-sion orders at which spectral functions are completelyconverged with respect to expansion order and systemsize. The linear discretization was feasible in the case ofthe VBDMFT studied in the previous section, as there,we faced a smaller entanglement entropy production dur-ing Chebyshev iterations.Using a logarithmic discretization is not necessary forCheMPS. But as it leads to exponentially decaying hop-ping constants, it gives rise to three advantages: (i) Onecan use smaller scaling constants a as the many-bodybandwidth is considerably reduced due to the exponen-tially small value of most hopping constants in the sys-tem. (ii) One faces a smaller entanglement entropy pro-duction: at the edges of the bath chains (far away fromthe impurity), hopping constants are exponentially small,and application of H therefore creates much less entan-glement than in the case in which a linear discretizationis used. In (14a), the action of H (cid:48) on | t n − (cid:105) is thenonly a small perturbation for most parts of the system,and the recursion is therefore dominated by the second A () D (a) this workNRGQMC A () D (b) this workNRGQMC FIG. 8. (Color online) Spectral function for the two-bandHubbard model. Panel (a): U/D = 1 . n = 2 (half fill-ing). Panel (b) U/D = 3 . n = 1 (quarter filling). In bothcases J = U , U (cid:48) = U − J . We fixed a truncation error ε compr = 10 − , used a scaling a = 25 D , computed N c = 150moments and used linear prediction. To represent the twobaths, we used two chains of length L = 20 each, obtainedwith a logarithmic discretization parameter of Λ = 2, leadingto grid energies Λ − n (see e.g. Ref. 10). The NRG calculationwas done for temperature T /D = 0 . T /D = 0 . 01. Both should be almost indistinguishablefrom a T = 0 calculation. NRG data from K. Stadler com-puted with a code of A. Weichselbaum. QMC data from M.Ferrero. term | t n − (cid:105) . Entanglement therefore builds up only inthe region where it is relevant, that is, in the vicinity ofthe impurity. Hence, matrix dimensions grow consider-ably more slowly when using a logarithmic discretizationas compared to a linear discretization. (iii) One faces afaster speed of convergence of the Chebyshev moments asin the linear case: The complexity of the spectral func-tion is considerably reduced when averaging over possiblepeaks in the high-energy structure of the spectral func-tion, as is done when using a logarithmic grid. The as-3 A () D loglin FIG. 9. (Color online) Spectral function for the two-bandHubbard model. The system parameters U/D = 1 . J/U = , U (cid:48) = U − J , n = 2 are very similar to the one in Fig. 8(a).We performed a calculation with linear (“lin”, L = 40 perbath) and one with logarithmic discretization (“log”, L =20 per bath). We fixed a truncation error ε compr = 10 − .For the calculation with logarithmic discretization, we useda scaling a = 25 D and computed N c = 300 moments. Forthe calculation with linear discretization, we used a scaling of a = 125 D and computed N c = 1250 moments. We used linearprediction in all cases. The logarithmic discretization used adiscretization parameter Λ = 2, leading to grid energies Λ − n (see e.g. Ref. 10). sociated Chebyshev expansion therefore converges morequickly than in the case of a linear grid.When using a logarithmic discretization, one hasto convolute the resulting spectral function with aGaussian to average over the finite-size features thatoriginate from the coarse log resolution at high energies.In Fig. 8, we compare exemplary calculations for thetwo-band Hubbard model with NRG and analyticallycontinued QMC data. We find good agreement in theregions around the Fermi energy, where the pinning crite-rion is respected to high accuracy without being enforced.We explain the observed disagreement far away from theFermi energy with a different specific implementation ofthe broadening convolution. One DMFT iteration for ourcalculations took around 20 min running on two 2 . A () D band 1band 2 FIG. 10. (Color online) Results for the spectral densities inthe two orbitals for J = 0. Our results are for U = 2 . D , U (cid:48) = 1 . D , n = 2 (half filling) and depicted by the solidlines. The reference NRG results are for U = 2 . D , U (cid:48) =1 . D and depicted by the dashed lines. We had to choose aslightly smaller interaction for a meaningful comparison, asfor the parameters of Greger et al. , we converged, thoughvery slowly, into an insulating solution without central peak.The non-interacting single-particle half-bandwidth of the firstband is D , and the one of the second band is 1 . D . Weused two chains of length L = 20 each, and a logarithmicdiscretization parameter of Λ = 2, leading to grid energies Λ − n (see e.g. Ref. 10). criterion is not fulfilled. The additional structure in theHubbard band, which is not visible in the calculationwith the logarithmic discretization, is seen to be similarto the one observed in Ref. 35. One DMFT iterationfor the computation that uses a logarithmic grid took20 min running on two 2 . ε compr = 2 × − . VII. CONCLUSIONS We solved several DMFT problems with two baths onthe real frequency axis with unbiased energy-resolutionbased on an DMRG impurity solver using Chebyshevpolynomials for the representation of spectral functionsat moderate numerical effort. DMRG is thereby seen tobe a viable alternative for DMFT impurity solvers also4beyond the well-understood single-impurity single-bandcase.Technically, it was crucial to apply the adaptive trun-cation scheme of Sec. III to maintain a modest numericaleffort: in all cases, the new scheme gave much better re-sults than the previously employed scheme based on fixedmatrix dimensions. Another important way of tuning thecalculation is provided by the mapping of the spectrum tothe convergence interval of Chebyshev polynomials: Thedifferent options to set up a CheMPS calculation can besummarized to yield two alternatives. (i) One uses the b = 0 setup and post-processes moments with linear pre-diction. (ii) One uses the b = − a setup and avoids linearprediction, using simple Jackson damping. Dependingon the problem, the first or the second method can bemore efficient. The second alternative is computationallymuch more efficient for cases in which linear prediction isa non-controlled extrapolation scheme, but has problemsto resolve sharp peaks at the Fermi edge.The method presented in this paper can in principlebe extended to the case of more than two baths with-out major changes to the DMFT-DMRG interface andthe Chebyshev-based impurity solver as such. However,while two baths can still be modeled by a single chainwith the impurity at the center (instead of at the end, asin single-band DMFT), this is no longer possible for threeand more baths. This will necessitate a new setup of theDMRG calculation replacing the chain-like by a star-likegeometry with the impurity at the center of the star,hence a generalization from a matrix-based to a tensor-based representation at the location of the impurity. Itremains to be seen at which numerical cost reliable re-sults on the real frequency axis will be obtainable. VIII. ACKNOWLEDGEMENTS FAW acknowledges discussions with C. Hubig and K.Stadler. US acknowledges discussions with M. Ganahland H.-G. Evertz. FAW and US acknowledge discussionswith P. Werner and support by the research unit FOR1807 of the DFG. O. P. acknowledges support from theERC Starting Grant 278472–MottMetals. We acknowl-edge K. Stadler and M. Ferrero for providing the data oftheir NRG and QMC calculations. Appendix A: Scaling of Chebyshev Moments withrespect to energy scaling The Chebyshev moments obtained by using two differ-ent scalings H (cid:48) = H/a and H (cid:48) = H/a are from (13) µ a n = (cid:80) i W i T n (( E i − E ) /a ) and µ a n = (cid:80) i W i T n (( E i − E ) /a ). As we consider one-particle operators c † theweights W i = |(cid:104) E i | c † | E (cid:105)| fulfill W i = 0 for E i with | E i − E | (cid:38) W single , (A1) where W single is the single-particle bandwidth. If thescalings a = min( a , a ) are chosen large enough, W single /a (cid:28) 1, then µ a n = µ a n if a n ∈ N and a n ∈ N . (A2)Proof: If these requirements are met, the eigenvalues E i with W i (cid:54) = 0 are close to the groundstate energy: x = ( E i − E ) /a (cid:28) 1. The Taylor expansion arccos( x ) = π/ − x − x / . . . becomes reliable already when x (cid:46) ,which is fulfilled if a is at least twice the single-particlebandwidth as in all hitherto known applications .Consider a particular energy E = E i − E for which W i > 0. It holds T a n ( E/a ) = T a n ( E/a )cos( a n arccos( E/a )) = cos( a n arccos( E/a ))cos( a n ( π/ − E/a )) (cid:39) cos( a n ( π/ − E/a )) a n ( π/ − E/a ) mod 2 π (cid:39) a n ( π/ − E/a ) mod 2 πa nπ/ π (cid:39) a nπ/ πa n/ (cid:39) a n/ . A sufficient condition for the last line to hold is that both a n/ a n/ Appendix B: Chebyshev Moments of Lorentzian andGaussian If we fix the shift to be b = 0, equation (24) is obtainedas follows. As µ n = (cid:80) i α i µ lin we only have to computethe moments for a single Lorentzian, which allows to dropthe index iµ ln = ηπ (cid:90) − dω cos( n arccos( ω ))( ω − ω ) + η (cid:39) ηπ (cid:90) − dω cos( n ( π − ω ))( ω − ω ) + η = ηπ (cid:90) − dω cos( n ( ω + ω (cid:48) )) ω + η , ω (cid:48) = ω − π ηπ Re (cid:90) − dω exp( in ( ω + ω (cid:48) )) ω + η = ηπ πi Res (cid:16) cos( in ( ω + ω (cid:48) )) ω + η (cid:17)(cid:12)(cid:12)(cid:12) ω = iη = cos( n ( ω − π e − nη . When closing the integral in the complex plane, we as-sumed that the Lorentzian concentrates almost all of itsweight within [ − , A Gauss ( ω ) = (cid:88) i α i √ πσ i e − ( ω − ωi )22 σ i , ⇒ µ gn (cid:39) (cid:88) i α i cos( n ( ω i − π e − ( σ i n ) / , (B1)as shown by a similar calculation: µ gn = 1 √ πσ (cid:90) − dω e − ( ω − ω σ cos( n arccos( ω ))= 1 √ πσ (cid:90) − dω e − ω σ cos( n ( ω + ω (cid:48) )) , ω (cid:48) = ω − π 2= 1 √ πσ Re (cid:90) − dω e − ω σ + inω + inω (cid:48) = Re e − σ n + inω (cid:48) = cos( n ( ω − π e − σ n . From the third to the fourth line, the extension of theintegral limits to ±∞ in order to apply the Gaussianintegral formula is well justified, as the Gaussian concen-trates all its weight within [ − , Appendix C: Single-bath impurity calculations1. Single-impurity Anderson Model The single impurity Anderson model (SIAM) in itstruncated chain representation is H = L − (cid:88) n =0 ,σ t n ( c † nσ c n +1 σ + h.c.) + L (cid:88) n =0 ,σ ε i n σ + U n ↓ n ↑ , (C1)with hybridization function ∆ ( z ) = t z − ε − t z − ε − · · · z − ε L − − t L − z − ε L . (C2)For an infinitely long chain, the continuous version ofthe SIAM is recovered. The bath density of states is Γ ( ω ) = − π Im ∆ ( ω + i + ). For an infinite homogeneoussystem with t i = t = D/ ε i = 0, Γ ( ω ) is the semiellipticdensity of states at half bandwidth D Γ ( ω ) = 2 πD (cid:112) − ( ω/D ) . (C3)In the non-interacting case, also the spectral function A ( ω ) is semielliptic.The computation of the spectral function A ( ω ) for theSIAM is much less demanding than for most DMFT ap-plications: A ( ω ) has only few sharp features, which in ad-dition are well approximated by Lorentzians (Sec. V 1). A () D U=2DU=D FIG. 11. (Color online) Single Impurity Anderson Model withsemi-elliptic density of states of half-bandwidth D . We com-pute the spectral function with CheMPS allowing a cumula-tive truncated weight of ε compr = 7 × − and post-processmoments with linear prediction (solid lines). These resultsare compared to data obtained with dynamic DMRG (dashedlines) by Raas et al. . We used a fermionic representationof the SIAM on a chain with length L = 80. Hence, linear prediction can be applied and we observevery good agreement with DDMRG data of Raas et al. in Fig. 11, confirming results of Ref. 31. For the case U = D , we observe a slight disagreement in the re-gion of the shoulders, where the linear prediction predictstwo small peaks, whereas DDMRG shows a perfectly flatshoulder. This might point out a failure of linear pre-diction for the description of this feature. Although thisshould be of minor importance here, it could matter inother cases. 2. Single-site single-orbital DMFT The single-site DMFT of the one orbital Hubbardmodel H = (cid:88) k σ ε k n k σ + U (cid:88) iν n i ↑ n i ↓ (C4)is well established and amounts to the determination ofthe self-consistent parameters { t i , ε i } of a SIAM (C1).We give a derivation of the DMFT equations only for themore complicated case of the cluster DMFT (Sec. D),which can easily be reduced to the single site case.Fig. 12 shows our results for which we fixed a maxi-mum cumulative truncated weight of ε compr = 5 × − .For the quite featureless spectral function of Fig. 12(a)( U = D ), the thermodynamic limit is already obtainedfor L = 40 and one DMFT iteration took 0 . U = 2 D ), we needed L = 80 and one DMFTiteration took around 3 h. For Fig. 12(b) ( U = 2 . D )we obtained converged DMFT loops, which violate the6 A () D (a) this workKarski (2008) A () D (b) this workKarski (2008) A () D (c) L =80L =120Karski (2008) FIG. 12. (Color online) Local density of states within DMFTfor the single-band Hubbard model on the Bethe lattice.Computed using CheMPS with an allowed cumulative trun-cated weight of ε compr = 5 × − . Panel (a): U = D . Panel(b): U = 2 D . Panel (c): U = 2 . D . We compare our resultswith data from Karski et al. . pinning criterion A (0) = 2 π/D , though. When employ-ing large bath sizes of L = 100 and more, we could notreach sufficiently high numbers of Chebyshev momentswithin reasonable computation times of up to 12 h perDMFT iteration; the linear prediction then overestimatesthe height of the central peak. Appendix D: Technical details of VBDMFT In this appendix, we provide the technical details forthe VBDMFT calculation. 1. Self-consistency loop The Green’s function for a patch κ has been introducedin Sec. VI A and reads G κ ( z ) = 1 | P κ | (cid:88) k ∈ P κ z + µ − ε k − Σ κ ( z ) . (D1)Within the DCA, one obtains an estimate for Σ κ ( z ) bysolving an auxiliary impurity-bath system, the Green’sfunction of which is G imp κ ( z ) − = z + µ − ∆ κ ( z ) − Σ κ ( z ) , (D2)where the bath is completely characterized by the hy-bridization function ∆ κ ( z ).The problem is then to determine ∆ κ ( z ) such that theimpurity-bath system best approximates the actual lat-tice environment, which amounts to the self-consistency condition G κ ( z ) = G imp κ ( z ) . (D3)This equation constitutes a fixed-point problem for thehybridization function ∆ ( z ) and can hence be solved it-eratively, starting with some initial guess, e.g. the non-interacting solution.Solving the impurity problem for the initial guess of ∆ ( z ), one obtains G imp κ ( z ). From that one obtains theestimate for the self-energy as Σ κ ( z ) = G imp0 κ ( z ) − − G imp κ ( z ) − , or by the method of Bulla et al. (we foundthe latter not to yield advantages for the CheMPS setup).The self-energy is then inserted into (D1) to obtain a newvalue for G κ ( z ). Using self-consistency, this defines a newhybridization function by inserting (D3) in (D2): ∆ κ ( z ) = − G κ ( z ) − + z + µ − ε − Σ κ ( z ) . (D4)In QMC calculations, one defines all quantities on theimaginary axis. In this work as in NRG calculations, wedefine all quantities on the real axis: the spectral densityof the bath is Γ ( ω ) = − π Im ∆ ( ω + i + ) , (D5)7which leads to slightly modified version of (D4) Γ κ ( ω ) = 1 π Im( G κ ( ω ) − + Σ κ ( ω )) . (D6)If one considers ordinary single-site DMFT, all equa-tions remain the same and the momentum patch index κ can be dropped. In a multi-band calculation, the in-dex κ plays the role of the band index. For DMFTcarried out for the Bethe lattice, self-consistency canbe written as Γ ( ω ) = D A imp ( ω ), where A imp ( ω ) = − π Im G imp ( ω + i + ). An iterative solution is partic-ularly simple in this case, as only the spectral functionhas to be computed and summations over k space are notnecessary. In the general case, also the real part of theGreen’s function is needed. This can either be accessedfrom the spectral function by the Kramers-Kronig rela-tion or directly from the Chebyshev moments through G imp ( ω ) = − ia (cid:88) n w n ( ω (cid:48) ) µ n exp( − in arccos( ω (cid:48) )) (D7)where ω (cid:48) ≡ ω (cid:48) ( ω ) is the rescaled frequency definedin (11). The preceding equation should be evaluatedslightly away from the real axis ω (cid:48) → ω (cid:48) + i + .In our computations, we parallelized the independentcomputations for the particle and the hole part of theGreen’s (spectral) function, as well as those for differentimpurity sites. 2. Bath discretization In order to represent the continuous hybridizationfunction ∆ ( z ) using a discrete chain, we use the gen-eral procedure of Bulla et al. (in the notation of Ref.47) adding details for the special case of the linear dis-cretization.If we know the hybridization function Γ ( ω ) (D5) onthe real axis, the bath and coupling Hamiltonian can bewritten as H b = (cid:90) − dε εa † ε a ε + (cid:90) − dε (cid:112) Γ ( ε )( d † a ε + h.c.) (D8)We discretize the Hamiltonian using a linear discretiza-tion of the bath energies I n = [ (cid:15) n , (cid:15) n +1 ] , (D9) (cid:15) n = n∆(cid:15) + (cid:15) for n ∈ { , , . . . , L b } . For a given bath size L b , we fix the free parameters (cid:15) and ∆(cid:15) by requiring (cid:82) (cid:15) Lb (cid:15) dωΓ ( ω ) = 0 . (cid:82) ∞−∞ dωΓ ( ω ). Thisleads to outer interval borders (cid:15) and (cid:15) L b that are closeenough to minimize finite-size effects, and far enoughapart from each other, to contain almost the completesupport of Γ ( ω ). Starting with an interval [ (cid:15) init0 , (cid:15) init L b ]that contains the full integrated weight of Γ ( ω ), we re-peatedly shift the boundaries by a fixed small number to shrink it down to the required size. In a single step, wechoose the boundary, that can be shifted with a smallerreduction of the total integral weight. The boundary thatleads to a higher reduction is left unchanged in this step.When using a logarithmic discretization, we defined thediscretization intervals via energies (cid:15) m ∝ ± Λ − m , where m ∈ [1 , ..., L b / The specific choice of boundaries ofthe support is not of much importance in this case.The discretized SIAM then couples to L b bath statescreated by a † n each of which corresponds to a bath energyinterval I n . One approximates the continuous H b by thediscrete version H b (cid:39) L b (cid:88) n =1 ξ n a † n a n + L b (cid:88) n =1 γ n ( d † a n + h.c.) .γ n = (cid:90) I n dε Γ ( ε ) , ξ n = 1 γ n (cid:90) I n dε εΓ ( ε ) . In order to use an MPS representation, one has tomap the preceding Hamiltonian on a chain Hamiltonian.This is done using the Lanczos algorithm with high-precision arithmetics for the diagonal quadratic matrix( ξ n δ nm ) L b n,m =1 applied to the initial vector ( γ n ) L b n =1 . After L b Lanczos iterations one obtains the site potentials ε i asthe diagonal of the tridiagonal Lanczos matrix, and thehopping terms as the side-diagonal entries t i . The hop-ping term from the impurity site to the first bath chainsite is the square root of the total hybridization magni-tude t = (cid:80) n γ n = (cid:82) dεΓ ( ε ). With these definitions, thefinal chain Hamiltonian reads H b (cid:39) L b (cid:88) i =0 t i ( c † i +1 c i + h.c.) + L b (cid:88) i =1 ε i c † i c i , (D10)where the impurity site is the first site of the chain c † ≡ d † .An alternative method to directly obtain the bath pa-rameters by truncating the continued fraction expansionof the hybridization function as put forward by Karski et al. , did not show any advantages but led to equiv-alent results. As the method of Karski et al. leads tohopping energies that converge to a constant far awayfrom the impurity, while the linear discretization schemeleads to polynomially decreasing hopping energies, thelinear discretization method leads to a smaller many-body bandwidth. This allows to use smaller rescalingvalues in CheMPS. 3. Finding the ground-state The first problem to solve is finding the ground stateof the model Hamiltonian. a. Initializing the wave function For the two-chain layout (29) of the model, the follow-ing problem arises: the chemical potential of both chains8can be strongly different, in which case the particle num-bers on the left N κ =+ and the right N κ = − chain maybe strongly different. Note that the Hamiltonian of (29)commutes with N κ =+ and N κ = − , as the chains are merelycoupled by an interaction, not a hopping term. If startinga DMRG groundstate search with a global random statefor such a system, convergence can be expected to bevery slow, as the local optimization does not pick up theglobal potential variation. Even worse, the absence of anhopping term between the two chains prevents that dur-ing minimization the particle numbers in the left N κ =+ and the right N κ = − chain change. This can in principlebe compensated by choosing White’s mixing factor tobe large when starting to sweep, reducing it when beingclose to convergence. But still we found it impossibleto implement a reliable automatized groundstate searchunder these circumstances.The problem can be solved by using a U = 0 solutionas initial guess for the groundstate search. One shouldrealize that the partition between N − and N + = N − N − (where N is the total particle number) only weakly de-pends on the interaction U : The total potential andhopping energies scale with the bath length, whereasthe interaction energy is a single-site quantity. Giventhe system parameters { ε κi } and { t κi } for each chain κ , we diagonalize the L = L b + 1 dimensional tridiago-nal single-particle representation of a single chain withits associated impurity site. This gives us the particlesectors N ± of the groundstate of each subsystem. The U = 0 estimate for the total particle number sector is N = N + + N − , as in this case both subsystems are un-coupled. Given an initial guess for the chemical potential µ , one should initialize a wave function that fulfills the U = 0 estimates for N and N + /N − . b. Finding the correct symmetry sector As the DMFT is grand-canonical, one still needs tosolve the problem of finding the correct particle numbersector for the DMRG calculation. This can be greatlyaccelerated using the U = 0 estimate for N , which con-stitutes a rigorous upper bound for the particle numberin the interacting system. For a given µ one can there-fore use a bisection search, starting with N , N − ∆N and N − ∆N . In case N − ∆N yields the lowest energy esti-mate, one has to extend the search regime to lower valuesof N . If N or N − ∆N yield the lowest energy, one cancontinue the ordinary bisection search. For typical inter-action values, ∆N/N = 0 . 05 is a meaningful choice. Ifsearching for the maximum energy state, which is neces-sary if one wants to determine the full many-body band-width W = E max − E , one searches for the groundstateof − H . In this case the interaction between electrons be-comes attractive, and the U = 0 solution for the particlenumber sector of | E max (cid:105) becomes a rigorous lower boundfor the interacting system.Having found the correct symmetry sector together with its groundstate for a given value of µ , one has tocheck whether the requirements for the local impuritydensities are fulfilled n − (cid:88) κ (cid:104) c † κ c κ (cid:105) ? = 0 . (D11)To find the correct value of the chemical potential, a sim-ple update of the chemical potential µ with the residuumof (D11) is usually not sufficient to achieve convergence.Instead, we use this method until we found a lower andupper bound for µ and then use a bisection again.In some cases, the algorithm has to break its searchbefore reaching the required tolerance. This is when thedesired chemical potential lies directly on the boundarywhich separates two different particle number sectors. Ifthis is the case, due to the discrete nature of our model,no solution can be found. Such a case is typically de-tected by observing oscillations in the residuum of (D11).When setting up the groundstate search naively, it caneasily take most of the computation time of the calcu-lation. Using the procedures just described, it usuallytakes only a negligible few percent of the total computa-tion time. 4. Definition of the model Hamiltonian In the following, we outline the standard procedurethat eliminates the high-energy contributions in the hy-bridization function.We want to represent the non-interacting patchGreen’s function G κ ( z ) = 1 | P κ | (cid:88) κ ∈ P κ z + µ − ε k , (D12)by an impurity model with Green’s function G imp0 κ ( z ) = z + µ − ∆ ( z ) , such that G κ ( z ) = G imp0 κ ( z ) . (D13)When defining the bath hybridization function naivelyvia ∆ κ ( z ) = z + µ − G − κ ( z ) , (D14)one observes that ∆ κ ( z ) → t κ for | z | → ∞ , when ex-panding for high values of | z | , as G κ ( z ) = 1 z + µ (cid:18) t κ z + µ + O ( z − ) (cid:19) , (D15a) G − κ ( z ) = z + µ − t κ + O ( z − ) , (D15b)where t κ = | P κ | (cid:80) k ∈ P κ ε k .This means that the corresponding spectral density ofthe bath Γ ( ω ) = − π Im ∆ ( ω + i + ) has contributions atarbitrarily high energies and the discretization procedure9that maps Γ ( ω ) onto the discrete bath Hamiltonian H b must fail.This problem is solved by defining an impurity modelat a shifted chemical potential µ → µ − t κ . In the hy- bridization function of this shifted impurity model ∆ κ ( z ) = z + µ − t κ − G − κ ( z ) , (D16)the constant t κ in the high-energy expansion of G − κ ( z )(D15b) cancels out. It therefore approaches zero for | z | →∞ while still fulfilling (D13) for G imp0 κ ( z ) = z + µ − t κ − ∆ ( z ) .As t κ is a simple constant shift of the chemical poten-tial, one can as well incorporate it into the Hamiltoniandescription of the impurity model, as done in (29). W. Metzner and D. Vollhardt, Physical Review Letters ,324 (1989). A. Georges and G. Kotliar, Phys. Rev. B , 6479 (1992). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). G. Kotliar, S. Savrasov, K. Haule, V. Oudovenko, O. Par-collet, and C. Marianetti, Reviews of Modern Physics ,865 (2006). T. Maier, M. Jarrell, T. Pruschke, and M. Hettler, Rev.Mod. Phys. , 1027 (2005). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,M. Troyer, and P. Werner, Rev. Mod. Phys. , 349(2011). A. Rubtsov, V. Savkin, and A. Lichtenstein, Phys. Rev.B , 035122 (2005). E. Gull, P. Werner, O. Parcollet, and M. Troyer, EPL(Europhysics Letters) , 57003 (2008). P. Werner, A. Comanac, L. de’ Medici, M. Troyer,and A. Millis, Physical Review Letters (2006),10.1103/physrevlett.97.076405. R. Bulla, T. Costi, and T. Pruschke, Rev. Mod. Phys. ,395 (2008). M. Caffarel and W. Krauth, Phys. Rev. Lett. , 1545(1994). M. Granath and H. U. R. Strand, Phys. Rev. B , 115111(2012). Y. Lu, M. H¨oppner, O. Gunnarsson, and M. W. Haverkort,ArXiv , 1402.0807 (2014), 1402.0807. A. K. Mitchell, M. R. Galpin, S. Wilson-Fletcher, D. E.Logan, and R. Bulla, Phys. Rev. B , 121105 (2014). S. R. White, Phys. Rev. Lett. , 2863 (1992). U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005). U. Schollw¨ock, Annals of Physics , 96 (2011). K. A. Hallberg, Phys. Rev. B , R9827 (1995). D. Garc´ıa, K. Hallberg, and M. Rozenberg, Phys. Rev.Lett. , 246403 (2004). P. E. Dargel, A. W¨ollert, A. Honecker, I. P. McCulloch,U. Schollw¨ock, and T. Pruschke, Phys. Rev. B , 205119(2012). F. A. Wolf, unpublished (2013). T. K¨uhner and S. White, Phys. Rev. B , 335 (1999). E. Jeckelmann, Physical Review B (2002),10.1103/PhysRevB.66.045114. M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B ,075116 (2008). M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B ,113110 (2005). S. Nishimoto and E. Jeckelmann, J. Phys.: Condens. Mat- ter , 613 (2004). A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev.Mod. Phys. , 275 (2006). A. Holzner, A. Weichselbaum, I. P. McCulloch,U. Schollw¨ock, and J. von Delft, Phys. Rev. B , 195115(2011). A. Braun and P. Schmitteckert, ArXiv , 1310.2724 (2013),1310.2724. A. C. Tiegel, S. R. Manmana, T. Pruschke, and A. Ho-necker, ArXiv , 1312.6044 (2013), 1312.6044. M. Ganahl, P. Thunstr¨om, F. Verstraete, K. Held, andH. G. Evertz, ArXiv (2014), 1403.1209. A. J. Daley, C. Kollath, U. Schollw¨ock, and G. Vidal, J.Stat. Mech.: Theor. Exp. , P04005 (2004). G. Vidal, Phys. Rev. Lett. , 040502 (2004). S. R. White and F. A. E., Phys. Rev. Lett. , 076401(2004). M. Ganahl, M. Aichhorn, P. Thunstr¨om, K. Held, H. G.Evertz, and F. Verstraete, ArXiv , 1405.6728 (2014),1405.6728. L. Lin, Y. Saad, and C. Yang, ArXiv (2013), 1308.5467. J. B. Boyd, Chebyshev and Fourier Spectral Methods (Dover Publications, Mineola, New York, 2001). F. Verstraete and J. Cirac, Phys. Rev. B , 094423 (2006). C. Raas, G. S. Uhrig, and F. B. Anders, Phys. Rev. B ,041102 (2004). W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery, Numerical Recipes 3rd Edition: The Art of Sci-entific Computing , 3rd ed. (Cambridge University Press,New York, NY, USA, 2007). S. White and I. Affleck, Phys. Rev. B , 134437 (2008). T. Barthel, U. Schollw¨ock, and S. R. White, Phys. Rev.B , 245101 (2009). R. N. Silver and H. R¨oder, Phys. Rev. E , 4822 (1997). M. Ferrero, P. Cornaglia, L. De Leo, O. Parcollet,G. Kotliar, and A. Georges, Physical Review B , 064501(2009). M. Ferrero, O. Parcollet, A. Georges, G. Kotliar, and D. N.Basov, Physical Review B , 054502 (2010). A. Weichselbaum and J. von Delft, Phys. Rev. Lett. ,076402 (2007). K. Stadler, Towards exploiting non-abelian symmetriesin the Dynamical Mean-Field Theory using the Numeri-cal Renormalization Group , Master’s thesis, LMU Munich(2013). A. Weichselbaum, Annals of Physics , 2972 (2012). M. Ferrero, Private Communication. M. Greger, M. Kollar, and D. Vollhardt, Physical Review Letters , 046403 (2013). C. Raas, Dynamic Density-Matrix Renormalization for theSymmetric Single Impurity Anderson Model , Ph.D. thesis,University of Cologne (2005). R. Bulla, A. C. Hewson, and T. Pruschke, Journal ofPhysics: Condensed Matter , 8365 (1998). S. R. White, Phys. Rev. B72