Imaginary-time matrix product state impurity solver for dynamical mean-field theory

We present a new impurity solver for dynamical mean-field theory based on imaginary-time evolution of matrix product states. This converges the self-consistency loop on the imaginary-frequency axis and obtains real-frequency information in a final real-time evolution. Relative to computations on the real-frequency axis, required bath sizes are much smaller and less entanglement is generated, so much larger systems can be studied. The power of the method is demonstrated by solutions of a three band model in the single and two-site dynamical mean-field approximation. Technical issues are discussed, including details of the method, efficiency as compared to other matrix product state based impurity solvers, bath construction and its relation to real-frequency computations and the analytic continuation problem of quantum Monte Carlo, the choice of basis in dynamical cluster approximation, and perspectives for off-diagonal hybridization functions.


I. INTRODUCTION
Dynamical mean-field theory (DMFT) in its singlesite [1][2][3] and cluster [4,5] variants is among the most widely employed computational techniques for solving quantum many-body problems.At the core of a numerical solution of DMFT is an impurity solver : an algorithm for solving a quantum impurity problem.The most prominent examples of impurity solvers are continuous time quantum Monte Carlo (CTQMC) [6][7][8], exact diagonalization (ED) [9][10][11], the numerical renormalization group (NRG) [12], and the density matrix renormalization group (DMRG) [13].While all methods have their strengths, key limitations mean that fundamentally important classes of problems have not yet been adequately addressed.Other recent suggestions for impurity solvers [14][15][16][17][18][19] including in particular the computationally inexpensive density matrix embedding theory (DMET) [20], are promising but have not been tested in detail.
CTQMC is widely employed but its application to situations involving low point symmetry, non-Hubbard interactions or multiple relevant orbitals is limited by the fermionic sign problem.Reaching low temperatures becomes highly computationally expensive while calculating real-frequency information requires analytical continuation, a numerically ill-posed procedure fraught with practical difficulties.
ED makes no assumption on the interaction and does not have a sign problem.It is limited by the size of the Hilbert space that can be studied, meaning in practice that it is restricted to a small number of correlated sites to which only a small number of bath sites can be attached.Recently, improvements have been achieved by considering only restricted subspaces of the Hilbert space [21][22][23][24], but the size of problem remains a significant limitation.
NRG converges the DMFT loop on the real-frequency axis and very effectively obtains real-frequency information in the low frequency limit.Current applications have been to relatively small problems (the most recent achievement is a solution of the single-site DMFT approximation to a three-band model [25]) and it remains to be seen how far the method can be extended.
DMRG [26] is a set of algorithms operating on the space of matrix product states (MPS) [27].It has been found to be extremely powerful for the calculation of ground states of one-dimensional quantum systems [27,28]; it was very successfully extended to the calculation of spectral functions which, in contrast to NRG, it obtains with equal resolution across the spectrum, see e.g.[29,30].In pioneering work the method was applied as a DMFT solver by García, Hallberg, and Rozenberg [13] and Nishimoto, Gebhard, and Jeckelmann [31] with important further work done by these and other authors [32][33][34][35][36][37][38].However, the method has not been widely accepted, perhaps because high-quality data were presented only for the single-site approximation to the single-band Hubbard model.Recently the method was used to reliably solve a two-site DCA [39], and insights into the entanglement of the impurity problem made it even more powerful [40].In view of these advances, DMRG now is a candidate for a highly flexible low-cost impurity solver, which can in addition be efficiently employed in the nonequilibrium formulation of DMFT [40][41][42].
This paper reformulates the DMRG method as an impurity solver for DMFT on the imaginary-time axis (Sec.II).As we will show, this strongly reduces entanglement and necessary bath sizes.The price to be paid is a reduced resolution on the real-frequency axis, which we study in detail by comparing with calculations that converge the DMFT loop on the real-frequency axis (Sec.III).The reformulation allows a much larger class of problems to be addressed, including some that are unreachable by other methods, due e.g. to the sign problem, the size of the correlated cluster or the number of bands.We illustrate this with calculations for three band models in the single-site and two-site DMFT approximation (Sec.IV).We append discussions of the optimization of typical DMFT Hamiltonians (Appendix A) and of the entanglement in different representations of the DCA including a discussion of off-diagonal hybridization functions (Appendix B).

II. METHOD
A. Overview: Green's functions in DMRG The computational key challenge in DMFT is the computation of the full frequency dependence of the Green's function of a quantum impurity model involving an essentially arbitrary bath.The "size" (number of correlated sites L c ) of the impurity model should be as large as possible and the kinds of interaction that can be treated should be as general as possible.The Green's function is used in a self-consistency loop, which may require many iterations for convergence.The solution should be as inexpensive as feasible, and must run automatically, without need for manual optimization of parameters or procedures.In this subsection we present a qualitative discussion of the issues involved in computing the Green's function using DMRG methods, to motivate the work described in detail below.
Within DMRG one computes Green's functions by first representing the system ground state |E 0 as an MPS.One then generates a one electron (one hole) excitation is at most as entangled as the ground state |E 0 [30], in order to compute a Green's function, one has to perform further operations on |ψ ≷ 0 .These operations typically increase entanglement and by that the bond dimension of an MPS, which ultimately limits all computations.
Let us be more concrete and consider a general MPS of bond dimension m for a system with L sites and open boundary conditions.Defining A σi , B σi ∈ C m×m for i = 1, L and A σ1 ∈ C 1×m , B σ L ∈ C m×1 , where σ i ∈ {0, ↑, ↓ , ↑↓} labels a local basis state of the Hilbert space, any MPS can be represented as [27] where S = diag(s 1 , ..., s m ) is a diagonal matrix, and A σi are left-normalized and B σi are right-normalized, respectively: Here, I are identity matrices.Left-and rightnormalization make Eq.( 1) the Schmidt decomposition of |ψ MPS that is associated with partitioning the system at bond (l, l + 1).The bond entanglement entropy for the associated reduced density matrix can therefore simply be read of from Eq. ( 1) [27] When subsequently we refer to an entanglement growth associated with repeated operations on |ψ MPS this implies the need to adjust the bond dimension m such that |ψ MPS still faithfully represents a physical state.If entanglement in the physical state becomes too large, we have to choose m so large that computations with MPS become impractical.Since the first suggestion for computing spectral functions within DMRG [43], the field has evolved by the important development of the correction vector method [44,45].The subsequent understanding of the connection between DMRG and MPS [27] opened the door to many further approaches to computing spectral and Green's functions, in particular, time evolution and subsequent Fourier transform [46,47], an improved Lanczos algorithm [48], and the Chebyshev recursion [29,30,49].All of these are formulated for the calculation of spectral functions at T = 0, as considered in the present paper, and came at much cheaper computational cost than the correction vector method [29,30].We note that for T > 0, there are perspectives for even more powerful algorithms: it was recently demonstrated that the numerically exact spectral function of a molecule consisting of several hundreds of interacting spins could be computed [50].
These developments (see Appendix C for more details) make MPS-based solvers an attractive possibility for dynamical mean-field theory.However, the growth of entanglement arising in all calculations of Green's function has limited the systems sizes that have been addressed to date.Also, in MPS computations manual adjustments, for example choosing optimal broadening [29] or combining results of different systems sizes [48], are still common practice.In the rest of this section, we show that these problems can to a large degree be circumvented by computing Matsubara Green's functions using imaginarytime evolution.The imaginary-time framework naturally extends existing approaches based on real-time evolution [38,40], which have recently been shown to provide currently the most efficient algorithmic approach to compute real-frequency spectral functions [30].

B. Imaginary-time computation
The central objects of technical interest in this paper are the greater and the lesser correlation functions G ≷ , which we define for imaginary time τ In the second line, we evaluate G ≷ (τ )| τ =it and by that obtain a correlation function for real time t, which will  12) for β eff = 315/D and α = 0.For the ground state optimization, we enforced a maximal bond dimension of m = 300.The Krylov time evolution used a time step of ∆t = 0.1/D and allowed for a maximal global truncation error of 10 −4 at each time step adjusting bond dimensions automatically.This leads to an immediate decay of m at τ 0 from m = 300 down to m 110, as seen in panel (b).We used the global SU(2) symmetry of the Hamiltonian to reach these low values of the bond dimension.
be useful later on.The functions G ≷ carry spin and orbital indices associated with the spin and orbital indices of the single-particle (hole) excitation |ψ ≷ 0 , but these indices are not explicitly written here.We will discuss the relationship of G ≷ to the physical Green's functions (which we denote by G) below.
While it is not essential in principle, we evaluate Eq. (4) using a Krylov algorithm [52], which represents the time evolution operator in a local Krylov space and is able to treat Hamiltonians with long-ranged interactions.Before performing a time evolution computation, one has to compute the initial state |ψ ≷ 0 using an MPS optimization of the ground state.As impurity models come with open boundary conditions, this is well suited for DMRG.We discuss this optimization for typical DMFT Hamiltonians in Appendix A 1.
Figure 1 presents representative results based on parameters obtained from a two-site DMFT solution of the Hubbard model.Figure 1(a) shows the time evolution of G ≷ (τ ) out to times as long as 350 times the basic timescale (inverse half-bandwidth D) of the model, which suffices to converge G ≷ (τ ) to a precision of 5•10 −4 .Figure 1(b) demonstrates the key advantage that makes this computation possible: the lack of growth of maximal bond dimensions m with time of the associated imaginary-time evolved states |ψ ≷ (τ ) = e −(H−E0)τ |ψ ≷ 0 .The imaginary time-evolution operator does not create entanglement as it projects on the lowly entangled ground state.
Figure 1(a) reveals additional information about the nature and rate of convergence of G ≷ (τ ).In the insulating phase, H has a gap and G ≷ (τ ) decays exponentially irrespective of whether one considers a finite system or the thermodynamic limit.In the metallic phase, G ≷ (τ ) decays algebraically in the thermodynamic limit.For a finite system though, there always remains a small gap, and even though the decay resembles an algebraic decay for short times, it always becomes exponential at long times.The exponential decay can be exploited to speed up computations considerably by a simple technique known as linear prediction [30,47,53,54].This technique is an efficient formulation of the fitting problem for the ansatz function f (τ ) = n α n e βnτ , α n , β n ∈ C, τ ∈ R, which can then be used to reliably extrapolate functions with an exponentially decaying envelope.This is illustrated by the dashed black line in Fig. 1(a), which has been fitted to match G ≷ (τ ) for τ D ∈ [150, 200] and was then extrapolated for higher times.The solid green line, by contrast, is the result of the MPS computation.Agreement can be seen to be perfect.

C. Physical Green's functions
Of particular interest in the rest of this paper are the imaginary time Green's functions G mat (τ ) defined via whose Fourier transform gives the Matsubara Green's function at zero temperature, where ω n = (2n + 1)π/β and β → ∞.We shall also be interested in the retarded real-time Green's function from which the retarded frequency-dependent Green's function is obtained as In numerical practice, we evaluate the Fourier transforms leading to ( 6) and ( 8) approximately as with cutoff times τ max and t max .This approximation is only controlled if we are able to reach long enough times τ max and t max , such that G ≷ (τ ) and G ≷ (it) have converged to zero to any desired accuracy.
In contrast to a computation on the imaginary axis, reaching arbitrarily long times t max on the real axis is prohibited by a logarithmic growth of entanglement, which comes with a power-law growth of bond dimensions.In addition, finite-size effects are a severe source of errors because the long-time behavior is determined by the bath size.For a numerically exact computation, one has to choose the system large enough to observe exponential "pseudo-convergence" of G ≷ (it) to zero before finite-size effects are resolved [30].In the context of the present paper, we deal with small system sizes and will never observe pseudo-convergence.In particular there is no exponential pseudo-convergence, so that linear prediction cannot be employed [30].Therefore, when computing the real-frequency spectral function after converging the DMFT loop, one has to use the further approximation of damping the finite-size effects that emerge at long times by computing, instead of G ret (ω) in (9), which yields the broadened spectral function . Instead of a Gaussian damping and broadening, one could also use an exponential damping leading to Lorentzian broadening, which damps out the original time-evolution information more strongly, though.
Before presenting detailed benchmark results for the solution of DMFT using imaginary-time evolution of MPS, let us clarify which price we have to pay for profiting from the great advantage of not facing entanglement growth.We do this by comparing the imaginary-time approach (itMPS) to approaches that solve the DMFT loop on the real axis.

III. COMPARISON OF IMAGINARY-AXIS WITH REAL-AXIS COMPUTATIONS
The self-consistency equation in DMFT relates an impurity model specified by a hybridization function and a 0 0.5 1 1.5 2 2.5 3 3.  self energy to a lattice model specified by a lattice Hamiltonian and the same self energy.We discuss the issues using the example of the dynamical cluster approximation to the single-band Hubbard model Here ε k denotes the single-particle dispersion of the lattice and µ is the chemical potential.In the dynamical cluster approximation, the Brillouin zone, consisting in N momentum vectors k, is covered by N c (for singleband L c = N c ) equal-area tiles (patches), labelled here by P K and the self energy is taken to be piecewise constant, with being a potentially different function of frequency in each tile.The impurity model is specified by the on-site energy ε K and the hybridization function Λ K (z), which is to be determined using the fixed point iteration referred to as the DMFT loop.The loop works as follows: make an initial guess for Λ K (z); compute ] −1 and repeat this procedure until convergence.
We discuss two aspects of the comparison of realand imaginary-frequency solutions of the DMFT selfconsistency equation (11).The first has to do with the number of bath sites needed to obtain a solution of the self-consistency equation.The second is the accuracy to which the spectral functions of physical interest can be reproduced.
The DMFT self-consistency equation (11), defines the hybridization function Λ K as a continuous function in terms of the difference between the computed self energy and the inverse of the lattice Green's function.In DMRG-type methods, the hybridization function Λ K is approximated as the hybridization function Λ discr K of a discrete impurity model, which is the sum of a number of poles.If the number of poles is too small one can-not construct a meaningful approximation on the real axis [55] and a DMFT loop cannot be converged.For this reason DMRG-based solutions of DMFT up to now [13,[31][32][33][34][35][36][37][38][39][40], all of which were real axis computations, have been performed using numbers of bath sites of at least L b /L c 30, and in the case of the single-band Hubbard model, even much more L b /L c 120. Use of such a large number of bath sites means that with modest broadening the hybridization function can be reasonably approximated as a continuum, enabling a stable solution of Eq. (11).
By contrast, formulating the problem on the imaginary axis (as is typically done in standard ED solvers where the number of bath sites is strictly limited) automatically smoothens the hybridization function Λ discr K and permits a stable solution.From the imaginary axis solution one must then determine the discrete set of bath parameters to represent Λ discr K .This is typically done [9,11,56] by numerical minimization of a cost function defined as Here, α defines a weighting function ω −α n .Choosing α > 0, e.g.α = 1, attributes more weight to smaller frequencies [11,56,57], which we find helpful when using small bath sizes L b /L c < 5. To define the frequency grid for the fit ω n = (2n + 1)π/β eff , one defines a fictitious inverse temperature β eff , which has no physical significance.We further employ a cutoff frequency ω c , which implies a finite number N fit of fitted Matsubara frequencies.
If one tries to define an analogous cost function for the real axis, the result is useless as then Λ discr K (ω + i0 + ) is a sum of poles, whereas the hybridization function Λ K (ω + i0 + ), as encountered in (11), is continuous [55].One can overcome this problem only when using a Lindbladt formalism [58], which increases the complexity of the problem substantially.
The minimization of ( 12) is done using using standard numerical optimization.The optimization in the initial DMFT iteration should be done using a global optimization scheme [59], and in subsequent iterations using a local optimization scheme (e.g.conjugate gradient), which takes as an initial guess for the new bath parameters the values of the previous iteration.Figure 2 shows the convergence of the fit of the hybridization function with the number of bath sites L b /L c .For L b /L c = 7 one already obtains errors as little as 10 −3 and for values L b /L c 9 the quality of the fit already stops improving.It is at this point, where we (and all ED-like techniques) face the problem of "analytic continuation" encountered in CTQMC, namely that Green's functions on the imaginary axis encode information in a much less usable form than on the real axis.
Consider again the example of the two-site DCA for the single-band Hubbard model on the square lattice.In Ref. 39, this problem has been solved entirely on the real axis using L b /L c = 39 bath orbitals.Here, we converge the DMFT loop on the imaginary axis and compute the spectral function in a final real-time evolution using L b /L c = 3, 5, 7 bath orbitals.We compare both solutions in Fig. 3. Whereas for the (central) momentum patch "+" shown in Fig. 3(a), we find satisfactory agreement of the imaginary-axis with the real-axis calculation, this is not the case for the (outer) momentum patch "−" shown in Fig. 3(b), even though the corresponding imaginaryaxis Green's function is well reproduced, see Fig. 3(d).Evidently, in Fig. 3(b), the central peak and the pseudogap at the Fermi edge are smeared out by a broadening η = 0.2D that hides finite-size effects to a large degree.Reducing the broadening to η = 0.05D as shown in Fig. 3(c), again reveals the pseudo-gap and the central peak; but together with unphysical finite-size effects.
We observe that the nature of these finite-size effects is qualitatively comparable when using different numbers of bath sites L b /L c = 3, 5, 7. Already for L b /L c = 3, we obtain reasonable results.On the imaginary axis, by contrast, L b /L c = 5, 7 still improve over L b /L c = 3 and almost agree with the numerically exact QMC data for β = 200/D of Ref. 51, see Fig. 3(d).However we emphasize that even with the modest number of bath sites used here the basic features of the spectral function are reproduced (for example the areas in given frequency ranges).

IV. THREE-BAND CALCULATIONS
A. Three-band model in single-site DMFT We now demonstrate the power of the method by applying it to three-band problems in the single-site approximation (where comparison to existing calculations can be made) and the two-site approximation.Both was hitherto not accessible to DMFT+DMRG computations.
We study the three-band Hubbard-Kanamori model with Hamiltonian (omitting the site index i in the following definition of H loc,i ) where i labels sites in a lattice and k wave vectors in the first Brillouin zone, n i,a,σ = d † i,a,σ d i,a,σ is the density of electrons of spin σ in orbital a on site i, µ is the chemical potential, ∆ a is a level shift for orbital a, ε ab k is the band  13) using the semi-elliptic density of states Eq. ( 14) and U = 12t, obtained from single-site DMFT approximation evaluated using imaginary MPS (crosses) and CTQMC data (circles, Figure 1 of Ref. 61, inverse temperature β = 50/t).In the DMRG computations the bath fitting was performed using β eff = 100/t, ωc = 6t and α = 1 with three bath sites per correlated site (L b /Lc = 3).The maximal matrix dimensions was m = 300 for the ground state calculation, exploiting the SU(2) symmetry, which leads to the high precision (H − E) 2  10 −14 .For the time evolution, we computed G ≷ a (τ ) in (4a) in steps of ∆τ = 0.1/t allowing a global truncation error of 5 • 10 −4 per step, up to imaginary time τmax = 100/t and used linear prediction for higher times.dispersion, U is the intra-orbital and U the inter-orbital Coulomb interaction, and J is the coefficient of the Hund coupling and pair-hopping terms.We adopt the conventional choice of parameters, U = U − 2J, which follows from symmetry considerations for d-orbitals in free space and holds (at least for reasonably symmetric situations) for the t 2g manifold in solids [60].
We study the orbital-diagonal and orbital-degenerate case (∆ a = 0) on the Bethe lattice, i.e. the noninteracting density of states is semi-elliptic, In the single-site approximation, the impurity Hamiltonian used within DMFT is given by where c † l,a,σ creates a fermion in the bath orbital l, V l,a,σ describes the coupling of the impurity to the orbital l, and ε l,a.σ denotes the potential energy of orbital l.The hybridization function is then given by Figure 4 compares the dependence of the particle density n on the chemical potential µ obtained by the MPS methods used here to those obtained by numerically exact CT-QMC methods [61].The plateaus in n(µ) are the Mott insulating regimes of the phase diagram.The agreement is very good in general, confirming the reliability of our new procedure even with only three bath sites per correlated site.This leads to an extremely cheap computation, for which a single iteration of the DMFT loop took about 30 min on two 2.8 GHz cores (see Appendix A 2 for more details).
In panel (a) of Fig. 5 we show a more stringent test, namely the dependence of the self energy on Matsubara frequency, in a parameter regime where the self energy was previously found [62] to exhibit an anomalous ω frequency dependence and (in some regimes) a non-zero intercept as ω → 0. These phenomena are associated with a spin-freezing transition [61,62].
Panel (a) of Figure 5 shows that the known low frequency ω t self energy is already accurately reproduced even for the computationally inexpensive choice of L b /L c = 3, although one observes deviations for the highfrequency behavior.The deviations at high frequency decrease as the number of bath sites is increased, although full convergence at all frequencies has not been demonstrated.Panel (b) of Fig. 5 shows that the deviations are linked to the impossibility of fitting the hybridization function equally well for all frequencies using only a small number of bath sites.The large deviations at high frequencies are due to the choice α = 1 in (12), which enforces good agreement for low frequencies and allowed to successfully reproduce the MIT transition in Fig. 4. Increasing the number of bath sites to L b /L c = 5 leads to a much better approximation of the hybridization function also for high frequencies, with concomitant improvement in the self-energy (Fig. 5(a)).

B. Three-band model in two-site DCA
We now present results obtained using a two-site DCA approximation to the three-band model of the previous section.For this problem there are no low-temperature results available in the literature.The size of the problem is beyond the scope of standard ED.The truncated configuration interaction (CI) impurity solver [20] allows to access a relatively high number of bath sites but is limited in the number of correlated sites: e.g. in Ref. 24 15) for U = 8t and J = U/6.Crosses (color online) represent itMPS data and black circles depict CTQMC data from Figure 3 of Ref. 62, computed at inverse temperature β = 100/t.We choose all parameters as described in the caption of Fig. 4, in particular, for the bath fitting (12), we use β eff = 100/t, ωc = 6t and α = 1.Choosing α = 1 enforces agreement for low frequencies at the price of disagreement at high frequencies, which is observed in both panels (a) and (b).In panel (b), Λ denotes the hybridization function that is fitted with the hybridization function Λ discr of the discrete impurity model.
the CI solver.The problem is also challenging for standard CTQMC.Recent technical improvements on mitigating the sign problem [63] enabled Ref. 64 to treat this model at the temperature of T = 0.025D with D the half bandwidth, which required large computational resources.This temperature is relatively high as the authors stated that in the study of a simpler two-band twosite model, where they reached T = 0.0125D, it was "intractable" to reach low enough temperatures to clarify whether bad metal/spin freezing behavior was intrinsic or not [63].
We study the model on the two-dimensional square lattice, i.e. using ε ab k = −2t (cos k x + cos k y ) δ ab .We use the momentum patching of Ref. 51; this definition is also used in the single-band computations of Figs. 1 and 3. We note that this model is not directly relevant to layered materials where the t 2g orbitals are relevant, because in the physical situation the two dimensionality will break the three-fold orbital degeneracy.However the system is well defined as a theoretical model and useful to demonstrate the power of our methods.As is the case for the CI method, the DMRG method used here is easily able to treat a large number of bath sites if the number of correlated sites is small: for L c = 1, DMRG has already often be proven to treat L b > 120 bath sites, and for L c = 2, L b > 80 is easily accessible [39,40].However for more correlated sites, the number of bath sites that can be added at given computational cost decreases.For L c = 6, we use L b = 18, i.e.L b /L c = 3, which we have shown to be sufficient to produce reliable results in the single-site case without requiring overly large computational resources (computation time of several hours per DMFT iteration on two cores).
We have tested the two-site calculation by converg-ing the DMFT loop for the three-band Hubbard model Eq. ( 15) with U and J = 0 and comparing the results with a corresponding two-site single-band DCA.Perfect agreement is obtained (not shown).Non-zero values of U and J create additional entanglement and make computations more costly.It is then a decisive question whether a real-space or a momentum-space representation of the impurity-cluster is less entangled.We discuss this in Appendix B finding that for the single-band Hubbard interaction, both representations yield similar entanglement, whereas for the Hubbard-Kanamori interaction, the realspace representation is much less entangled.Computational cost is therefore tremendously reduced by using the real-space representation, which comes with an offdiagonal hybridization function.This is the opposite behavior as observed for QMC, where the off-diagonal hybridization function creates a severe sign problem.We further note that in the real-space representation, strong interactions yield a less and less entangled impurity problem, as electrons become more and more localized.We now present results for the more physically relevant case U = U − 2J with J = U/4.For these parameters, at half filling the critical interaction for the MIT in the single-site DMFT approximation is U c 1.3D [65].Fig. 6(a) shows that our results are consistent with this estimate: the dashed lines depict the single-site (1s) results, showing a metallic solution (spectral function nonzero at ω = 0) for U = D, and an insulating solution (spectral function zero at ω = 0) for U = 2D.In the two-site (2s) DCA (solid lines), by contrast, the critical value U c for the MIT is lowered.Even at U = D the ω = 0 spectral function is zero (the small non-zero value in panel (a) is an effect of broadening, as seen in panel (b)).The different nature of the metallic and insulating solutions is also visible on the imaginary axis in the different nature of the decay of the imaginary-time Green's function.This is plotted in Fig. 6(c) for U = D; clearly, a power-law decay is observed for the metallic solution obtained in the single-site DMFT, whereas an exponential decay is obtained for the insulating solution obtained within the two-site DCA.

V. CONCLUSION
This paper introduces an imaginary-time MPS (itMPS) solver for DMFT and shows that it can treat complex models, not easily accessible with other methods, at modest computational cost.This development establishes DMRG as a flexible low-coast impurity solver for realistic problems, such as encountered in the study of strongly-correlated materials.The crucial advance stems from the fact that imaginary-time evolution does not create entanglement, and hence allows to compute Green's functions numerically exactly, provided a ground state calculation is feasible.
The method can be improved in many ways.In particular, different representations of the impurity problem exhibit different degrees of entanglement, so optimizing the representation of the impurity problem is a promising route.Ideas from ED approaches for constructing relevant subspaces [21][22][23][24] of the Hilbert space may lead to further improvements.Such techniques have been successfully combined with MPS [66].Another route to reduce computational effort and by that reach even more complex models could consist in performing computations for the reduced dynamics of the impurity [67], or making use of extremely cheap algorithms for computing the Green's function at elevated temperatures [50].Finally, we note that using MPS as an impurity solver makes accessing entanglement as a quantity for understanding the properties of the embedded impurity cluster very easily accessible.Proposals in this direction have been made for cellular DMFT [68] and for impurity models generally [69].The main challenge in solving the ground state problem of a typical cluster-bath Hamiltonian as encountered in DMFT, stems from the fact that DMRG is a variational procedure that is initialized with a random state, which is then optimized locally.A local optimization procedure is slow when optimizing a global energy landscape.In addition, the local optimization is prone to getting stuck in local minima, if no "perturbation steps that mix symmetry sectors" are applied.The standard perturbation techniques for single-site DMRG [70][71][72] rely on "perturbation terms" that are produced by contracting the Hamiltonian with the MPS.If the Hamiltonian itself does not contain terms that mix the symmetry sector, these methods do not work.
A typical cluster-bath Hamiltonian has both features, a global variation of the potential energy and parts that are not connected with symmetry mixing terms, such as in the three-band Hubbard Kanamori model at J = 0.This situation is sketched in Fig. 7.
In Ref. 39, the models under study allowed to solve this problem using the non-interacting solution.For the general models studied in the present paper, an unbiased numerical technique has to be employed.What we do in practice, is to first find the ground state of a system with additional symmetry-mixing couplings (denoted as red solid lines in Fig. 7) that are then adiabatically switched off.In practice, we sweep 5 to 10 times with additional hoppings of 10% magnitude of the physical hoppings, and another 5 to 10 times with additional hoppings of 1% magnitude.After these preliminary sweeps, the quantum number (e.g.particle number) distribution has globally converged, and we can continue with converging the ground state of the exact Hamiltonian.

Convergence of DMFT iteration
The calculations for the three-band single-site DMFT in Sec.IV A are only trivially parallelized using one core to compute the imaginary time evolution of each the particle (>) and the hole (<) Green's functions G ≷ a (τ ).In Fig. 8, we show the convergence DMFT loop for the single-site DMFT for the three-band Hubbard-Kanamori model as studied in Fig. 5. Fig. 8 DMRG is able to eliminate these superpositions as potential energy is separated locally, i.e. in the same basis in which DMRG optimizes the reduced density matrix in order to discard irrelevant contributions.In principle, as mentioned in Appendix C, ideas from basis selective approaches in exact diagonalization are a different method to account for the fact that many states in the Hilbert space have a negligible weight for the computation of the Green's function and only few physically relevant states occupy a small fraction of the Hilbert space.Among these are the truncated configuration interaction (CI) [20,23,24,56], the basis-selective ED [21] or the coupled cluster methods in quantum chemistry.As these methods can be combined with DMRG [66], they might be a further route to construct efficient representations of the impurity-cluster problem In the present paper, the question of the least entangled representation of the impurity problem shall be restricted to the question of which basis to choose in a DCA calculation.This is of high relevance also in another context: In the real-space representation, the hybridization function becomes off-diagonal.For CTQMC, this generates a sign problem.In our approach, this doesn't affect computational cost much in the single-band Hubbard model.It even leads to a tremendous reduction of computational cost for the three-band Hubbard Kanamori interaction.

DCA in momentum or real space
The complexity of the interaction determines whether the real-or the momentum-space representation of the cluster-bath Hamiltonian is less entangled.In real space, the interaction has a simple form, but the hybridization function has off-diagonal contributions, which result in additional couplings of cluster and bath sites.In momentum space, the hybridization function is diagonal but the interaction becomes off-diagonal.The additional couplings induced by that depend on the complexity of the interaction.
Let us be more concrete.For the two-site case, the discrete Fourier transform yields the even and odd superposition of the real-space cluster.
where the index of d † K labels momentum patches K and the index of d † i labels real space cluster sites i.There might be further indices labeling spin or orbital.
In real-space, the hybridization function has the form where the symmetry of the real-space cluster imposes Λ ij (z) = Λ ji (z).In momentum space, the hybridization function is diagonal and symmetry is reflected in the reduced number of bath sites per patch L b = L b /L c , where L c = 2 is the number of momentum patches.We choose to use the momentum representation for the bath discretization, as was done for the real-axis in Ref. 39.While on the real-frequency axis this is the only viable option, the bath fitting on the imaginaryfrequency axis via (12) is possible also for the off-diagonal real-space case.In real space, e.g. , particle hole symmetry can be easily imposed in the fitting procedure, while this is not possible in momentum space.
Given the parameters of the momentum space representation obtained by performing a bath fit via (12), we define the parameters of the equivalent real-space representation as follows: In momentum space, bath parameters are indexed by l K = 1, ..., L b , L b = L b /L c and in real space, bath parameters are indexed by l = 1, ..., L b , then From this one could naively expect that the real-space representation is less entangled if Numerical experiments show that the real-space representation is much more favorable than this estimate.For a single-band Hubbard model, we find about the same entanglement in the real space and the momentum space representation, with slight advantages for momentum space.In the three-band Hubbard Kanamori model, the real-space representation is considerably less entangled and leads to a tremendous reduction of computational cost.In particular, we were not able to obtain the results of Fig. 6 in the momentum space representation when using L b /L c = 3, only for L b /L c = 2 but then at much higher computational cost.
Appendix C: Green's functions from matrix product states Even though the following discussion is not needed to set up the imaginary-time real-time impurity solver, it is of general interest in this context and might stimulate further advancements.
A computation of A(ω) = ψ 0 |δ(ω−(H −E 0 ))|ψ 0 via a computation of all eigen states of H is extremely redundant as only a tiny neighborhood N = {|ψ ψ|H|ψ 0 = 0} of a the single-particle excitation |ψ 0 contributes in the sum (inserting identities n |E n E n |) in A(ω).In ED, this is exploited by systematically constructing the subspace N by spanning it using particle-hole excitations [20,21], which might also be a viable route for further developments within DMRG [66].In DMRG, one needs to make a statement about the entanglement of the states in the subspace N : one might note that these are in general more strongly entangled than the single-particle excitation |ψ 0 , but should still be much less entangled than the rest of the Hilbertspace.This is illustrated in the sketch Fig. 9.
In Ref. 30, some of us argued that expanding the spectral function in a family of orthogonal functions is a natural way to construct a basis for N , starting from the lowly entangled |ψ 0 and successively increasing entanglement of states and thereby computational complexity in a sequence of basis states |ψ n .Ref. 30 discussed the expansion of A(ω) in Chebyshev polynomials  [30].This is due to the observation that error accumulation in the Chebyshev recursion is worse conditioned than in time propagation [30], which necessitates to keep the error in a single step of the Chebyshev recursion much smaller than in the equivalent time evolution step, which in turn requires to use higher bond dimensions in the Chebyshev recursion making it less efficient.In addition to the statements of Ref. 30, we note here that the sequence produced by the Lanczos algorithm, can be associated with an expansion of the spectral function in polynomials that are orthogonal with respect to an inner product weighted by w(x) = A(x) [74].This is very efficient but numerically unstable.
In contrast to the previous methods, which generate an increasingly complex basis when determining the spectral function to a higher and higher precision, correctionvector DMRG aims to optimize a state in frequency space, which a priori contains contributions that have undergone an infinitely long time evolution.As time evolution creates entanglement, these states are much too strongly entangled for an efficient treatment.They are "far away" from the controlled, lowly entangled singleparticle excitation |ψ 0 .In order to still perform a meaningful computation in frequency space, one introduces a so-called (Lorentzian) broadening parameter η that damps out contributions from an infinite time evolution.One does then not obtain the exact spectral function but a broadened version as in Eq. (10).The broadening parameter has to be guessed a priori : If it is chosen too small, high entanglement prevents convergence of the calculation.If it is chosen too large, one will be far from the exact version of the spectral function.In the expansion methods discussed above, by contrast, one can stop the computation simply when it becomes too costly.If one has not recovered the exact A(ω) at this point, a broadened version can be systematically constructed with an a posteriori determined η as in Eq. (10).

FIG. 1 .
FIG. 1. (Color online) Panel (a): Imaginary time correlation functions G ≷ (τ ) defined in Eq. (4a) for an impurity model arising in the context of the two-site dynamical cluster approximation to the single-band Hubbard model on the square lattice with next-nearest neighbor hopping t /t = 0.3, halfband width D = 4t, interaction U = 2.5D and band filling n = 0.96 in the paramagnetic phase.See Ref. 51 for definition of the model and the meaning of the orbital (patch) quantum number "K = ±".The dashed line is obtained using linear prediction for times τ D ≥ 200.Panel (b): Maximal bond dimension m of time evolved states.The MPS computation used a Hamiltonian representation of the impurity model with Lc = 2 correlated sites and L b = 14 bath sites.The hybridization function of the impurity model Λ discr was fitted using Eq.(12) for β eff = 315/D and α = 0.For the ground state optimization, we enforced a maximal bond dimension of m = 300.The Krylov time evolution used a time step of ∆t = 0.1/D and allowed for a maximal global truncation error of 10 −4 at each time step adjusting bond dimensions automatically.This leads to an immediate decay of m at τ 0 from m = 300 down to m 110, as seen in panel (b).We used the global SU(2) symmetry of the Hamiltonian to reach these low values of the bond dimension.
FIG.2.(Color online) Fit of the hybridization function in the two-site DCA problem studied in Figs.1 and 3, but here for the case U = 0.The minimization(12) was done using α = 0 and a frequency grid defined by β eff = 100/D and a cutoff frequency of ωc = 6D.Evidently, the quality of the fit does not improve any more for L b /Lc 9.

FIG. 3 .
FIG. 3. (Color online)Real-and imaginary-frequency Green's functions computed by converging the DMFT self-consistency equation Eq. (11) for the two-site dynamical cluster approximation to the single-band Hubbard model on the square lattice with next-nearest neighbor hopping t /t = 0.3, half-band width D = 4t, interaction U = 2.5D, and band filling n = 0.96 in the paramagnetic phase (as in Fig.1).See Ref.51 for definition of the model and the meaning of the orbital (patch) quantum number "K = ±".Panels (a-c): Electron spectral function A(ω) = − 1 π ImG ret (ω) obtained by converging on imaginary-frequency axis (itMPS) using number of bath sites and broadenings as specified in the figure, and compared to unbroadened (η = 0) real-frequency axis computation using L b /Lc = 39 bath sites per correlated site of Ref. 39. Panel (d): Converged Matsubara Green's function for number of bath sites shown, compared to numerically exact quantum Monte-Carlo result of Ref. 51, computed at β = 200/D.

FIG. 4 .
FIG. 4. (Color online) Density per orbital as function of chemical potential for three-band Hubbard-Kanamori model Eq.(13) using the semi-elliptic density of states Eq. (14) and U = 12t, obtained from single-site DMFT approximation evaluated using imaginary MPS (crosses) and CTQMC data (circles, Figure1of Ref. 61, inverse temperature β = 50/t).In the DMRG computations the bath fitting was performed using β eff = 100/t, ωc = 6t and α = 1 with three bath sites per correlated site (L b /Lc = 3).The maximal matrix dimensions was m = 300 for the ground state calculation, exploiting the SU(2) symmetry, which leads to the high precision (H − E)2  10 −14 .For the time evolution, we computed G FIG. 5. (Color online) Imaginary part of Matsubara axis self energy Σ (panel (a)) and imaginary part of hybridization function Λ (panel (b)) for densities shown obtained from converged itMPS solution of single-site DMFT for three-band Hubbard-Kanamori model (15) for U = 8t and J = U/6.Crosses (color online) represent itMPS data and black circles depict CTQMC data from Figure3of Ref. 62, computed at inverse temperature β = 100/t.We choose all parameters as described in the caption of Fig.4, in particular, for the bath fitting(12), we use β eff = 100/t, ωc = 6t and α = 1.Choosing α = 1 enforces agreement for low frequencies at the price of disagreement at high frequencies, which is observed in both panels (a) and (b).In panel (b), Λ denotes the hybridization function that is fitted with the hybridization function Λ discr of the discrete impurity model.

FIG. 6 .
FIG. 6. (Color online) Comparison of results obtained using MPS methods with for L b /Lc = 3 for single site (1s) and twosite (2s) DMFT approximations to the Hubbard-Kanamori model (13) on the two-dimensional square lattice with half bandwidth D = 4t, ε ab k = −2t (cos kx + cos ky) δ ab , U = U − 2J, J = U/4 and n = 3 (µ = 5U/2 − 5J), that is, in the particle-hole symmetric case.Panel (a): spectral functions for broadening η = 0.2D; panel (b) broadening η = 0.05D.In panel (c), we show the imaginary-time evolution of G > (τ ) as defined in (4a), confirming by comparison to a calculation for a smaller bath L b /Lc = 2 that this quantity converged with respect to size.The maximal bond dimension for the ground state search was m = 1000.
VI. ACKNOWLEDGEMENTS FAW thanks G. K.-L.Chan for stressing the relevance of converging the DMFT loop on the imaginaryfrequency axis, and N.-O.Linden for helpful discussions.AJM and US acknowledge the hospitality of the Aspen Center for Physics NSF Grant 1066293 during the inception of this work.FAW and US acknowledge funding by FOR1807 of the DFG.AJM and AG were supported by the US Department of Energy under grant ER-046169.

Appendix A: Further technical details 1 .
Ground state optimization

FIG. 7 .
FIG. 7. (Color online) Sketch of a typical cluster-bath Hamiltonian (Lc = 3, L b = 6) when it's mapped to a onedimensional chain.Dashed lines depict couplings that do not mix symmetry sectors, and solid lines depict couplings that mix symmetry sectors.

FIG. 8 .
FIG. 8. (Color online) Single-site DMFT for three-band Hubbard-Kanamori model as studied in Fig. 5.Here for the case n = 1.77 (µ = 5.0) and L b /Lc = 3.To obtain the solution for n = 1.79 as shown in Fig. 5, we chose µ = 5.1 and started from the n = 1.77 solution.Panel (a): Convergence of Matsubara Green's function in the DMFT loop, starting from the non-interacting solution.Panel (b): Convergence of the density in the DMFT loop.Panel (c): Convergence of the ground state energy per particle in the DMFT loop.Panel (d): Computation time.An iteration on the Matsubara axis takes about 30 min.The final real-axis computation (iteration 31) is considerably more expensive, but can still be optimized.
l2=l−L b for l = L b + 1, ..., L b .(B5) Whereas the momentum-space Hamiltonian has L b non-zero couplings V Kl K , the real space Hamiltonian has L c × L b couplings V il .On the other hand, the interaction part generates L c ×(L c −1) additional non-local couplings in the momentum-space representation as compared to the real-space Hamiltonian.

FIG. 9 .
FIG. 9. (Color online) Single particle excitation |ψ0 of the ground state |E0 and the subspace N = {|ψ | ψ|H|ψ0 = 0} of the Hilbert space that is relevant for the computation of a single-particle spectral function of the form ψ0|δ(ω − H)|ψ0 .The single particle excitation is very lowly entangled, the subspace is more strongly entangled, but still in general more lowly entangled than the rest of the Hilbert space.