Fork Tensor Product States - Efficient Three Orbital Real Time DMFT Solver

We present a tensor network especially suited for multi-orbital Anderson impurity models and as an impurity solver for multi-orbital dynamical mean-field theory (DMFT). The solver works directly on the real-frequency axis and yields very high spectral resolution at all frequencies. We use a large number $\left(\mathcal{O}(100)\right)$ of bath sites, and therefore achieve an accurate representation of the bath. The solver can treat full rotationally invariant interactions with reasonable numerical effort. We show the efficiency and accuracy of the method by a benchmark for the testbed material SrVO$_3$. There we observe multiplet structures in the high-energy spectrum which are almost impossible to resolve by other multi-orbital methods. The resulting structure of the Hubbard bands can be described as a broadened atomic spectrum with rescaled interaction parameters. Additional features emerge when $U$ is increased. The impurity solver offers a new route to the calculation of precise real-frequency spectral functions of correlated materials.


I. INTRODUCTION
Strongly correlated systems are among the most fascinating objects solid-state physics has to offer. The interactions between constituents of such systems lead to emergent phenomena that cannot be deduced from the properties of non-interacting particles [1].
One of the most widely used methods to describe strongly-correlated electrons is the dynamical mean-field theory (DMFT) [2,3]. DMFT treats local electronic correlations by a self-consistent mapping of the lattice problem onto an effective Anderson impurity model (AIM). Calculating the single particle spectral function of this impurity model in an accurate and efficient way is at the heart of every DMFT calculation. To this end, many numerical methods have been developed or adapted. These are based for instance on continuoustime quantum Monte Carlo (CTQMC) [4,5], exact diagonalization (ED) [6][7][8], the numerical renormalization group (NRG) [9,10], configuration interaction (CI) based solver [11,12], and also the density-matrix renormalization group (DMRG) with matrix-product states (MPS) [13,14].
Every algorithm has strengths and weaknesses: CTQMC is exact apart from statistical errors on the imaginary axis and can deal with multiple orbitals, but it is in some cases plagued by the fermionic sign problem. Additionally, an ill-posed analytic continuation is necessary to obtain real-frequency spectra, which therefore become broadened, especially at high energies. ED directly provides spectra on the real axis, but it is severely limited in the size of the Hilbert space, i.e. in the number of bath sites. Quite recently, NRG was shown to be a * daniel.bauernfeind@tugraz.at viable three-band solver by exploiting non-abelian quantum number conservation [15][16][17]. NRG works on the real axis and captures the low-energy physics well, but it has by construction a poor resolution at higher energies. Another interesting route that has been proposed recently are solvers that tackle the problem of exponential growth of the Hilbert space using ideas from quantum chemistry, i.e. the configuration interaction [11,12]. They allow to go beyond the small bath sizes of ED, keeping all the advantages such as absence of fermionic sign problems. However, in multi-orbital applications (see Appendix of Ref. [12]), the spectral resolution has so far been restricted by the restricted number of bath sites (O (20)).
MPS based techniques like DMRG, finally, do not suffer from a sign problem and can be used on the real-as well as on the imaginary-frequency axis. The price to pay for the absence of the sign problem is an, in general, very large growth of bond dimension with the number of orbitals.
Dynamical properties and spectral functions can be calculated within DMRG and have been used for impurity solvers, e.g. with the Lanczos-like continuedfraction expansion [18,19]. Other solvers using the more stable correction vector [20] and dynamical DMRG (DDMRG) [21] methods were developed [22][23][24][25]. Both algorithms produce very accurate spectral functions, but have the disadvantage that a separate calculation for each frequency has to be performed. The Chebyshev expansion [26] with MPS [27], supplemented by linear prediction [28], was used for impurity solvers in the single band case [29] and for two bands [30]. Recently, some of us introduced a method based on real-time evolution [31][32][33][34] and achieved a self consistent DMFT solution for a twoband model [35]. In such calculations, the physical orbitals for each spin direction are usually combined to one large site in the MPS. Three or more orbitals have not been feasible with this approach, because of a large increase in computational cost with the number of orbitals. Another MPS-based solver, which works on the imaginary axis, was recently introduced [36] and it was applied as a solver for three bands in two-site cluster DMFT. It was supplemented by a single real-time evolution to compute the spectral function, avoiding the analytic continuation. However, this method is restricted by the number of bath sites which can be employed. In the calculation mentioned, only three bath sites per orbital were used, limiting the energy resolution for real-frequency spectral functions.
In the present paper, we introduce a novel impurity solver which works directly on the real-frequency axis.
To this end, we use a tensor network that captures the geometry of the interactions in the Anderson model better than a standard MPS. Our approach is to some extent inspired by the work of Ref. [37], which used a similar network for a two orbital NRG ground state calculation. We are not restricted to a small number of bath sites. This is imperative for exploiting the spectral resolution achievable with real-time calculations. We emphasize that (i) our method is by construction free of any fermionic sign problem, (ii) one can fully converge the DMFT self-consistency loop on the real-frequency axis and (iii) we can achieve an almost exact representation of the bath spectral function. We apply this method to multi-orbital DMFT for the testbed material SrVO 3 and show that one can resolve a multiplet structure in the Hubbard bands, keeping at the same time a good description of the low-energy quasi-particle excitations.
The paper is structured as follows. First we show how impurity solvers with tensor networks work in general and introduce our new tensor network approach which we call fork tensor-product states (FTPS) (Sec. II). Next we explain in detail how our solver is used in the context of multi-orbital DMFT (Sec. III). In Sec. IV, we apply our approach to SrVO 3 and discuss the multiplet structure that the FTPS solver allows to resolve. In order to check the accuracy of the method, we also compare the FTPS results to CTQMC for SrVO 3 . Finally, we show the efficiency of the FTPS solver by applying it to a five-orbital model.

II. TENSOR NETWORK IMPURITY SOLVERS
The Anderson impurity model (AIM) describes an impurity (with Hamiltonian H loc ) coupled to a bath of noninteracting fermions hybridized with it. A typical AIM Hamiltonian is given by: where c † mlσ (c mlσ ) creates (annihilates) an electron in band m (m ∈ {1, 2, 3} for a three-orbital model) with spin σ at the l-th site of the system (the impurity has index l = 0, the bath degrees of freedom have l ≥ 1), and n mlσ are the corresponding particle number operators. H DD describes density-density (DD) interactions between all orbitals and H SF-PH are the spin-flip and pair-hopping terms. This three-orbital Hamiltonian is not only important in the context of real-material calculations. It has also been studied extensively on the model level, most importantly because it hosts unconventional correlation phenomena. For a selection of recent work, see for instance Refs. [15,[38][39][40][41].
An impurity solver calculates the retarded impurity Greens function G(t) of the interacting problem (1), either in real or imaginary time t. In the present paper, we introduce a new tensor network similar to an MPS, which can be used as a real-time impurity solver for three orbitals. We first introduce MPS before moving on to what we call fork tensor-product states (FTPS) in Sec. II B.
Each site i has a local Hilbert space of dimension d i spanned by the states |s i . Through repeated use of singular-value decompositions (SVDs), it is possible to factorize every coefficient c s1,··· ,s N into a product of matrices [14], i.e. into an MPS, Each A si i is a rank-3 tensor, except the first and last ones (A s1 1 , A s N N ), which are of rank two. The index s i is called physical index, and the matrix indices, which are summed over, are the so called bond indices. A general state of the full Hilbert space is unfeasible to store, but it can be shown that ground states are well described by an MPS with limited bond dimension m (dimension of the bond index) [42].
In complete analogy to the states, one can factorize an operator into what is called a matrix-product operator (MPO) [14], where each W si,s i i is a rank-4 tensor. Tensor networks in general have a very useful graphical representation, which is shown for an MPS and an MPO in Fig. 1. Note that when we use the term MPS we always mean a onedimensional chain of tensors as shown in Fig. 1. To calculate Greens functions within the MPS formalism, one usually first applies the DMRG [13,14], which acts on the space of MPS and finds a variational ground state |ψ 0 and ground-state energy E 0 . It minimizes the expectation value Graphical representation of a fork tensor product state (FTPS) for multi-orbital AIM. The idea to separate bath degrees of freedom leads to a fork-like structure. In this picture, a two-orbital model with four bath sites each is shown. Orbitals are labeled A and B and the arrows denote the spin. Each spin-orbital combination has its own bath sticking out to the right. As in Fig. 1, the vertical lines are the physical degrees of freedom (all of dimension two, for empty, resp. occupied bath sites). All other lines are bond indices and like in the MPS they are summed over. As mentioned in the text, the bath is represented in star geometry due to the smaller bond dimensions needed. The bath sites are ordered according to their on-site energies. Two example hoppings V1 and V2 are drawn.
by updating usually two neighboring MPS tensors before moving on to the next bond. This procedure also yields the Schmidt decomposition of the state at the current bond on the fly. The DMRG approximation is to keep only those states with the largest Schmidt coefficient. It is important to note that one can perform a DMRG calculation for any tensor network, as long as one can generate a Schmidt decomposition [43]. For obtaining the Greens function, we employ an evolution in real time. Eq. (2) is split into the greater G > and lesser Greens function G < : which are calculated in two separate time evolutions. This is done by first applying c † (or c) and then time evolving this state and calculating the overlap with the state at time t = 0. The time evolution is the most computationally expensive part, since time evolved states are not ground states anymore, and the needed bond dimensions usually grow very fast with time.

B. Fork Tensor Product States (FTPS)
So far, the usual way of dealing with Hamiltonians like Eq. (1) using MPS [29,30,35] has been to place the impurity in the middle of the system and the up-and down-spin degrees of freedom to its left and right, resp. The local state space of each bath site then consists of M spinless-fermion degrees of freedom, with dimension 2 M , where M is the number of orbitals in the Hamiltonian Eq. (1). This exponential growth is usually accompanied by a very fast growth in bond dimension when using the above arrangement. We did indeed encounter this very fast growth upon calculating the ground state of some one-two-and three-orbital test cases.
For treatment by MPS, the general Hamiltonian Eq. (1) with hopping terms from the impurity to each bath site is usually transformed into a Wilson chain with nearest-neighbor hoppings only, i.e. of the form [10]. This was thought to be necessary since long-range interactions look problematic for MPSbased algorithms. Quite recently, though, it was discovered that MPS can deal with the original form of H bath in Eq. (1) better [44]. Because all hopping terms in H bath originate from the impurity, this is called the star geometry. The reason for the better performance is that in the star geometry one has many nearly fully occupied (empty) bath sites with very low (high) on-site energies l .
Since basis states with many unoccupied low-energy sites have a very low Schmidt coefficient, these states are discarded from the MPS. The same holds for occupied high energy sites. However, when dealing with multiorbital models, the star geometry is not enough to be able to calculate Greens functions using MPS. The growth of the bond dimensions still makes those calculations unfeasible.
The key idea of the present work is to construct a tensor network which is beyond a standard MPS, but similar enough to be able to use established methods like DMRG and time evolution. From Hamiltonian (1) one can immediately notice that there are no terms coupling bath sites of different orbitals. Hence, it might not be advantageous to combine those, not directly interacting, degrees of freedom into one large physical index in the MPS.
Our proposed tensor network, therefore, separates the bath degrees of freedom as much as possible. It consists of separate tensors for every orbital-spin combination, each connected to bath tensors as shown in Fig. 2. This tensor network is no MPS anymore, since there are some tensors (labeled A ↓ and B ↑ in the example of Fig. 2) that have three bond indices and one physical index, i.e. which are of rank 4. Cutting any bond splits the network into two separate parts. Therefore, one can calculate the Schmidt decomposition in a way very similar to an MPS, which means that also DMRG is possible. The main bottleneck of calculations with FTPS is to perform SVDs of the rank-4 tensors representing the impurities. When all bond indices have the same dimension χ, it is necessary to do a SVD for a χ 2 d × χ matrix with computational complexity O(χ 4 d). However, as we show below, this operation does not pose a substantial problem for calculations using FTPS. Since the impurity tensors pose the biggest challenge, our tensor network would likely also allow us to deal with the chain geometry without a drastic increase in computational cost. In the present paper we will only use FTPS with baths in star geometry. The proposed FTPS are similar to the tensor network used by Holzner et al. [37] to perform NRG calculations for ground state properties of an AIM with two orbitals.
The three-legged tensors in our network ( Fig. 2) can also be interpreted as two coupled junctions with three legs in the language of Ref. [45], where it has been shown that DMRG is possible on such junctions. Furthermore, our approach has similarities with the so called Tree Tensor Networks (TTN) [43,[46][47][48].

Time Evolution
Time evolution with the Hamiltonian Eq. (1) is not straightforward, since it features long-range hoppings. Possible methods include Krylov approaches [49], the time-dependent variational principle [50,51] and the series expansion of e iHt proposed by Zaletel et al. [52]. In this work, however, we use a much simpler approach.
First, we split the Hamiltonian into the following terms: (i) the spin-flip and pair-hopping terms (1)), (ii) the density-density interaction terms H DD , and (iii) all other terms H free = H bath + 0 mσ n m0σ . With these terms, we write the timeevolution operator for a small time step ∆t using a second-order Suzuki-Trotter decomposition [53], Note that in this decomposition, the order of the spinflip and pair-hopping terms is important. The order of operators in the second product must be opposite to the one in the first. We see that Eq. (8) involves three different operators H SF-PH m,m , H DD and H free , each of which will be treated differently.
Time evolution of the density-density interactions is performed with an MPO-like representation of the timeevolution operator e −i ∆t 2 HDD . For a three-orbital model, first the full matrix (4 3 × 4 3 ) of e −i ∆t 2 H DD is created, which is then decomposed into MPO-form by repeated SVDs. Since H DD only consists of density-density interactions, no fermionic sign appears in e −i∆tHDD .
Time evolution of the spin-flip and pair-hopping terms is more involved than the density-density interactions, since the operators change the particle numbers on the impurity sites. Therefore, it can be difficult to deal with the fermionic sign of the time evolution operator when the impurities are not next to each other in the fermionic order. It turns out that the spin-flip and the pair-hopping terms have the propertyÂ 3 =Â individually, withÂ being either the spin-flip or the pair-hopping operator, resp. Furthermore they commute with each other allowing us to separate them without Trotter error. The time-evolution operator of JÂ is then given by: For this operator, an MPO can be found for which the fermionic sign can easily be determined [54].
To time evolve the bath terms we use an iterative second-order Suzuki-Trotter breakup for each term in H free . Neglecting orbital (m) and spin (σ) indices, the first step in this breakup is the following: . Next we split off H 2 and iterate this procedure until we end up with with N b the number of bath sites and H mlσ = l n mlσ + V l c † m0σ c mlσ + h.c. . In the above equation, we neglected the term 0 n m0σ that we add to H m1σ . Eq. (10) is a product of two-site gates (an operator acting nontrivially only on two sites) with one of the two sites always being the impurity. This means that those two sites are not nearest neighbors in the tensor network. To overcome this problem, we use so called swap gates [14,55]. The two-site operator swaps the state of site i (s i with occupation n i ) with the state of site j (s j with occupation n j ). The factor (−1) ninj gives the correct fermionic sign and is negative if an odd number of particles on site i gets swapped with an odd number of particles on site j. To be more precise, the matrix representation of the swap gates used in this work is: It turns out that every swap gate can be combined with an actual time evolution gate without additional computational time. For example, the first step in this time evolution would be to apply e −i ∆t 2 Hm1σ . Immediately afterwards, even before the SVD (to separate the tensors again), the swap gate is applied so that the impurity and the first bath sites are swapped. By repeating this process one moves the impurity along its horizontal arm in Fig. 2. Because a second-order decomposition is used, now all time evolution gates except the one at site N b have to be applied again. But now, the impurity and bath site needs to be swapped before time evolution.
Note that the algorithm presented above cannot only be used to perform real-time evolutions, but it is applicable also to evolution in imaginary time simply by replacing idt by dτ .

III. MULTI ORBITAL DMFT WITH FTPS
In this section we present details of our impurity solver. We refer to Refs. [2,56] for DMFT in general, and to Refs. [57,58] for DMFT in the context of realistic abinitio calculations for correlated materials.
In the latter approach, called density-functional theory (DFT)+DMFT, the correlated subspace is described by a Hubbard-like Hamiltonian. Within DMFT, this model is mapped onto the AIM Hamiltonian (1). This mapping defines the bath hybridization function ∆(ω) describing the influence of the surrounding electrons.
Since FTPS provide the Greens function of the AIM on the real-frequency axis, also the self-consistency loop is performed directly for real frequencies. For calculating the bath hybridization, we use retarded Greens functions with a finite broadening η SC in order to avoid numerical difficulties with the poles of the Greens function. Throughout this work, we use η SC = 0.005 eV [59]. The impurity solver calculates the self energy Σ(ω) of the AIM, given the bath hybridization function ∆(ω) and the interaction Hamiltonian on the impurity. To this end, our solver performs the following steps, which are explained in more detail in the text below: 1. Obtain bath parameters l and V l by a deterministic approach based on integration of the bath hybridization function ∆(ω).
2. Calculate the ground state |ψ 0 and ground-state energy E 0 of the interacting problem.
3. Apply impurity creation or annihilation operators, and time evolve these states to determine the interacting Greens function ( Eq. (2) ).

Fourier transform Eq.
(2) to obtain G(ω) and calculate the local self-energy, To perform step 1 we use similar to Refs. [10] (NRG) and [44]. Each interval I l corresponds to a bath site. This discretization can be interpreted as representing each interval I l as a delta peak at position l and weight V 2 l . Sum rules for such discretization parameters can be found analytically [60]. In this work we choose the length of each interval such that the area of the bath spectral function − 1 π Im ∆(ω) is approximately constant for each interval [61]. For the case at hand, this discretization was found to be numerically more stable than using intervals of constant length. Unless stated otherwise, we use N b = 109 bath sites per orbital and spin. We note that this scheme is not restricted to diagonal hybridizations. In the general case of off-diagonal hybridizations the hybridization function is a matrix ∆. Therefore, instead of taking the imaginary part we can use the bath spectral function i 2π (∆ − ∆ † ). Similarly to Eq. (14), we represent each interval by one delta-peak for each orbital. For instance, fixing l to the center of the interval, the hopping parameters V l can be found systematically from the Cholesky factorization of Most importantly, this scheme does not involve any fitting procedure on the Matsubara axis. A very similar approach was developed independently in Ref. [62].
In step 2 we use a DMRG approach with the following parameters, unless specified otherwise. The truncated weight t w (sum of all discarded singular values of each SVD) is kept smaller than 10 −8 . When spin-flip and pair-hopping terms are neglected, we use an even smaller cutoff of 10 −9 . Note that, except in the five-band calculation, we do not restrict the bond dimensions by some hard cutoff (see Appendix 2).
During time evolution (step 3), we use a truncated weight of t w = 2 · 10 −8 , or 10 −8 with density-density interactions only. We time evolve to t = 16 eV −1 , with a time step of ∆t = 0.01 eV −1 . Greens functions are measured every fifth time step. The time-evolution operator of H loc is applied using the zip-up algorithm [55]. Afterwards the Greens functions are extrapolated in time using the linear prediction method [28,35] up to t = 250 eV −1 . Time evolution is split into two runs one forward and one backwards in time [63] to be able to reach longer times.
In the Fourier transform to ω-space (step 4), we use a broadening in the kernel e iωt−η F T |t| of η F T = 0.02 eV to avoid cutoff effects remaining after the linear prediction. The influence of the linear prediction on our results is discussed in Appendix 1. We want to stress that although a calculation with full rotational symmetry is more demanding, the computational effort is still very feasible. With the parameters mentioned above one full DMFTcycle takes about five hours on 16 cores.
To verify that our implementation of DMRG and time evolution produces correct results when used with our tensor network, we first compared Greens functions and ground-state energies for U = J = 0 for several bath parameter sets. The next step of our testing was to include density-density interactions, one term at a time. For example, we only included (U − J)n 10↑ n 30↑ and compared energy and Greens function to a standard one-orbital MPS solver. Finally, we also compared our method to the MPS two-band solver used in Ref. [35]. Indeed all tests performed produced correct energies and Greens functions.

IV. RESULTS
We performed DMFT calculations based on a band structure obtained from density functional theory (DFT) for the prototypical compound SrVO 3 , using the approximation of the Kanamori Hamiltonian (Eq. (1)). It has a cubic crystal structure with a nominal filling of one electron in the V-3d shell [64]. Due to the crystal symmetry, the five orbitals of the V-3d shell split into two e g and three t 2g orbitals. The latter form the correlated subspace. We performed the DFT calculation with Wien2k [65], and used 34220 k-points in the irreducible Brillouin zone in order to reach an energy resolution comparable with the η SC = 0.005 eV broadening.
The TRIQS/DFTTools package (v1.4) [66][67][68], which is based on the TRIQS library (v1.4) [69], was used to generate the projective Wannier functions and to perform the DMFT self-consistency cycle. Fig. 3 shows the main results of this paper, the DMFT spectral function A(ω) for SrVO 3 , (i) in the approximation of density-density interactions only and (ii) with full rotational invariance including spin-flip and pair-hopping terms. Overall, both cases show the well known features of the SrVO 3 spectral function [70,71]. We see a hole excitation at around −2 eV, and the quasi-particle peak at zero energy whose shape and position does not depend on the inclusion of full rotational invariance. In the upper Hubbard band, a distinctive three-peak structure can be seen. This structure has not been resolved in other exact methods like CTQMC (problem with analytic continuation, see below) or NRG (logarithmic discretization problem). In our real time approach, high energies correspond to short times, where the calculations are particularly precise [72]. Most methods allow to resolve structures in the Hubbard bands only in special cases (see Ref. [73] for an example using ED). Of course, atomic-limit based algorithms such as the Hubbard-I approximation or non-crossing approximation (NCA) show atomic-like features, but they have very limited accuracy for the description of the low-energy quasi-particle excitations in the metallic phase [74]. Thus, our FTPS solver combines the best of the two worlds, with atomic multiplets at high energy and excellent low-energy resolution at the same time.
The energies of the three peaks in the upper Hubbard band differ depending on whether SF-PP terms are taken into account or not. Details of this peak structure, as well as additional excitations visible at higher energies, will be discussed below in Sec. IV C.
First we examine the convergence of our results with respect to the number of bath sites and compare our spectrum to CTQMC. The following discussion is mostly based on calculations without spin-flip and pair-hopping terms. In this case, the calculations can be done faster and with higher precision, since there is no particle exchange between impurities. In all subsequent plots, we show results from calculations with DD interactions only.

A. Effect of Bath Size
In order to achieve a reliable high resolution spectrum on the real-frequency axis, it is imperative to have a good representation of the hybridization function ∆(ω) in terms of the bath parameters, for which a sufficient number of bath sites is needed. Fig. 4 shows how well a hybridization function can be represented with our approach (Eq. (14)) using a certain number of bath sites. We see that for N b = 9 bath sites (we always denote sites per orbital), ∆(ω) can be reconstructed only very roughly, which in turn gives an incorrect spectral function (Fig. 4 bottom). To some extent, the difference in the spectrum is due to the shorter time evolution and therefore a higher broadening η F T we were forced to use. For such a small bath, the finite size effects from reflections at the bath ends appear much earlier in the time evolution.
Increasing the number of bath sites to N b = 29, we observe that the reconstructed bath spectral function already shows the relevant features of ∆(ω). The spectrum is well converged for the largest bath sizes N b = 59 and N b = 109.

B. Comparison to CTQMC
In Fig. 5 we compare the converged spectral function of our approach (FTPS) with a spectrum obtained from CTQMC and analytic continuation. In both calculations, we used the same interaction Hamiltonian with densitydensity interactions only. The CTQMC calculation was performed with the TRIQS CTHYB-solver (v1.4) [75,76] with 3.2 · 10 7 measurements and at inverse temperature β = 200 eV −1 . For the analytic continuation we applied the ΩMaxEnt method [77].
The three-peak structure in the upper Hubbard band is not present in the CTQMC spectrum. We will show below in an example that even for a Greens function that does contain these peaks the analytic continuation does not resolve this structure. For another comparison, we consider the imaginary time Greens functions G(τ ) in Fig. 6. Apart from the effect of statistical errors, CTQMC provides an exact self consistent solution of DMFT on the imaginary-frequency axis. As mentioned above, when we use the FTPS solver, we formulate the DMFT self-consistency equations on the real-axis. To obtain an approximate finite temperature imaginary-time Greens function from FTPS that we can compare to the CTQMC result, we need to take the finite temperature of the CTQMC calculation into account. Therefore, we use the FTPS spectrum A(ω) and assume that we would obtain the same spectrum for a finite (but high enough) inverse temperature β, and use: The results in Fig. 6 show very good agreement on a logarithmic scale. Another important indication of the validity of our results is the value of A(ω = 0). To get a comparable number, we use the CTQMC imaginary time Greens function G(τ ) and Fourier transform it to get G(iω n ): Looking at the last few DMFT-cycles, we estimate it to be around A(ω = 0) = − 1 π lim iωn→0 G(iω n ) ≈ 0.272 eV −1 with fluctuations in the last digit. For the FTPS, the exact height of A(ω = 0) of the FTPS spectrum changes a little for each DMFT iteration, mainly due to slight variations in the linear prediction. Using the same prescription as for CTQMC, we estimate it to be A(ω = 0) = 0.28 eV −1 , with fluctuations of about 0.01 eV −1 . This agreement is very good considering that linear prediction has its strongest influence at small energies. Further benchmarks concerning the linear prediction can be found in Appendix 1. For better visibility of the τ > 0 data, the value of 9 · 10 −3 at τ = 0 is not shown.
Finally, we show that the ill-posedness of the analytic continuation is the most likely explanation for the missing peak structure in the upper Hubbard band of the spectral function obtained from the CTQMC data. To do so, we take the FTPS spectrum A(ω), calculate G(τ ) as described above, and perform the same analytic continuation that we did for the G(τ ) from CTQMC. We added noise of the order of the CTQMC error to the FTPS data. The resulting spectrum is shown in Fig. 7, and indeed the peak structure in the upper Hubbard band vanishes.

C. Discussion of Peak Structure -Effective Atomic Physics
In order to understand the peak structure observed in the spectral functions, we take a look at the underlying atomic problem, where for simplicity we start with density-density interactions only. We will show that the same arguments hold for full rotationally invariant interactions.
Tab. I shows the relevant atomic states and their cor- responding energies. The atomic model has a hole excitation at energy − 0 and three single electron excitations with energies U + 0 , U − 2J + 0 and U − 3J + 0 relative to the ground state. If we measure the energy differences between the three peaks of the upper Hubbard band in our results, we find values of 1.27 eV and 0.69 eV, which is close to the atomic energy differences of 1.2 eV and 0.6 eV (J = 0.6 eV). We also find the hole excitation at −2.0 eV. This indicates that we can describe the positions of the observed peaks approximately by atomic physics with effective parameters¯ 0 ,Ū andJ and widened peaks. Furthermore, the heights of the peaks roughly correspond to the degeneracy of the states in the atomic model (see Tab. I).
We can determineŪ = 5.97 eV (where U = 4.00 eV) from the energy difference of the peak highest in energy to the hole excitation. This increase ofŪ compared to U is plausible considering the following. When coupling the impurity to the bath, particles have the possibility to avoid each other by jumping into unoccupied sites of the bath. This results in a decrease of n ↑ n ↓ . To model this situation using atomic physics, one needs to increase the interaction strength. Finally, it is well known that J is much less affected by the surrounding electrons than U , since the latter is screened significantly stronger [78].
Tab. II shows how bare atomic parameters change when adding a bath and we see that our qualitative arguments give a correct idea of how parameters are rescaled.
Further evidence that the observed three-peaked structure is indeed a result of atomic physics can be seen in Fig. 8. It shows a closeup of the upper Hubbard band for three different values of J. The corresponding effective parametersJ are shown in Tab. II. We observe that also J is rescaled slightly, but the rescaling gets smaller for higher J. Furthermore, increasing J also increases the total width of the Hubbard band, which scales mostly linearly with J. At the same time, measuring the quasi- particle spectral weight as a function of J at constant U shows that it increases with increasing J, implying also an increasing critical U c for the metal-to-insulator transition [40]. Upon a careful inspection of the spectral function in Fig. 5, we observe small peaks at energies around 8 eV. A closeup of this energy region for different values of J is shown in Fig. 8. The energy difference between the peaks is close to 2J and can, again, be well explained by atomic physics, namely excitations into states with 3 electrons on the impurity (Tab. I) [79]. These excitations originate from small admixtures of N = 2 states to the ground state.
With atomic physics in mind, let us take a look again at the spectrum using full rotational symmetry (Fig. 3). The spin-flip and pair-hopping terms only contribute if there are two or more particles present. Thus, the quasi-  0.70 0.72 (2) particle peak and the hole excitation do not change. The atomic N = 2 sector does change, however. Diagonalizing the Hamiltonian, we find eigenstates with three different energies and differences of 3J = 1.8 eV and 2J = 1.2 eV, resp. Measuring the energy differences in Fig. 3, we find 3J = 1.75 eV and 2J = 1.32 eV. Estimat-ingŪ = 5.81 (5) we see that it does not change much compared to DD only [80]. Again, we can describe the spectrum approximately by atomic physics with effective parameters. Like in the case with density-density terms only, we also see the tiny excitations to states belonging to the atomic N = 3 sector.

D. Beyond Atomic Physics
The previous section showed that at U = 4.0 eV the spectral features in the Hubbard bands can be well described by atomic physics with effective parameters and widened peaks. It is not clear whether this picture is valid for higher interaction strengths U in the metallic regime. In Fig. 9 we show results with U = 5.5 eV at constant J = 0.6 eV. We see a shift of the upper Hubbard band to higher energies, but little shift of the holeexcitation. Also the central peak is shifted and gets slimmer since more weight is transferred into the Hubbard bands. Most importantly, as we approach the stronglycorrelated metallic regime, we clearly leave the realm where atomic physics can describe all the spectral features.
We find that the three-peak structure of the upper With the help of these lines one can discern a three-peaked structure again, but extended by a feature at the inside of the Hubbard band.
Hubbard bands smears out, and even vanishes. The closeup of the upper Hubbard band in Fig. 9 shows that with the help of the bare energy differences all three atomic peaks can be discerned again, accompanied by an additional structure at the low-energy side of the Hubbard band, which is reminiscent of the Hubbard side peaks found in the one-and two-band Hubbard model on the Bethe lattice [35] upon increasing U . We leave further investigation of this feature to future work. It might at first seem counter-intuitive that increasing U makes the physics less atomic-like. Indeed, at very high interaction strengths, in the insulating regime, the spectrum must become atomic-like again. Here, however, we identify an intermediate regime where additional structures appear when increasing U , since we get closer to the Mott metal-to-insulator transition.

E. Solution of a five-band AIM
In this section we show that FTPS can not only deal with three-band models, but also work in the case of five orbitals. To do so, we use the bath parameters k and V k from the converged N b = 59 calculation for SrVO 3 and construct an artificial degenerate five-band AIM. Interaction parameters are U = 4.0 eV and J = 0.6 eV. We decrease the on-site energy 0 to get a similar occupation of each impurity orbital as for SrVO 3 ( n m,0,σ ≈ 1 6 ). Note that in doing so we have a model with, in total, 5 3 electrons on the impurity. We only use densitydensity interactions and carry out the time evolution to t = 16 eV −1 . We set the truncated weight to t w = 10 −8 , but restrict the bond dimension of the impurity-impurity links to χ max = 200.
In Fig. 10 we compare the results obtained for this fiveband model to results from CTQMC, where we used the same discretized bath hybridization as input to CTQMC. We again see excellent agreement, even on a logarithmic scale. The spectrum A(ω) (not shown) again exhibits strong structure in the upper Hubbard band. Of course, the computational complexity is larger than in the threeorbital case and it grows during time evolution. Calculating the Greens function took about 190 hours on the processors specified in Appendix 2. We want to stress tough that the resulting spectrum (as well as Fig. 10) was already fully converged at t = 12 eV −1 (70 hours). We note that even with only one CPU hour (t = 6 eV −1 ) the resulting spectrum is almost converged and barely distinguishable from the final result. The benchmark therefore shows that with our FTPS approach a full five-orbital DMFT calculation is well within reach.

V. CONCLUSIONS
We have presented a novel multi-orbital impurity solver which uses a fork-like tensor network whose geometry resembles that of the Hamiltonian. The network structure is simple enough to generate Schmidt decompositions, allowing us to truncate the tensor network safely and to use established methods like DMRG and real time evolution. The solver works on the real frequency axis, and hence allows to formulate the full DMFT self consistency procedure for real frequencies. Therefore, results are not plagued by an ill-conditioned analytic continuation. Our approach exhibits no sign problem, tough it does become more involved for larger numbers of orbitals.
We tested the solver within DMFT on a Hamiltonian typically used for the testbed material SrVO 3 and investigated the influence of full rotational invariance on the results. We found clear spectral structures in particular in the upper Hubbard band that have not been accessible by CTQMC, for which the necessary analytic continuation prohibits the resolution of fine structures in the spectral function at higher energies. For our calculations with U = 4.0 eV, each peak in the spectrum corresponds to an atomic excitation. Even excitations into states with three particles on the impurity are resolved, as tiny spectral peaks at high energies. Furthermore, upon increasing U , an additional structure appears on the inside of the Hubbard bands, similar to the precursors of the sharp Hubbard side peaks found for the one-and two-band Hubbard models on the Bethe lattice [29,35]. We have also shown that our approch is feasible for five-orbital models, by comparing results from the FTPS solver to CTQMC for an artificial five-band model. . Spectrum A(ω) using different linear prediction (LP) parameters for a calculation without spin-flip and pairhopping terms. The calculations with LP were performed with a broadening of ηF T = 0.02 eV. Except for small changes around ω = 0, the effect of the various LP parameters is minor. The blue line directly lies below the red and green line. We also show a DMFT calculation without any LP. Even then the main features are still present.

ACKNOWLEDGMENTS
The authors acknowledge financial support by the Austrian Science Fund (FWF) through SFB ViCoM F41 (P04 and P03), through project P26220, and through the START program Y746, as well as by NAWI-Graz. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. We are grateful for stimulating discussions with F. Verstraete, K. Held, F. Maislinger and G. Kraberger. The computational resources have been provided by the Vienna Scientific Cluster (VSC). All calculations involving tensor networks were performed using the ITensor library [81].

APPENDIX
In this appendix we show that our results are very stable over a wide range of computational parameters. First we focus on the linear prediction method (Sec. 1). Then we show that the results are converged with respect to the usual MPS-approximation (Sec. 2).

Linear Prediction
In order to obtain smooth and sharp spectra, we used linear prediction (LP) to extrapolate the Greens function in time [28,35]. Without going into detail, we state the fact that linear prediction has two parameters, the pseudo inverse cutoff p inv and the order N of the linear prediction. Fig. 11 shows that the results are converged in these parameters. We also show a DMFT run without any linear prediction, which is only possible if we increase the broadening parameter of the Fourier transform to η F T = 0.1 eV, since otherwise we would get oscillations due to the hard cutoff of the time series. Except for a shift towards the right, omitting the linear prediction only changes the height (and width) of the peaks, but not the overall structure. This is a strong indicator of the stability of these features.

Truncation of the Tensor Network
One, if not the most important, parameter in any MPS-like calculation is the sum of discarded singular values in each SVD (truncated weight t w ). We want to emphasize that this parameter is the only approximation in the representation of a state as a tensor-product state, as we do not impose any hard cutoff on the bond dimensions. Fig. 12 shows that the spectrum is well converged with respect to the truncation error during time evolution.
Finally, we want to comment on the required computational effort. In the calculation without full rotational symmetry, the size of the largest tensor to represent the ground state was [82] 35×22×9×2 (t w = 10 −9 ) and at the end of time evolution 127 × 79 × 30 × 2 (t w = 10 −8 ). For a truncated weight of t w = 10 −7 , calculating the Greens function took about 17 minutes on a node with two processors (Intel Xeon E5-2650v2, 2.6 GHz with 8 cores, and G > and G < each calculated on one processor). This time increases to five hours for the lowest truncated weight of t w = 5·10 −9 . Using the full rotationally invariant Hamiltonian, the biggest tensor in the ground-state search was 90 × 40 × 10 × 2 (t w = 10 −8 ) and at the end of time evolution 79 × 46 × 21 × 2 (t w = 2 · 10 −8 ). The Greens function takes about three hours, and we need five hours for one full DMFT-cycle on the same two processors as above.