Fork Tensor-Product States: Efficient Multiorbital Real-Time DMFT Solver

Fork Tensor-Product States: Efficient Multiorbital Real-Time DMFT Solver PHYSICAL REVIEW X 7, 031013 (2017) 1,* 1 1 1 1,2 Daniel Bauernfeind, Manuel Zingl, Robert Triebl, Markus Aichhorn, and Hans Gerd Evertz Institute of Theoretical and Computational Physics, Graz University of Technology, 8010 Graz, Austria Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA (Received 23 December 2016; revised manuscript received 14 April 2017; published 20 July 2017) 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 high spectral resolution at all frequencies. We use a large number (Oð100Þ) 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 three-orbital test-bed material SrVO . 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. Finally, we show that our solver can be applied even to models with five orbitals. This impurity solver offers a new route to the calculation of precise real-frequency spectral functions of correlated materials. DOI: 10.1103/PhysRevX.7.031013 Subject Areas: Condensed Matter Physics I. INTRODUCTION real-frequency spectra, which therefore become broadened, especially at high energies. ED directly provides spectra on Strongly correlated systems are among the most fasci- the real axis, but it is severely limited in the size of the nating objects that solid-state physics has to offer. The Hilbert space, i.e., in the number of bath sites. Quite interactions between constituents of such systems lead to recently, NRG was shown to be a viable three-band solver emergent phenomena that cannot be deduced from the by exploiting non-Abelian quantum number conservation properties of noninteracting particles [1]. [15–17]. NRG works on the real axis and captures the low- One of the most widely used methods to describe energy physics well, but by construction, it has a poor strongly correlated electrons is the dynamical mean-field resolution at higher energies. Another interesting route that theory (DMFT) [2,3]. DMFT treats local electronic corre- has been proposed recently is the one of solvers that tackle lations using a self-consistent mapping of the lattice the problem of exponential growth of the Hilbert space problem onto an effective Anderson impurity model using ideas from quantum chemistry, i.e., the configuration (AIM). Calculating the single-particle spectral function interaction [11,12]. They allow one to go beyond the small of this impurity model in an accurate and efficient way bath sizes of ED, keeping all the advantages, such as the is at the heart of every DMFT calculation. To this end, absence of fermionic sign problems. However, in multi- many numerical methods have been developed or adapted. orbital applications (see the Appendix of Ref. [12]), the These are based, for instance, on continuous-time quantum spectral resolution has so far been restricted by the number Monte Carlo (CTQMC) [4,5], exact diagonalization (ED) of bath sites [Oð20Þ]. [6–8], the numerical renormalization group (NRG) [9,10], Finally, MPS-based techniques like DMRG do not suffer the configuration interaction (CI) [11,12], and also the from a sign problem and can be used on the real-frequency density-matrix renormalization group (DMRG) with axis as well as on the imaginary-frequency axis. The absence matrix-product states (MPS) [13,14]. of the sign problem, in general, comes at the cost of a very Every algorithm has strengths and weaknesses: CTQMC large growth of bond dimension with the number of orbitals. is exact, apart from statistical errors on the imaginary axis, Dynamical properties and spectral functions can be and can deal with multiple orbitals, but it is, in some cases, calculated within DMRG and have been used for impurity plagued by the fermionic sign problem. Additionally, solvers, e.g., with the Lanczos-like continued-fraction an ill-posed analytic continuation is necessary to obtain expansion [18,19]. Other solvers using the more stable correction vector [20] and dynamical DMRG (DDMRG) [21] methods were developed [22–25]. Both algorithms daniel.bauernfeind@tugraz.at produce very accurate spectral functions but have the Published by the American Physical Society under the terms of disadvantage that a separate calculation for each frequency the Creative Commons Attribution 4.0 International license. has to be performed. The Chebyshev expansion [26] with Further distribution of this work must maintain attribution to MPS [27], supplemented by linear prediction [28],was the author(s) and the published article’s title, journal citation, and DOI. used for impurity solvers in the single-band case [29] and 2160-3308=17=7(3)=031013(15) 031013-1 Published by the American Physical Society DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) for two bands [30]. Recently, some of us introduced a real-time calculations. We emphasize that (i) our method method based on real-time evolution [31–34] and achieved is, by construction, free of any fermionic sign problem, a self-consistent DMFT solution for a two-band model (ii) one can fully converge the DMFT self-consistency loop [35]. In such calculations, the physical orbitals for each on the real-frequency axis, and (iii) we can achieve an almost spin direction are usually combined into one large site in exact representation of the bath spectral function. We apply the MPS. Three or more orbitals have not been feasible with this method to multi-orbital DMFT for the test-bed material this approach because of a large increase in computational SrVO and show that one can resolve a multiplet structure in cost with the number of orbitals. Another MPS-based the Hubbard bands, at the same time keeping a good solver, which works on the imaginary axis, was recently description of the low-energy quasiparticle excitations. introduced [36], and it was applied as a solver for three The paper is structured as follows. First, we outline how bands in two-site cluster DMFT. It was supplemented by a impurity solvers with tensor networks work and introduce single real-time evolution to compute the spectral function, our new tensor network approach, which we call fork avoiding the analytic continuation. However, this method is tensor-product states (FTPS) (Sec. II). Next, we explain in restricted by the number of bath sites that can be employed. detail how our solver is used in the context of multi-orbital In the calculation mentioned, only three bath sites per DMFT (Sec. III). In Sec. IV, we apply our approach to orbital were used, limiting the energy resolution for real- SrVO and discuss the multiplet structure that the FTPS frequency spectral functions. solver allows us to resolve. In order to check the accuracy In the present paper, we introduce a novel impurity of the method, we also compare the FTPS results to solver that works directly on the real-frequency axis. To this CTQMC for SrVO . Finally, we show the efficiency of end, we use a tensor network that captures the geometry of the FTPS solver by applying it to a five-orbital model. the interactions in the Anderson model better than a standard MPS. Our approach is, to some extent, inspired II. TENSOR NETWORK IMPURITY SOLVERS by the work of Ref. [37], which used a similar network for a The AIM describes an impurity (with Hamiltonian H ) two-orbital NRG ground-state calculation. We are not loc coupled to a bath of noninteracting fermions hybridized restricted to a small number of bath sites. This is imperative for exploiting the spectral resolution achievable with with it. A typical AIM Hamiltonian is given by H ¼ H þ H ; loc bath H ¼ ϵ n þ H þ H ; loc 0 m0σ DD SF-PH mσ X X X 0 0 H ¼ U n n þðU − 2JÞ n n þðU − 3JÞ n n ; DD m0↑ m0↓ m0σ m 0¯σ m0σ m 0σ 0 0 m >m;σ m >m;σ X X † † † † 0 0 0 H ¼ J ðc c c c þ H:c:Þ − J ðc c c c þ H:c:Þ; SF-PH m0↓ m 0↑ m 0↑ m 0↓ m0↑ m 0↓ m0↑ m0↓ 0 0 m >m m >m H ¼ ϵ n þ V ðc c þ H:c:Þ; ð1Þ bath l mlσ l mlσ m0σ mlσ where c (c ) creates (annihilates) an electron in band GðtÞ¼ −iθðtÞhψ j½c ðtÞ;cð0Þjψ i; ð2Þ mlσ mlσ 0 0 m (m ∈ f1; 2; 3g for a three-orbital model) with spin σ at the lth site of the system (the impurity has index l ¼ 0; the bath degrees of freedom have l ≥ 1), and n are the mlσ of the interacting problem (1), either in real- or imaginary- corresponding particle-number operators. Note that H DD time t. In the present paper, we introduce a new tensor describes density-density (DD) interactions between all network similar to a MPS, which can be used as a real-time orbitals, and H are the spin-flip and pair-hopping SF-PH impurity solver for three orbitals. We first introduce MPS terms. This three-orbital Hamiltonian is not only important before moving on to what we call fork tensor-product states in the context of real-material calculations. It has also in Sec. II B. been studied extensively on the model level—most im- portantly, because it hosts unconventional correlation phe- A. MPS and DMRG nomena. For a selection of recent work, see, for instance, Refs. [15,38–41]. MPS are a powerful tool to efficiently encode quantum- An impurity solver calculates the retarded impurity mechanical states. Consider a state jψi of a system Green’s function GðtÞ, consisting of N sites: 031013-2 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) state jψ i and ground-state energy E . It minimizes the jψi¼ c js ;s ;…;s i: ð3Þ 0 0 s ;…;s 1 2 N 1 N s ;s ;…;s expectation value 1 2 N Each site i has a local Hilbert space of dimension d hψjHjψi E ¼ min ; ð6Þ spanned by the states js i. Through repeated use of i jψi hψjψi singular-value decompositions (SVDs), it is possible to factorize every coefficient c into a product of usually by updating two neighboring MPS tensors before s ;…;s 1 N matrices [14], i.e., into a MPS, moving on to the next bond. This procedure also yields the Schmidt decomposition of the state at the current bond on s s s 1 2 N jψi¼ A · A  A js ;s ;…;s i: ð4Þ the fly. The DMRG approximation keeps only those states 1 2 N 1 2 N s ;s ;…;s with the largest Schmidt coefficient. It is important to note 1 2 N that one can perform a DMRG calculation for any tensor Each A is a rank-three tensor, except the first and last ones network, as long as one can generate a Schmidt decom- s s 1 N (A , A ), which are of rank two. The index s is called a i position [43]. 1 N physical index, and the matrix indices, which are summed For obtaining the Green’s function, we employ an over, are the so-called bond indices. A general state of the evolution in real-time. Equation (2) is split into the greater > < full Hilbert space is not feasible to store, but it can be shown (G ) and lesser Green’s function (G ): that ground states are well described by a MPS with limited > < bond dimension m (dimension of the bond index) [42]. GðtÞ¼ −iΘðtÞðG ðtÞþ G ðtÞÞ; In complete analogy to the states, one can factorize an > −iHt † iE t G ðtÞ¼hψ jce c jψ ie ; 0 0 operator into what is called a matrix-product operator < † iHt −iE t (MPO) [14], 0 G ðtÞ¼hψ jc e cjψ ie ; ð7Þ 0 0 0 0 s ;s s ;s 1 N 1 N 0 0 which are calculated in two separate time evolutions. This H ¼ W  W js ;…;s ihs ;…;s j; ð5Þ 1 N 1 N 1 N s ;…;s 1 N is done by first applying c (or c) and then time evolving 0 0 s ;…;s 1 N this state and calculating the overlap with the state at time s ;s t ¼ 0. The time evolution is the most computationally where each W is a rank-four tensor. Tensor networks, in expensive part since time-evolved states are not ground general, have a very useful graphical representation, which states anymore, and the needed bond dimensions usually is shown for a MPS and a MPO in Fig. 1. Note that when grow very fast with time. we use the term MPS, we always mean a one-dimensional chain of tensors, as shown in Fig. 1. B. Fork tensor-product states To calculate Green’s functions within the MPS formal- ism, one usually first applies the DMRG [13,14], which So far, the usual way of dealing with Hamiltonians like acts on the space of MPS and finds a variational ground 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, respectively. The local state space of each bath site then (a) consists of M spinless-fermion degrees of freedom, with dimension 2 , 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 encoun- ter this very fast growth upon calculating the ground state (b) 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 FIG. 1. (a) Graphical representation of a MPS. Every circle usually transformed into a Wilson chain with nearest- corresponds to a tensor A and each line to an index of this tensor. neighbor hoppings only, i.e., of the form t ðc c þ i iþ1 In this picture, the physical indices are the vertical lines, while the H:c:Þ [10]. This was thought to be necessary since long- horizontal lines show the bond indices. Connected lines mean range interactions look problematic for MPS-based algo- that the corresponding index is summed over. Fixing all the rithms. Quite recently, though, it was discovered that MPS physical indices s for each site results in a tensor of rank zero can deal with the original form of H in Eq. (1) better bath with the value of the coefficient c . (b) Graphical repre- s ;…;s 1 N [44]. Because all hopping terms in H originate from the bath sentation of a MPO. The difference from a MPS is that a MPO has impurity, this is called the star geometry. The reason for the incoming indices s and outgoing indices s corresponding to the bra- and ket vectors of the operator. better performance is that in the star geometry, one has 031013-3 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) many nearly fully occupied (empty) bath sites with very tensors representing the impurities. When all bond indices low (high) on-site energies ϵ . have the same dimension χ, it is necessary to do a SVD for 2 4 Since basis states with many unoccupied low-energy sites a χ d × χ matrix with computational complexity Oðχ dÞ. have a very low Schmidt coefficient, these states are However, as we show below, this operation does not pose a discarded from the MPS. The same holds for occupied substantial problem for calculations using FTPS. Since the high-energy sites. However, when dealing with multi-orbital impurity tensors pose the biggest challenge, our tensor models, the star geometry is not enough to be able to network would likely also allow us to deal with the chain calculate Green’s functions using MPS. The growth of the geometry without a drastic increase in computational cost. In the present paper, we only use FTPS with baths in star bond dimensions still makes those calculations unfeasible. geometry. The proposed FTPS are similar to the tensor The key idea of the present work is to construct a tensor network used by Holzner et al. [37] to perform NRG network which is beyond a standard MPS but similar calculations for ground-state properties of an AIM with two enough to be able to use established methods like DMRG orbitals. and time evolution. From Hamiltonian (1), one can The three-legged tensors in our network (Fig. 2) can also immediately notice that there are no terms coupling bath be interpreted as two coupled junctions with three legs in sites of different orbitals. Hence, it might not be advanta- the language of Ref. [45], where it has been shown that geous to combine those (not directly interacting) degrees of DMRG is possible on such junctions. Furthermore, our freedom into one large physical index in the MPS. approach has similarities to the so-called tree tensor net- Therefore, our proposed tensor network separates the works (TTN) [43,46–48]. 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 1. Time evolution network is not a MPS anymore since there are some tensors Time evolution with the Hamiltonian Eq. (1) is not (labeled A↓ and B↑ in the example of Fig. 2) that have straightforward since it features long-range hoppings. three bond indices and one physical index, i.e., which are of Possible methods include Krylov approaches [49], the rank four. Cutting any bond splits the network into two time-dependent variational principle [50,51], and the series separate parts. Therefore, one can calculate the Schmidt iHt expansion of e proposed by Zaletel et al. [52]. In this decomposition in a way very similar to a MPS, which work, however, we use a much simpler approach. means that DMRG is also possible. The main bottleneck of First, we split the Hamiltonian into the following terms: calculations with FTPS is to perform SVDs of the rank-four SF-PH (i) the spin-flip and pair-hopping terms H for each P m;m SF-PH orbital combination, with 0 H ¼ H [see SF-PH m >m m;m Eq. (1)], (ii) the density-density interaction terms H , DD and (iii) all other terms, H ¼ H þ ϵ n . With free bath 0 m0σ mσ these terms, we write the time-evolution operator for a small time step Δt using a second-order Suzuki-Trotter decomposition [53], SF-PH −i½ðΔtÞ=2H −iΔtH 0 −i½ðΔtÞ=2H m;m DD e ≈ e e · m >m SF-PH −i½ðΔtÞ=2H −iΔtH −i½ðΔtÞ=2H 0 free DD m;m ×e e e : ð8Þ m >m Note that in this decomposition, the order of the spin-flip and pair-hopping terms is important. The order of operators in FIG. 2. Graphical representation of a FTPS for multi-orbital AIM. Separating bath degrees of freedom leads to a forklike the second product must be opposite to the one in the first. structure. In this picture, a two-orbital model is shown, with four We see that Eq. (8) involves three different operators, SF-PH bath sites in each orbital. Orbitals are labeled A and B, and the H , H , and H , each of which will be treated DD free m;m arrows denote the spin. Each spin-orbital combination has its own differently. bath sticking out to the right. As in Fig. 1, the vertical lines are the Time evolution of the density-density interactions is physical degrees of freedom [all of dimension two, for empty performed with a MPO-like representation of the time- (respectively, occupied) bath sites]. All other lines are bond −i½ðΔtÞ=2H DD evolution operator e . For a three-orbital model, indices, and like in the MPS, they are summed over. As 3 3 −i½ðΔtÞ=2H DD first the full matrix (4 × 4 )of e is created, mentioned in the text, the bath is represented in star geometry which is then decomposed into MPO form by repeated because of the smaller bond dimensions needed. The bath sites SVDs. Since H only consists of density-density inter- are ordered according to their on-site energies. Two example DD −iΔtH DD hoppings V and V are drawn. actions, no fermionic sign appears in e . 1 2 031013-4 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) Time evolution of the spin-flip and pair-hopping terms is SVD(toseparatethe tensorsagain),the swapgateisappliedso more involved than the density-density interactions since that the impurity and the first bath sites are swapped. By the operators change the particle numbers on the impurity repeating this process, one moves the impurity along its sites. Therefore, it can be difficult to deal with the fermionic horizontal arm in Fig. 2. Because a second-order decompo- sign of the time-evolution operator when the impurities are sitionisused,nowalltime-evolutiongatesexcepttheoneatsite not next to each other in the fermionic order. It turns out N have to be applied again. But now, the impurity and bath that the spin-flip and the pair-hopping terms have the sitesneedtobeswapped before time evolution. ˆ ˆ ˆ Note that not only can the algorithm presented above be property A ¼ A individually, with A being either the spin- flip or the pair-hopping operator, respectively. Furthermore, used to perform real-time evolutions, but it is also appli- they commute with each other, allowing us to separate them cable to evolution in imaginary-time simply by replacing without Trotter error. The time-evolution operator of JA is idt by dτ. then given by ˆ III. MULTI-ORBITAL DMFT WITH FTPS −iΔtJA ˆ ˆ e ¼ 1 þ A ðcos ðΔtJÞ − 1Þ − iA sin ðΔtJÞ: ð9Þ In this section, we present details of our impurity solver. For this operator, a MPO can be found for which the We refer to Refs. [2,56] for DMFT in general, and to fermionic sign can easily be determined [54]. Refs. [57,58] for DMFT in the context of realistic ab initio To time evolve the bath terms, we use an iterative second- calculations for correlated materials. order Suzuki-Trotter breakup for each term in H . free In the latter approach, called density-functional theory Neglecting orbital (m) and spin (σ) indices, the first (DFT)+DMFT, the correlated subspace is described by a −iΔt H l¼1 Hubbard-like Hamiltonian. Within DMFT, this model is step in this breakup is the following: e ≈ −iΔt H mapped onto the AIM Hamiltonian (1). This mapping −i½ðΔtÞ=2H −i½ðΔtÞ=2H 1 l¼2 1 e e e . Next, we split off H defines the bath hybridization function ΔðωÞ describing the and iterate this procedure until we end up with influence of the surrounding electrons. N −1 b Since FTPS provide the Green’s function of the AIM on Y Y Δt −iΔtH −i H free mlσ e ≈ e the real-frequency axis, the self-consistency loop is also mσ l¼1 performed directly for real frequencies. For calculating the bath hybridization, we use retarded Green’s functions with Δt −iΔtH mN σ −i H mlσ b 2 × e e ; ð10Þ a finite broadening η in order to avoid numerical SC l¼N −1 difficulties with the poles of the Green’s function. Throughout this work, we use η ¼ 0.005 eV [59]. SC with N thenumberofbathsitesand H ¼ b mlσ The impurity solver calculates the self-energy ΣðωÞ of ϵ n þ V ðc c þ H:c:Þ. In the above equation, we l mlσ l mlσ m0σ the AIM, given the bath hybridization function ΔðωÞ and neglected the term ϵ n that we add to H . 0 m0σ m1σ the interaction Hamiltonian on the impurity. To this end, Equation (10) is a product of two-site gates (an operator acting our solver performs the following steps, which are nontrivially only on two sites) with one of the two sites always explained in more detail in the text below: being the impurity. This means that those two sites are not (1) Obtain bath parameters ϵ and V by a deterministic l l nearest neighbors in the tensor network. To overcome this approach based on integration of the bath hybridi- problem, we use so-called swap gates [14,55]. The two-site zation function ΔðωÞ. operator (2) Calculate the ground state jψ i and ground-state energy E of the interacting problem. n n i j 0 0 0 S ¼ δ δ ð−1Þ ð11Þ ij s ;s s ;s i j j i (3) Apply impurity creation or annihilation operators, and time evolve these states to determine the swaps the state of site i (s with occupation n ) with the i i interacting Green’s function [Eq. (2)]. n n i j state of sitej (s with occupationn ).The factor ð−1Þ gives j j (4) Fourier transform Eq. (2) to obtain GðωÞ and the correct fermionic sign and is negative if an odd number of calculate the local self-energy, particlesonsiteigetsswappedwithanoddnumberofparticles −1 −1 on site j. To be more precise, the matrix representation of the ΣðωÞ¼ G ðωÞ − GðωÞ : ð13Þ swap gates used in this work is To perform step 1 we use S ¼j00ih00jþj10ih01 þj01ih10j − j11ih11j: ð12Þ ij h i V ¼ − ImΔðωÞ dω; Itturnsoutthateveryswapgatecanbecombinedwithanactual l time-evolution gate without additional computational time. h i 1 1 For example, the first step in this time evolution would be to ϵ ¼ ω − ImΔðωÞ dω; ð14Þ −i½ðΔtÞ=2H V π m1σ l l apply e . Immediately afterwards, even before the 031013-5 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) similar to Refs. [10] (NRG) and [44]. Each interval I parameter sets. The next step of our testing was to include corresponds to a bath site. This discretization can be density-density interactions, one term at a time. For interpreted as representing each interval I as a delta peak example, we only included ðU − JÞn n and compared l 10↑ 30↑ at position ϵ and weight V . Sum rules for such discre- energy and Green’s functions to a standard one-orbital MPS tization parameters can be found analytically [60]. In this solver. Finally, we also compared our method to the MPS two-band solver used in Ref. [35]. Indeed, all tests work, we choose the length of each interval such that the performed produced the correct energies and Green’s area of the bath spectral function −ð1=πÞImΔðωÞ is approximately constant for each interval [61]. For the case functions. at hand, this discretization was found to be numerically more stable than using intervals of constant length. Unless IV. RESULTS stated otherwise, we use N ¼ 109 bath sites per orbital We performed DMFT calculations based on a band and spin. We note that this scheme is not restricted to structure obtained from DFT for the prototypical com- diagonal hybridizations. In the general case of off-diagonal pound SrVO , using the approximation of the Kanamori hybridizations, the hybridization function is a matrix Δ. Hamiltonian [Eq. (1)]. It has a cubic crystal structure with a Therefore, instead of taking the imaginary part, we can use nominal filling of one electron in the V-3d shell [64]. the bath spectral function ½i=ð2πÞðΔ − Δ Þ. Similarly to Because of the crystal symmetry, the five orbitals of the Eq. (14), we represent each interval by one delta peak for V-3d shell split into two e and three t orbitals. The latter each orbital. For instance, fixing ϵ to the center of g 2g form the correlated subspace. We performed the DFT the interval, the hopping parameters V can be found calculation with Wien2k [65] and used 34220 k points systematically from the Cholesky factorization of in the irreducible Brillouin zone in order to reach an ½i=ð2πÞðΔ − Δ Þdω. Most importantly, this scheme energy resolution comparable with the η ¼ 0.005 eV does not involve any fitting procedure on the Matsubara SC broadening. axis. A very similar approach was developed independently The TRIQS/DFTTools package (v1.4) [66–68], which is in Ref. [62]. based on the TRIQS library (v1.4) [69], was used to In step 2, we use a DMRG approach with the following generate the projective Wannier functions and to perform parameters, unless specified otherwise. The truncated the DMFT self-consistency cycle. weight t (sum of all discarded singular values of each −8 Figure 3 shows the main results of this paper: the DMFT SVD) is kept smaller than 10 . When spin-flip and pair- spectral function AðωÞ for SrVO , (i) in the approximation hopping terms are neglected, we use an even smaller cutoff −9 of density-density interactions only and (ii) with full of 10 . Note that, except in the five-band calculation, we rotational invariance, including spin-flip and pair-hopping do not restrict the bond dimensions by some hard cutoff terms. Overall, both cases show the well-known features of (see Appendix A2). the SrVO spectral function [70,71]. We see a hole During time evolution (step 3), we use a truncated weight 3 −8 −8 excitation at around −2 eV, as well as the quasiparticle of t ¼ 2 × 10 ,or 10 with density-density interactions −1 peak at zero energy whose shape and position does not only. We time evolve to t ¼ 16 eV , with a time step of −1 Δt ¼ 0.01 eV . Green’s functions are measured every fifth time step. The time-evolution operator of H is loc applied using the zip-up algorithm [55]. Afterwards, the Green’s functions are extrapolated in time using the linear- −1 prediction method [28,35] up to t ¼ 250 eV . 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 iωt−η jtj FT broadening in the kernel e of η ¼ 0.02 eV to FT avoid cutoff effects remaining after the linear prediction. The influence of the linear prediction on our results is discussed in Appendix A1. We stress that although a calculation with full rotational symmetry is more demand- ing, the computational effort is still very feasible. With the parameters mentioned above, one full DMFT cycle takes FIG. 3. Spectral functions AðωÞ for density-density interactions about five hours on 16 cores. (DD) only (blue line), and with spin-flip and pair-hopping To verify that our implementation of DMRG and time terms included (red line). In both calculations, we used U ¼ evolution produces correct results when used with our 4.0 eV and J ¼ 0.6 eV. Both spectra show a three-peak structure tensor network, we first compared Green’s functions and in the upper Hubbard band and additional features at high ground-state energies for U ¼ J ¼ 0 for several bath energies (around 8 eV). 031013-6 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) depend on the inclusion of full rotational invariance. In the (a) upper Hubbard band, a distinctive three-peak structure can be seen. This structure has not been resolved in other exact methods like CTQMC (because of problems with analytic continuation, see below) or NRG (because of the loga- rithmic discretization problem). In our real-time approach, high energies correspond to short times, where the calcu- lations are particularly precise [72]. Most methods allow us 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 noncrossing approximation (NCA) show atomiclike features, but they have very limited accuracy for the description of the low-energy quasiparticle excitations in the metallic phase [74]. Thus, our FTPS solver combines the best of both worlds, with atomic multiplets at high (b) 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 FIG. 4. (a) We take the bath spectral function ΔðωÞ from the between impurities. In all subsequent plots, we show results DMFT self-consistent solution for N ¼ 109 and represent it from calculations with DD interactions only. b using various numbers of bath sites. It is obvious that N ¼ 9 is too small to represent the bath well. (b) Converged DMFT spectral function using the AIM with different numbers of bath A. Effect of bath size sites. Only the smallest bath shows a noticeable difference. This is mostly due to the fact that, in this case, a higher broadening of In order to achieve a reliable high-resolution spectrum η ¼ 0.1 eV had to be used in the Fourier transform and time on the real-frequency axis, it is imperative to have a good FT −1 evolution was only possible for t ¼ 14 eV . The additional representation of the hybridization function ΔðωÞ in small structure at ω ¼ 0 for N ¼ 59 bath sites is most likely a terms of the bath parameters, for which a sufficient linear-prediction artifact. number of bath sites is needed. Figure 4 shows how well a hybridization function can be represented with our B. Comparison to CTQMC approach [Eq. (14)] using a certain number of bath sites. We see that for N ¼ 9 bath sites (we always denote sites In Fig. 5, we compare the converged spectral function of per orbital), ΔðωÞ can be reconstructed only very our approach (FTPS) with a spectrum obtained from roughly, which in turn gives an incorrect spectral function CTQMC and analytic continuation. In both calculations, (Fig. 4, bottom). To some extent, the difference in the we used the same interaction Hamiltonian with density- spectrum is due to the shorter time evolution and there- density interactions only. The CTQMC calculation was fore a higher broadening, η , which we were forced to performed with the TRIQS CTHYB solver (v1.4) [75,76] FT use. For such a small bath, the finite-size effects from with 3.2 × 10 measurements and at inverse temperature −1 reflections at the bath ends appear much earlier in the β ¼ 200 eV . For the analytic continuation, we applied time evolution. the Ω MaxEnt method [77]. Increasing the number of bath sites to N ¼ 29,we The three-peak structure in the upper Hubbard band is observe that the reconstructed bath spectral function not present in the CTQMC spectrum. We show in an already shows the relevant features of ΔðωÞ. The spectrum example below that even for a Green’s function that does is well converged for the largest bath sizes N ¼ 59 contain these peaks, the analytic continuation does not and N ¼ 109. resolve this structure. 031013-7 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) FIG. 5. DMFT spectral functions AðωÞ from CTQMC þ −1 MaxEnt (blue line) at β ¼ 200 eV , and from FTPS (red line). The FTPS result shows a distinctive three-peak structure in the upper Hubbard band. For another comparison, we consider the imaginary-time FIG. 6. Comparison of imaginary-time Green’s functions GðτÞ Green’s functions GðτÞ in Fig. 6. Apart from the effect of from CTQMC (G , blue line) and FTPS using Eq. (15) CTQMC statistical errors, CTQMC provides an exact self-consistent (G , red squares). The agreement is also equally good at β ¼ AðωÞ solution of DMFT on the imaginary-frequency axis. As −1 −1 100 eV and β ¼ 400 eV (not shown). The difference be- mentioned above, when we use the FTPS solver, we tween the two Green’s functions is shown in the bottom panel. formulate the DMFT self-consistency equations on the Note that on both ends, G is smaller than G . The AðωÞ CTQMC real axis. To obtain an approximate finite-temperature normalization of the spectral function demands that imaginary-time Green’s function from FTPS that we can Gðτ ¼ 0Þþ Gðτ ¼ βÞ¼ −1. The CTQMC data deviate in the −2 compare to the CTQMC result, we need to take the finite order of 10 from this constraint due to statistical noise, while temperature of the CTQMC calculation into account. FTPS gives (by construction) the correct result to a precision of −8 10 . This explains the bigger differences of the Green’s Therefore, we use the FTPS spectrum AðωÞ and assume functions around τ ¼ 0 and τ ¼ β. For better visibility of the that we would obtain the same spectrum for a finite (but −3 τ > 0 data, the value of 9 × 10 at τ ¼ 0 is not shown. high enough) inverse temperature β; we use −ωτ energies. Further benchmarks concerning the linear pre- dω e GðτÞ¼ AðωÞ : ð15Þ −βω diction can be found in Appendix A1. 2π e þ 1 Finally, we show that the ill-posedness of the analytic continuation is the most likely explanation for the missing 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 Green’s function GðτÞ and Fourier transform it to get Gðiω Þ: iω τ Gðiω Þ¼ e GðτÞdτ: Looking at the last few DMFT cycles, we estimate it to be around Aðω ¼ 0Þ¼ −ð1=πÞlim ℑGðiω Þ ≈ iω →0 n −1 0.272 eV , 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 FIG. 7. Spectral functions from analytically continued −1 Aðω ¼ 0Þ¼ 0.28 eV , with fluctuations of about imaginary-time Green’s functions GðτÞ calculated by CTQMC −1 0.01 eV . This agreement is very good considering that (blue line) and by FTPS (red line). Clearly, the analytic continu- linear prediction has its strongest influence at small ation cannot resolve the peak structure in the upper Hubbard band. 031013-8 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) TABLE I. Relevant states of the atomic problem of Hamiltonian (1) without spin-flip and pair-hopping terms. Type States Energy difference to ground state Degeneracy N ¼ 1, ground state j↑; 0; 0ij↓; 0; 0ij0; ↑; 0i 06 N ¼ 0 j0; 0; 0i −ϵ 1 N ¼ 2, same spin j↑; ↑; 0ih↑; 0; ↑ij0; ↑; ↑i U − 3J þ ϵ 6 N ¼ 2, different spin j↑; ↓; 0ij↑; 0; ↓ij↓; ↑; 0i  U − 2J þ ϵ 6 N ¼ 2, double occupation j↑↓; 0; 0ij0; ↑↓; 0ij0; 0; ↑↓i U þ ϵ 3 N ¼ 3, all spins equal j↑; ↑; ↑ih↓; ↓; ↓i 3U − 9J þ 2ϵ 2 N ¼ 3, one spin different j↑; ↑; ↓ij↑; ↓; ↑ij↓; ↑; ↑i  3U − 7J þ 2ϵ 6 N ¼ 3, double occupation j↑↓; ↑; 0ij↑↓; ↓; 0ij↑↓; 0; ↑i 3U − 5J þ 2ϵ 12 peak structure in the upper Hubbard band of the spectral Table II shows how bare atomic parameters change when function obtained from the CTQMC data. To do so, we take adding a bath, and we see that our qualitative arguments the FTPS spectrum AðωÞ, calculate GðτÞ as described give a correct idea of how parameters are rescaled. Further evidence that the observed three-peak structure is above, and perform the same analytic continuation that we did for the GðτÞ from CTQMC. We added noise of the indeed a result of atomic physics can be seen in Fig. 8. This order of the CTQMC error to the FTPS data. The resulting figure shows a close-up of the upper Hubbard band for spectrum is shown in Fig. 7, and indeed, the peak structure three different values of J. The corresponding effective in the upper Hubbard band vanishes. parameters J are shown in Table II. We observe that J is also rescaled slightly, but the rescaling gets smaller for higher J. Furthermore, increasing J also increases the total C. Discussion of peak structure—Effective width of the Hubbard band, which scales mostly linearly atomic physics with J. At the same time, measuring the quasiparticle In order to understand the peak structure observed in the spectral weight as a function of J at constant U shows that it spectral functions, we take a look at the underlying atomic increases with increasing J, also implying an increasing problem, where, for simplicity, we start with density- critical U for the metal-to-insulator transition [40]. density interactions only. We show that the same arguments Upon a careful inspection of the spectral function in hold for full rotationally invariant interactions. Fig. 5, we observe small peaks at energies around 8 eV. A Table I shows the relevant atomic states and their close-up of this energy region for different values of J is corresponding energies. The atomic model has a hole shown in Fig. 8. The energy difference between the peaks is excitation at energy −ϵ and three single electron excita- close to 2J and can, again, be well explained by atomic tions with energies U þ ϵ , U − 2J þ ϵ and U − 3J þ ϵ 0 0 0 physics, namely, excitations into states with three electrons relative to the ground state. If we measure the energy on the impurity (Table I) [79]. These excitations originate differences between the three peaks of the upper Hubbard from small admixtures of N ¼ 2 states to the ground state. band in our results, we find values of 1.27 eV and 0.69 eV, With atomic physics in mind, let us take a look again at which are close to the atomic energy differences of 1.2 eV the spectrum using full rotational symmetry (Fig. 3). The and 0.6 eV (J ¼ 0.6 eV). We also find the hole excitation at spin-flip and pair-hopping terms only contribute if there −2.0 eV. This indicates that we can describe the positions are two or more particles present. Thus, the quasiparticle of the observed peaks approximately by atomic physics peak and the hole excitation do not change. The atomic ¯ ¯ with effective parameters ϵ , U, and J and widened peaks. N ¼ 2 sector does change, however. Diagonalizing the Furthermore, the heights of the peaks roughly correspond to the degeneracy of the states in the atomic model (see Table I). TABLE II. Atomic parameters and their effective values ob- We can determine U ¼ 5.97 eV (where U ¼ 4.00 eV) tained from the spectral functions shown in Figs. 5 and 8. For J, from the energy difference of the peak highest in energy to the values themselves were obtained from the energy difference the hole excitation. This increase of U compared to U is of the highest peak to the lowest peak, whereas the uncertainty is plausible, considering the following. When coupling the estimated from ω þ 2J and ω − J. M M impurity to the bath, particles have the possibility to avoid Parameter Atomic value (eV) Effective value (eV) each other by jumping into unoccupied sites of the bath. ϵ −0.86 −2.00 This results in a decrease of hn n i. To model this situation 0 ↑ ↓ U 4.00 5.97 using atomic physics, one needs to increase the interaction J 0.50 0.59(6) strength. Finally, it is well known that J is much less J 0.60 0.66(3) affected by the surrounding electrons than U since the latter J 0.70 0.72(2) is screened significantly stronger [78]. 031013-9 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) regime. In Fig. 9, we show results with U ¼ 5.5 eV at (a) constant J ¼ 0.6 eV. We see a shift of the upper Hubbard band to higher energies, but little shift of the hole excitation. Also, the central peak is shifted and gets slimmer since more weight is transferred into the Hubbard bands. Most importantly, as we approach the strongly correlated 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 Hubbard bands smears out and even vanishes. The close-up 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 addi- (b) tional 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 models on the Bethe lattice [35] upon increasing U. We leave further investigation of this feature to future work. (a) FIG. 8. (a) Close-up of the three-peak structure for various values of J. Additionally, we show vertical lines for the J ¼ 0.6 eV spectrum at energies ω (position of the middle peak) and at ω þ 2J and ω − J. We see that the width of the upper M M Hubbard band is close to 3J. (b) Close-up of the small spectral peaks at high energies. These correspond to excitations into the N ¼ 3 sector of the atomic model (see Table I). The height of each peak can be estimated by the degeneracy of the atomic states. Effective parameters J are 0.53 eV (J ¼ 0.5 eV), 0.59 eV (J ¼ 0.6 eV), and 0.68 eV (J ¼ 0.7 eV). They are obtained from the difference between the two peaks that are highest in energy. (b) Hamiltonian, we find eigenstates with three different energies and differences of 3J ¼ 1.8 eV and 2J ¼ 1.2 eV, respectively. Measuring the energy differences in ¯ ¯ Fig. 3, we find 3J ¼ 1.75 eV and 2J ¼ 1.32 eV. Estimating U ¼ 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 effec- tive 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. FIG. 9. (a) Increasing U results in a slimmer central peak and a shift of the upper Hubbard band. Also, the three-peak structure D. Beyond atomic physics gets smeared out. (b) Close-up of the upper Hubbard band. As in The previous section showed that at U ¼ 4.0 eV, the Fig. 8, additional vertical lines are plotted at ω (position of the spectral features in the Hubbard bands can be well middle peak) and at ω þ 2J and ω − J as a rough guide to M M described by atomic physics with effective parameters where the atomic peaks would be located. With the help of these and widened peaks. It is not clear whether this picture is lines, one can discern a three-peak structure again, but it is valid for higher interaction strengths U in the metallic extended by a feature at the low energy side of the Hubbard band. 031013-10 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) It might at first seem counterintuitive that increasing U Green’s function took about 190 hours on the processors makes the physics less atomiclike. Indeed, at very high specified in Appendix A2. We want to stress, though, that interaction strengths, in the insulating regime, the spectrum the resulting spectrum (as well as Fig. 10) was already fully −1 must become atomiclike again. Here, however, we identify converged at t ¼ 12 eV (70 hours). We note that even −1 an intermediate regime where additional structures appear with only one CPU hour (t ¼ 6 eV ), the resulting when increasing U since we get closer to the Mott metal-to- spectrum is almost converged and barely distinguishable insulator transition. from the final result. The benchmark therefore shows that with our FTPS approach, a full five-orbital DMFT calcu- lation is well within reach. E. Solution of a five-band AIM In this section, we show that FTPS not only applies to three-band models, but it also works in the case of five V. CONCLUSIONS orbitals. To do so, we use the bath parameters ϵ and V k k We have presented a novel multi-orbital impurity solver from the converged N ¼ 59 calculation for SrVO and b 3 which uses a forklike tensor network whose geometry construct an artificial degenerate five-band AIM. resembles that of the Hamiltonian. The network structure is Interaction parameters are U ¼ 4.0 eV and J ¼ 0.6 eV. simple enough to generate Schmidt decompositions, We decrease the on-site energy ϵ to get a similar allowing us to truncate the tensor network safely and to occupation of each impurity orbital as for SrVO use established methods like DMRG and real-time evolu- ðhn i ≈ Þ. Note that, in doing so, we have a model m;0;σ tion. The solver works on the real frequency axis and hence with, in total, electrons on the impurity. We only use allows us to formulate the full DMFT self-consistency density-density interactions and carry out the time evolu- procedure for real frequencies. Therefore, results are not −1 tion to t ¼ 16 eV . We set the truncated weight to plagued by an ill-conditioned analytic continuation. Our −8 t ¼ 10 , but we restrict the bond dimension of the approach exhibits no sign problem, though it does become impurity-impurity links to χ ¼ 200. max more involved for larger numbers of orbitals. In Fig. 10, we compare the results obtained for this five- We tested the solver within DMFT on a Hamiltonian band model to results from CTQMC, where we used the typically used for the test-bed material SrVO and inves- same discretized bath hybridization as input to CTQMC. tigated the influence of full rotational invariance on the We again see excellent agreement, even on a logarithmic results. We found clear spectral structures, in particular, in scale. The spectrum AðωÞ (not shown) again exhibits strong the upper Hubbard band that have not been accessible by structure in the upper Hubbard band. Of course, the CTQMC, for which the necessary analytic continuation computational complexity is larger than in the three-orbital prohibits the resolution of fine structures in the spectral case, and it grows during time evolution. Calculating the 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 approach is feasible for five-orbital models, by comparing results from the FTPS solver to CTQMC for an artificial five-band model. ACKNOWLEDGMENTS The authors acknowledge financial support by the Austrian Science Fund (FWF) through SFB ViCoM F41 (P04 and P03), through Project No. P26220, and through the START Program No. 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. FIG. 10. Comparison of imaginary-time Green’s functions GðτÞ Held, F. Maislinger, and G. Kraberger. The computational from CTQMC (G , blue line) and FTPS using Eq. (15) CTQMC resources have been provided by the Vienna Scientific (G , red squares). As in Fig. 6, they compare very well. AðωÞ 031013-11 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) Cluster (VSC). All calculations involving tensor networks were performed using the ITensor library [81]. APPENDIX A: INFLUENCE OF COMPUTATIONAL PARAMETERS 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. A1). Then, we show that the results are converged with respect to the usual MPS approximation (Sec. A2). 1. Linear prediction In order to obtain smooth and sharp spectra, we used FIG. 12. Different truncation values t in the time evolution do linear prediction to extrapolate the Green’s function in time not influence the shape of the spectrum AðωÞ. [28,35]. Without going into detail, we state the fact that linear prediction has two parameters, the pseudo-inverse cutoff p and the order N of the linear prediction. Figure 11 parameter is the only approximation in the representation of inv shows that the results are converged in these parameters. a state as a tensor-product state, as we do not impose any We also show a DMFT run without any linear prediction, hard cutoff on the bond dimensions. Figure 12 shows that which is only possible if we increase the broadening the spectrum is well converged with respect to the trunca- parameter of the Fourier transform to η ¼ 0.1 eV since FT tion error during time evolution. otherwise we would get oscillations due to the hard cutoff Finally, we comment on the required computational of the time series. Except for a shift towards the right, effort. In the calculation without full rotational symmetry, omitting the linear prediction only changes the height (and the size of the largest tensor to represent the ground state −9 width) of the peaks but not the overall structure. This is a was 35 × 22 × 9 × 2 (t ¼ 10 ) [82] and at the end of −8 strong indicator of the stability of these features. time evolution 127 × 79 × 30 × 2 (t ¼ 10 ). For a trun- −7 cated weight of t ¼ 10 , calculating the Green’s function 2. Truncation of the tensor network took about 17 minutes on a node with two processors (Intel Xeon E5-2650v2, 2.6 GHz with eight cores, and G and One of the most important parameters in any MPS-like G each calculated on one processor). This time increases calculation is the sum of discarded singular values in each to five hours for the lowest truncated weight of SVD (truncated weight t ). We want to emphasize that this −9 t ¼ 5 × 10 . Using the full rotationally invariant Hamiltonian, the biggest tensor in the ground-state search −8 was 90 × 40 × 10 × 2 (t ¼ 10 ), and at the end of time −8 evolution, 79 × 46 × 21 × 2 (t ¼ 2 × 10 ). The Green’s function takes about three hours, and we need five hours for one full DMFT cycle on the same two processors as above. [1] P. W. Anderson, More Is Different, Science 177, 393 (1972). [2] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical Mean-Field Theory of Strongly Correlated Fermion Systems and the Limit of Infinite Dimensions, Rev. Mod. Phys. 68, 13 (1996). [3] W. Metzner and D. Vollhardt, Correlated Lattice Fermions FIG. 11. Spectrum AðωÞ using different linear prediction in d ¼ ∞ Dimensions, Phys. Rev. Lett. 62, 324 (1989). parameters for a calculation without spin-flip and pair-hopping [4] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. terms. The calculations with linear prediction were performed Troyer, and P. Werner, Continuous-Time Monte Carlo with a broadening of η ¼ 0.02 eV. Except for small changes Methods for Quantum Impurity Models, Rev. Mod. Phys. FT around ω ¼ 0, the effect of the various linear prediction param- 83, 349 (2011). eters is minor. The blue line lies directly below the red and green [5] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. lines. We also show a DMFT calculation without any linear Millis, Continuous-Time Solver for Quantum Impurity prediction. The main features are still present. Models, Phys. Rev. Lett. 97, 076405 (2006). 031013-12 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) [6] M. Caffarel and W. Krauth, Exact Diagonalization [24] M. Karski, C. Raas, and G. S. Uhrig, Single-Particle Approach to Correlated Fermions in Infinite Dimensions: Dynamics in the Vicinity of the Mott-Hubbard Metal-to- Mott Transition and Superconductivity, Phys. Rev. Lett. 72, Insulator Transition, Phys. Rev. B 77, 075116 (2008). 1545 (1994). [25] M. Karski, C. Raas, and G. S. Uhrig, Electron Spectra Close [7] M. Capone, L. de’ Medici, and A. Georges, Solving the to a Metal-to-Insulator Transition, Phys. Rev. B 72, 113110 Dynamical Mean-Field Theory at Very Low Temperatures (2005). Using the Lanczos Exact Diagonalization, Phys. Rev. B 76, [26] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The 245116 (2007). Kernel Polynomial Method, Rev. Mod. Phys. 78, 275 [8] J. Kolorenč, A. B. Shick, and A. I. Lichtenstein, Electronic (2006). Structure and Core-Level Spectra of Light Actinide Diox- [27] A. Holzner, A. Weichselbaum, I. P. McCulloch, U. Schollwöck, and J. von Delft, Chebyshev Matrix Product ides in the Dynamical Mean-Field Theory, Phys. Rev. B 92, 085125 (2015). State Approach for Spectral Functions, Phys. Rev. B 83, [9] K. G. Wilson, The Renormalization Group: Critical Phe- 195115 (2011). nomena and the Kondo Problem, Rev. Mod. Phys. 47, 773 [28] S. R. White and I. Affleck, Spectral Function for the S ¼ 1 (1975). Heisenberg Antiferromagnetic Chain, Phys. Rev. B 77, [10] R. Bulla, T. A. Costi, and T. Pruschke, Numerical Renorm- 134437 (2008). [29] M. Ganahl, P. Thunström, F. Verstraete, K. Held, and H. G. alization Group Method for Quantum Impurity Systems, Rev. Mod. Phys. 80, 395 (2008). Evertz, Chebyshev Expansion for Impurity Models Using [11] Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort, Matrix Product States, Phys. Rev. B 90, 045144 (2014). Efficient Real-Frequency Solver for Dynamical Mean-Field [30] F. A. Wolf, I. P. McCulloch, O. Parcollet, and U. Theory, Phys. Rev. B 90, 085102 (2014). Schollwöck, Chebyshev Matrix Product State Impurity [12] D. Zgid, E. Gull, and G. K. L. Chan, Truncated Configu- Solver for Dynamical Mean-Field Theory, Phys. Rev. B ration Interaction Expansions as Solvers for Correlated 90, 115124 (2014). Quantum Impurity Models and Dynamical Mean-Field [31] G. Vidal, Efficient Classical Simulation of Slightly Entangled Theory, Phys. Rev. B 86, 165128 (2012). Quantum Computations, Phys. Rev. Lett. 91, 147902 (2003). [13] S. R. White, Density Matrix Formulation for Quantum [32] G. Vidal, Efficient Simulation of One-Dimensional Quan- Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992). tum Many-Body Systems, Phys. Rev. Lett. 93, 040502 [14] U. Schollwöck, The Density-Matrix Renormalization (2004). [33] S. R. White and A. E. Feiguin, Real-Time Evolution Using Group in the Age of Matrix Product States, Ann. Phys. (Amsterdam) 326, 96 (2011). the Density Matrix Renormalization Group, Phys. Rev. Lett. [15] K. M. Stadler, Z. P. Yin, J. von Delft, G. Kotliar, and A. 93, 076401 (2004). Weichselbaum, Dynamical Mean-Field Theory Plus [34] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Numerical Renormalization-Group Study of Spin-Orbital Time-Dependent Density-Matrix Renormalization-Group Using Adaptive Effective Hilbert Spaces, J. Stat. Mech. Separation in a Three-Band Hund Metal, Phys. Rev. Lett. 115, 136401 (2015). (2004) P04005. [16] A. Horvat, R. Žitko, and J. Mravlje, Low-Energy Physics of [35] M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Three-Orbital Impurity Model with Kanamori Interaction, Held, and F. Verstraete, Efficient DMFT Impurity Solver Phys. Rev. B 94, 165140 (2016). Using Real-Time Dynamics with Matrix Product States, [17] K. M. Stadler, A. K. Mitchell, J. von Delft, and A. Phys. Rev. B 92, 155132 (2015). [36] F. A. Wolf, A. Go, I. P. McCulloch, A. J. Millis, and U. Weichselbaum, Interleaved Numerical Renormalization Group as an Efficient Multiband Impurity Solver, Phys. Schollwöck, Imaginary-Time Matrix Product State Impurity Rev. B 93, 235101 (2016). Solver for Dynamical Mean-Field Theory, Phys. Rev. X 5, [18] K. A. Hallberg, Density-Matrix Algorithm for the Calcu- 041032 (2015). lation of Dynamical Properties of Low-Dimensional [37] A. Holzner, A. Weichselbaum, and J. von Delft, Matrix Product State Approach for a Two-Lead Multilevel Ander- Systems, Phys. Rev. B 52, R9827 (1995). [19] D. J. García, K. Hallberg, and M. J. Rozenberg, Dynamical son Impurity Model, Phys. Rev. B 81, 125126 (2010). Mean Field Theory with the Density Matrix Renormaliza- [38] P. Werner, E. Gull, M. Troyer, and A. J. Millis, Spin tion Group, Phys. Rev. Lett. 93, 246403 (2004). Freezing Transition and Non-Fermi-Liquid Self-Energy in [20] T. D. Kühner and S. R. White, Dynamical Correlation a Three-Orbital Model, Phys. Rev. Lett. 101, 166405 (2008). Functions Using the Density Matrix Renormalization [39] L. de’ Medici, J. Mravlje, and A. Georges, Janus-Faced Influence of Hund’s Rule Coupling in Strongly Correlated Group, Phys. Rev. B 60, 335 (1999). [21] E. Jeckelmann, Dynamical Density-Matrix Renormaliza- Materials, Phys. Rev. Lett. 107, 256401 (2011). tion-Group Method, Phys. Rev. B 66, 045114 (2002). [40] A. Georges, L. de’ Medici, and J. Mravlje, Strong Corre- [22] S. Nishimoto and E. Jeckelmann, Density-Matrix Renorm- lations from Hund’s Coupling, Annu. Rev. Condens. Matter alization Group Approach to Quantum Impurity Problems, Phys. 4, 137 (2013). J. Phys. Condens. Matter 16, 613 (2004). [41] A. J. Kim, H. O. Jeschke, P. Werner, and R. Valentí, [23] R. Peters, Spectral Functions for Single- and Multi-impurity J-Freezing and Hund’s Rules in Spin-Orbit-Coupled Multi- Models Using Density Matrix Renormalization Group, orbital Hubbard Models, Phys. Rev. Lett. 118, 086401 Phys. Rev. B 84, 075139 (2011). (2017). 031013-13 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) [42] F. Verstraete and J. I. Cirac, Matrix Product States [59] For stability reasons, a larger broadening of η ¼ 0.01 eV SC Represent Ground States Faithfully, Phys. Rev. B 73, was used in the first two DMFT cycles. 094423 (2006). [60] E. Koch, G. Sangiovanni, and O. Gunnarsson, Sum Rules [43] V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, and Bath Parametrization for Quantum Cluster Theories, Simulating Strongly Correlated Quantum Systems with Tree Phys. Rev. B 78, 115102 (2008). Tensor Networks, Phys. Rev. B 82, 205105 (2010). [61] I. de Vega, U. Schollwöck, and F. A. Wolf, How to [44] F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Solving Discretize a Quantum Bath for Real-Time Evolution, Phys. Nonequilibrium Dynamical Mean-Field Theory Using Rev. B 92, 155126 (2015). Matrix Product States, Phys. Rev. B 90, 235131 (2014). [62] J. G. Liu, D. Wang, and Q. H. Wang, Quantum Impurities [45] H. Guo and S. R. White, Density Matrix Renormalization in Channel Mixing Baths, Phys. Rev. B 93, 035102 Group Algorithms for Y-Junctions, Phys. Rev. B 74, 060401 (2016). (2006). [63] T. Barthel, Precise Evaluation of Thermal Response [46] Y. Y. Shi, L. M. Duan, and G. Vidal, Classical Simulation of Functions by Optimized Density Matrix Renormalization Quantum Many-Body Systems with a Tree Tensor Network, Group Schemes, New J. Phys. 15, 073010 (2013). [64] Indeed, model calculations performed for fillings of N ¼ 2 Phys. Rev. A 74, 022320 (2006). [47] L. Tagliacozzo, G. Evenbly, and G. Vidal, Simulation of and N ¼ 3 electrons (the latter in the insulating phase) show Two-Dimensional Quantum Systems Using a Tree Tensor that these calculations are of comparable computational Network that Exploits the Entropic Area Law, Phys. Rev. B effort. For any N, we do expect increased numerical effort 80, 235127 (2009). close to the Mott transition though. [48] I. Pižorn, F. Verstraete, and R. M. Konik, Tree Tensor [65] P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Networks and Entanglement Spectra, Phys. Rev. B 88, Luitz, WIEN2k, An Augmented Plane Wave +Local Orbitals 195102 (2013). Program for Calculating Crystal Properties (Techn. [49] P. Schmitteckert, Nonequilibrium Electron Transport Using Universitat Wien, Austria, 2001). the Density Matrix Renormalization Group Method, Phys. [66] M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola, M. Zingl, Rev. B 70, 121302 (2004). O. E. Peil, X. Deng, J. Mravlje, G. J. Kraberger, C. Martins, [50] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. M. Ferrero, and O. Parcollet, TRIQS/DFTTools: A TRIQS Verschelde, and F. Verstraete, Time-Dependent Variational Application for Ab Initio Calculations of Correlated Principle for Quantum Lattices, Phys. Rev. Lett. 107, Materials, Comput. Phys. Commun. 204, 200 (2016). 070601 (2011). [67] M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. [51] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and Parcollet, T. Miyake, A. Georges, and S. Biermann, F. Verstraete, Unifying Time Evolution and Optimization with Dynamical Mean-Field Theory within an Augmented Matrix Product States, Phys. Rev. B 94, 165116 (2016). Plane-Wave Framework: Assessing Electronic Correlations [52] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and in the Iron Pnictide LaFeAsO, Phys. Rev. B 80, 085101 F. Pollmann, Time-Evolving a Matrix Product State with (2009). Long-Ranged Interactions, Phys. Rev. B 91, 165112 [68] M. Aichhorn, L. Pourovskii, and A. Georges, Importance of (2015). Electronic Correlations for Structural and Magnetic Prop- [53] M. Suzuki, Fractal Decomposition of Exponential erties of the Iron Pnictide Superconductor LaFeAsO, Phys. Operators with Applications to Many-Body Theories Rev. B 84, 054529 (2011). and Monte Carlo Simulations, Phys. Lett. A 146, 319 [69] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. (1990). Krivenko, L. Messio, and P. Seth, TRIQS: A Toolbox for [54] Note that we use the term MPO a bit loosely here. What we Research on Interacting Quantum Systems, Comput. Phys. are referring to is an operator factorized in the same forklike Commun. 196, 398 (2015). structure as the state in Fig. 2. [70] A. Liebsch, Surface versus Bulk Coulomb Correlations in [55] E. M. Stoudenmire and S. R. White, Minimally Entangled Photoemission Spectra of SrVO and CaVO , Phys. Rev. 3 3 Typical Thermal State Algorithms, New J. Phys. 12, 055026 Lett. 90, 096401 (2003). (2010). [71] A. Sekiyama, H. Fujiwara, S. Imada, S. Suga, H. Eisaki, S. I. [56] A. Georges, Strongly Correlated Electron Materials: Uchida, K. Takegahara, H. Harima, Y. Saitoh, I. A. Dynamical Mean-Field Theory and Electronic Structure, Nekrasov et al., Mutual Experimental and Theoretical in American Institute of Physics Conference Series, edited Validation of Bulk Photoemission Spectra of by A. Avella and F. Mancini (AIP, New York, 2004), Sr Ca VO , Phys. Rev. Lett. 93, 156402 (2004). 1−x x 3 Vol. 715, pp. 3–74. [72] We note that high-energy peaks already appear in the first [57] F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. DMFT iteration, for which the bath does not have any Posternak, A. Yamasaki, and O. K. Andersen, Dynamical spectral weight at high energies. Mean-Field Theory Using Wannier Functions: A Flexible [73] G. Sangiovanni, A. Toschi, E. Koch, K. Held, M. Capone, Route to Electronic Structure Calculations of Strongly C. Castellani, O. Gunnarsson, S.-K. Mo, J. W. Allen, H.-D. Correlated Materials, Phys. Rev. B 74, 125120 (2006). Kim, A. Sekiyama, A. Yamasaki, S. Suga, and P. Metcalf, [58] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Static versus Dynamical Mean-Field Theory of Mott Parcollet, and C. A. Marianetti, Electronic Structure Antiferromagnets, Phys. Rev. B 73, 205121 (2006). Calculations with Dynamical Mean-Field Theory, Rev. [74] E. Gull, D. R. Reichman, and A. J. Millis, Bold-Line Mod. Phys. 78, 865 (2006). Diagrammatic Monte Carlo Method: General Formulation 031013-14 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) and Application to Expansion around the Noncrossing [78] L. Vaugier, H. Jiang, and S. Biermann, Hubbard U Approximation, Phys. Rev. B 82, 075109 (2010). and Hund Exchange J in Transition Metal Oxides: Screen- [75] P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, TRIQS/ ing versus Localization Trends from Constrained Random CTHYB: A Continuous-Time Quantum Monte Carlo Phase Approximation, Phys. Rev. B 86, 165105 (2012). Hybridisation Expansion Solver for Quantum Impurity [79] Nevertheless, the effective parameters J differ a little from Problems, Comput. Phys. Commun. 200, 274 (2016). those obtained from the main Hubbard band. [76] P. Werner and A. J. Millis, Hybridization Expansion Impu- [80] Note that the peak highest in energy has an atomic energy of rity Solver: General Formulation and Application to Kondo E ¼ U þ 2J þ 2ϵ . Therefore, U can only be determined Lattice and Two-Orbital Models, Phys. Rev. B 74, 155107 after J is found. (2006). [81] ITensor library, http://itensor.org. [77] D. Bergeron and A.-M. S. Tremblay, Algorithms for [82] The numbers 35 and 22 correspond to the impurity links, 9 Optimized Maximum Entropy and Diagnostic Tools for is the bond dimension to the first bath site, and 2 is the Analytic Continuation, Phys. Rev. E 94, 023303 (2016). physical bond dimension. 031013-15 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Physical Review X American Physical Society (APS)

Fork Tensor-Product States: Efficient Multiorbital Real-Time DMFT Solver

Free
15 pages

Loading next page...
 
/lp/aps_physical/fork-tensor-product-states-efficient-multiorbital-real-time-dmft-FJE8qC4Q2l
Publisher
American Physical Society (APS)
Copyright
Copyright © Published by the American Physical Society
eISSN
2160-3308
D.O.I.
10.1103/PhysRevX.7.031013
Publisher site
See Article on Publisher Site

Abstract

PHYSICAL REVIEW X 7, 031013 (2017) 1,* 1 1 1 1,2 Daniel Bauernfeind, Manuel Zingl, Robert Triebl, Markus Aichhorn, and Hans Gerd Evertz Institute of Theoretical and Computational Physics, Graz University of Technology, 8010 Graz, Austria Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA (Received 23 December 2016; revised manuscript received 14 April 2017; published 20 July 2017) 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 high spectral resolution at all frequencies. We use a large number (Oð100Þ) 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 three-orbital test-bed material SrVO . 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. Finally, we show that our solver can be applied even to models with five orbitals. This impurity solver offers a new route to the calculation of precise real-frequency spectral functions of correlated materials. DOI: 10.1103/PhysRevX.7.031013 Subject Areas: Condensed Matter Physics I. INTRODUCTION real-frequency spectra, which therefore become broadened, especially at high energies. ED directly provides spectra on Strongly correlated systems are among the most fasci- the real axis, but it is severely limited in the size of the nating objects that solid-state physics has to offer. The Hilbert space, i.e., in the number of bath sites. Quite interactions between constituents of such systems lead to recently, NRG was shown to be a viable three-band solver emergent phenomena that cannot be deduced from the by exploiting non-Abelian quantum number conservation properties of noninteracting particles [1]. [15–17]. NRG works on the real axis and captures the low- One of the most widely used methods to describe energy physics well, but by construction, it has a poor strongly correlated electrons is the dynamical mean-field resolution at higher energies. Another interesting route that theory (DMFT) [2,3]. DMFT treats local electronic corre- has been proposed recently is the one of solvers that tackle lations using a self-consistent mapping of the lattice the problem of exponential growth of the Hilbert space problem onto an effective Anderson impurity model using ideas from quantum chemistry, i.e., the configuration (AIM). Calculating the single-particle spectral function interaction [11,12]. They allow one to go beyond the small of this impurity model in an accurate and efficient way bath sizes of ED, keeping all the advantages, such as the is at the heart of every DMFT calculation. To this end, absence of fermionic sign problems. However, in multi- many numerical methods have been developed or adapted. orbital applications (see the Appendix of Ref. [12]), the These are based, for instance, on continuous-time quantum spectral resolution has so far been restricted by the number Monte Carlo (CTQMC) [4,5], exact diagonalization (ED) of bath sites [Oð20Þ]. [6–8], the numerical renormalization group (NRG) [9,10], Finally, MPS-based techniques like DMRG do not suffer the configuration interaction (CI) [11,12], and also the from a sign problem and can be used on the real-frequency density-matrix renormalization group (DMRG) with axis as well as on the imaginary-frequency axis. The absence matrix-product states (MPS) [13,14]. of the sign problem, in general, comes at the cost of a very Every algorithm has strengths and weaknesses: CTQMC large growth of bond dimension with the number of orbitals. is exact, apart from statistical errors on the imaginary axis, Dynamical properties and spectral functions can be and can deal with multiple orbitals, but it is, in some cases, calculated within DMRG and have been used for impurity plagued by the fermionic sign problem. Additionally, solvers, e.g., with the Lanczos-like continued-fraction an ill-posed analytic continuation is necessary to obtain expansion [18,19]. Other solvers using the more stable correction vector [20] and dynamical DMRG (DDMRG) [21] methods were developed [22–25]. Both algorithms daniel.bauernfeind@tugraz.at produce very accurate spectral functions but have the Published by the American Physical Society under the terms of disadvantage that a separate calculation for each frequency the Creative Commons Attribution 4.0 International license. has to be performed. The Chebyshev expansion [26] with Further distribution of this work must maintain attribution to MPS [27], supplemented by linear prediction [28],was the author(s) and the published article’s title, journal citation, and DOI. used for impurity solvers in the single-band case [29] and 2160-3308=17=7(3)=031013(15) 031013-1 Published by the American Physical Society DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) for two bands [30]. Recently, some of us introduced a real-time calculations. We emphasize that (i) our method method based on real-time evolution [31–34] and achieved is, by construction, free of any fermionic sign problem, a self-consistent DMFT solution for a two-band model (ii) one can fully converge the DMFT self-consistency loop [35]. In such calculations, the physical orbitals for each on the real-frequency axis, and (iii) we can achieve an almost spin direction are usually combined into one large site in exact representation of the bath spectral function. We apply the MPS. Three or more orbitals have not been feasible with this method to multi-orbital DMFT for the test-bed material this approach because of a large increase in computational SrVO and show that one can resolve a multiplet structure in cost with the number of orbitals. Another MPS-based the Hubbard bands, at the same time keeping a good solver, which works on the imaginary axis, was recently description of the low-energy quasiparticle excitations. introduced [36], and it was applied as a solver for three The paper is structured as follows. First, we outline how bands in two-site cluster DMFT. It was supplemented by a impurity solvers with tensor networks work and introduce single real-time evolution to compute the spectral function, our new tensor network approach, which we call fork avoiding the analytic continuation. However, this method is tensor-product states (FTPS) (Sec. II). Next, we explain in restricted by the number of bath sites that can be employed. detail how our solver is used in the context of multi-orbital In the calculation mentioned, only three bath sites per DMFT (Sec. III). In Sec. IV, we apply our approach to orbital were used, limiting the energy resolution for real- SrVO and discuss the multiplet structure that the FTPS frequency spectral functions. solver allows us to resolve. In order to check the accuracy In the present paper, we introduce a novel impurity of the method, we also compare the FTPS results to solver that works directly on the real-frequency axis. To this CTQMC for SrVO . Finally, we show the efficiency of end, we use a tensor network that captures the geometry of the FTPS solver by applying it to a five-orbital model. the interactions in the Anderson model better than a standard MPS. Our approach is, to some extent, inspired II. TENSOR NETWORK IMPURITY SOLVERS by the work of Ref. [37], which used a similar network for a The AIM describes an impurity (with Hamiltonian H ) two-orbital NRG ground-state calculation. We are not loc coupled to a bath of noninteracting fermions hybridized restricted to a small number of bath sites. This is imperative for exploiting the spectral resolution achievable with with it. A typical AIM Hamiltonian is given by H ¼ H þ H ; loc bath H ¼ ϵ n þ H þ H ; loc 0 m0σ DD SF-PH mσ X X X 0 0 H ¼ U n n þðU − 2JÞ n n þðU − 3JÞ n n ; DD m0↑ m0↓ m0σ m 0¯σ m0σ m 0σ 0 0 m >m;σ m >m;σ X X † † † † 0 0 0 H ¼ J ðc c c c þ H:c:Þ − J ðc c c c þ H:c:Þ; SF-PH m0↓ m 0↑ m 0↑ m 0↓ m0↑ m 0↓ m0↑ m0↓ 0 0 m >m m >m H ¼ ϵ n þ V ðc c þ H:c:Þ; ð1Þ bath l mlσ l mlσ m0σ mlσ where c (c ) creates (annihilates) an electron in band GðtÞ¼ −iθðtÞhψ j½c ðtÞ;cð0Þjψ i; ð2Þ mlσ mlσ 0 0 m (m ∈ f1; 2; 3g for a three-orbital model) with spin σ at the lth site of the system (the impurity has index l ¼ 0; the bath degrees of freedom have l ≥ 1), and n are the mlσ of the interacting problem (1), either in real- or imaginary- corresponding particle-number operators. Note that H DD time t. In the present paper, we introduce a new tensor describes density-density (DD) interactions between all network similar to a MPS, which can be used as a real-time orbitals, and H are the spin-flip and pair-hopping SF-PH impurity solver for three orbitals. We first introduce MPS terms. This three-orbital Hamiltonian is not only important before moving on to what we call fork tensor-product states in the context of real-material calculations. It has also in Sec. II B. been studied extensively on the model level—most im- portantly, because it hosts unconventional correlation phe- A. MPS and DMRG nomena. For a selection of recent work, see, for instance, Refs. [15,38–41]. MPS are a powerful tool to efficiently encode quantum- An impurity solver calculates the retarded impurity mechanical states. Consider a state jψi of a system Green’s function GðtÞ, consisting of N sites: 031013-2 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) state jψ i and ground-state energy E . It minimizes the jψi¼ c js ;s ;…;s i: ð3Þ 0 0 s ;…;s 1 2 N 1 N s ;s ;…;s expectation value 1 2 N Each site i has a local Hilbert space of dimension d hψjHjψi E ¼ min ; ð6Þ spanned by the states js i. Through repeated use of i jψi hψjψi singular-value decompositions (SVDs), it is possible to factorize every coefficient c into a product of usually by updating two neighboring MPS tensors before s ;…;s 1 N matrices [14], i.e., into a MPS, moving on to the next bond. This procedure also yields the Schmidt decomposition of the state at the current bond on s s s 1 2 N jψi¼ A · A  A js ;s ;…;s i: ð4Þ the fly. The DMRG approximation keeps only those states 1 2 N 1 2 N s ;s ;…;s with the largest Schmidt coefficient. It is important to note 1 2 N that one can perform a DMRG calculation for any tensor Each A is a rank-three tensor, except the first and last ones network, as long as one can generate a Schmidt decom- s s 1 N (A , A ), which are of rank two. The index s is called a i position [43]. 1 N physical index, and the matrix indices, which are summed For obtaining the Green’s function, we employ an over, are the so-called bond indices. A general state of the evolution in real-time. Equation (2) is split into the greater > < full Hilbert space is not feasible to store, but it can be shown (G ) and lesser Green’s function (G ): that ground states are well described by a MPS with limited > < bond dimension m (dimension of the bond index) [42]. GðtÞ¼ −iΘðtÞðG ðtÞþ G ðtÞÞ; In complete analogy to the states, one can factorize an > −iHt † iE t G ðtÞ¼hψ jce c jψ ie ; 0 0 operator into what is called a matrix-product operator < † iHt −iE t (MPO) [14], 0 G ðtÞ¼hψ jc e cjψ ie ; ð7Þ 0 0 0 0 s ;s s ;s 1 N 1 N 0 0 which are calculated in two separate time evolutions. This H ¼ W  W js ;…;s ihs ;…;s j; ð5Þ 1 N 1 N 1 N s ;…;s 1 N is done by first applying c (or c) and then time evolving 0 0 s ;…;s 1 N this state and calculating the overlap with the state at time s ;s t ¼ 0. The time evolution is the most computationally where each W is a rank-four tensor. Tensor networks, in expensive part since time-evolved states are not ground general, have a very useful graphical representation, which states anymore, and the needed bond dimensions usually is shown for a MPS and a MPO in Fig. 1. Note that when grow very fast with time. we use the term MPS, we always mean a one-dimensional chain of tensors, as shown in Fig. 1. B. Fork tensor-product states To calculate Green’s functions within the MPS formal- ism, one usually first applies the DMRG [13,14], which So far, the usual way of dealing with Hamiltonians like acts on the space of MPS and finds a variational ground 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, respectively. The local state space of each bath site then (a) consists of M spinless-fermion degrees of freedom, with dimension 2 , 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 encoun- ter this very fast growth upon calculating the ground state (b) 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 FIG. 1. (a) Graphical representation of a MPS. Every circle usually transformed into a Wilson chain with nearest- corresponds to a tensor A and each line to an index of this tensor. neighbor hoppings only, i.e., of the form t ðc c þ i iþ1 In this picture, the physical indices are the vertical lines, while the H:c:Þ [10]. This was thought to be necessary since long- horizontal lines show the bond indices. Connected lines mean range interactions look problematic for MPS-based algo- that the corresponding index is summed over. Fixing all the rithms. Quite recently, though, it was discovered that MPS physical indices s for each site results in a tensor of rank zero can deal with the original form of H in Eq. (1) better bath with the value of the coefficient c . (b) Graphical repre- s ;…;s 1 N [44]. Because all hopping terms in H originate from the bath sentation of a MPO. The difference from a MPS is that a MPO has impurity, this is called the star geometry. The reason for the incoming indices s and outgoing indices s corresponding to the bra- and ket vectors of the operator. better performance is that in the star geometry, one has 031013-3 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) many nearly fully occupied (empty) bath sites with very tensors representing the impurities. When all bond indices low (high) on-site energies ϵ . have the same dimension χ, it is necessary to do a SVD for 2 4 Since basis states with many unoccupied low-energy sites a χ d × χ matrix with computational complexity Oðχ dÞ. have a very low Schmidt coefficient, these states are However, as we show below, this operation does not pose a discarded from the MPS. The same holds for occupied substantial problem for calculations using FTPS. Since the high-energy sites. However, when dealing with multi-orbital impurity tensors pose the biggest challenge, our tensor models, the star geometry is not enough to be able to network would likely also allow us to deal with the chain calculate Green’s functions using MPS. The growth of the geometry without a drastic increase in computational cost. In the present paper, we only use FTPS with baths in star bond dimensions still makes those calculations unfeasible. geometry. The proposed FTPS are similar to the tensor The key idea of the present work is to construct a tensor network used by Holzner et al. [37] to perform NRG network which is beyond a standard MPS but similar calculations for ground-state properties of an AIM with two enough to be able to use established methods like DMRG orbitals. and time evolution. From Hamiltonian (1), one can The three-legged tensors in our network (Fig. 2) can also immediately notice that there are no terms coupling bath be interpreted as two coupled junctions with three legs in sites of different orbitals. Hence, it might not be advanta- the language of Ref. [45], where it has been shown that geous to combine those (not directly interacting) degrees of DMRG is possible on such junctions. Furthermore, our freedom into one large physical index in the MPS. approach has similarities to the so-called tree tensor net- Therefore, our proposed tensor network separates the works (TTN) [43,46–48]. 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 1. Time evolution network is not a MPS anymore since there are some tensors Time evolution with the Hamiltonian Eq. (1) is not (labeled A↓ and B↑ in the example of Fig. 2) that have straightforward since it features long-range hoppings. three bond indices and one physical index, i.e., which are of Possible methods include Krylov approaches [49], the rank four. Cutting any bond splits the network into two time-dependent variational principle [50,51], and the series separate parts. Therefore, one can calculate the Schmidt iHt expansion of e proposed by Zaletel et al. [52]. In this decomposition in a way very similar to a MPS, which work, however, we use a much simpler approach. means that DMRG is also possible. The main bottleneck of First, we split the Hamiltonian into the following terms: calculations with FTPS is to perform SVDs of the rank-four SF-PH (i) the spin-flip and pair-hopping terms H for each P m;m SF-PH orbital combination, with 0 H ¼ H [see SF-PH m >m m;m Eq. (1)], (ii) the density-density interaction terms H , DD and (iii) all other terms, H ¼ H þ ϵ n . With free bath 0 m0σ mσ these terms, we write the time-evolution operator for a small time step Δt using a second-order Suzuki-Trotter decomposition [53], SF-PH −i½ðΔtÞ=2H −iΔtH 0 −i½ðΔtÞ=2H m;m DD e ≈ e e · m >m SF-PH −i½ðΔtÞ=2H −iΔtH −i½ðΔtÞ=2H 0 free DD m;m ×e e e : ð8Þ m >m Note that in this decomposition, the order of the spin-flip and pair-hopping terms is important. The order of operators in FIG. 2. Graphical representation of a FTPS for multi-orbital AIM. Separating bath degrees of freedom leads to a forklike the second product must be opposite to the one in the first. structure. In this picture, a two-orbital model is shown, with four We see that Eq. (8) involves three different operators, SF-PH bath sites in each orbital. Orbitals are labeled A and B, and the H , H , and H , each of which will be treated DD free m;m arrows denote the spin. Each spin-orbital combination has its own differently. bath sticking out to the right. As in Fig. 1, the vertical lines are the Time evolution of the density-density interactions is physical degrees of freedom [all of dimension two, for empty performed with a MPO-like representation of the time- (respectively, occupied) bath sites]. All other lines are bond −i½ðΔtÞ=2H DD evolution operator e . For a three-orbital model, indices, and like in the MPS, they are summed over. As 3 3 −i½ðΔtÞ=2H DD first the full matrix (4 × 4 )of e is created, mentioned in the text, the bath is represented in star geometry which is then decomposed into MPO form by repeated because of the smaller bond dimensions needed. The bath sites SVDs. Since H only consists of density-density inter- are ordered according to their on-site energies. Two example DD −iΔtH DD hoppings V and V are drawn. actions, no fermionic sign appears in e . 1 2 031013-4 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) Time evolution of the spin-flip and pair-hopping terms is SVD(toseparatethe tensorsagain),the swapgateisappliedso more involved than the density-density interactions since that the impurity and the first bath sites are swapped. By the operators change the particle numbers on the impurity repeating this process, one moves the impurity along its sites. Therefore, it can be difficult to deal with the fermionic horizontal arm in Fig. 2. Because a second-order decompo- sign of the time-evolution operator when the impurities are sitionisused,nowalltime-evolutiongatesexcepttheoneatsite not next to each other in the fermionic order. It turns out N have to be applied again. But now, the impurity and bath that the spin-flip and the pair-hopping terms have the sitesneedtobeswapped before time evolution. ˆ ˆ ˆ Note that not only can the algorithm presented above be property A ¼ A individually, with A being either the spin- flip or the pair-hopping operator, respectively. Furthermore, used to perform real-time evolutions, but it is also appli- they commute with each other, allowing us to separate them cable to evolution in imaginary-time simply by replacing without Trotter error. The time-evolution operator of JA is idt by dτ. then given by ˆ III. MULTI-ORBITAL DMFT WITH FTPS −iΔtJA ˆ ˆ e ¼ 1 þ A ðcos ðΔtJÞ − 1Þ − iA sin ðΔtJÞ: ð9Þ In this section, we present details of our impurity solver. For this operator, a MPO can be found for which the We refer to Refs. [2,56] for DMFT in general, and to fermionic sign can easily be determined [54]. Refs. [57,58] for DMFT in the context of realistic ab initio To time evolve the bath terms, we use an iterative second- calculations for correlated materials. order Suzuki-Trotter breakup for each term in H . free In the latter approach, called density-functional theory Neglecting orbital (m) and spin (σ) indices, the first (DFT)+DMFT, the correlated subspace is described by a −iΔt H l¼1 Hubbard-like Hamiltonian. Within DMFT, this model is step in this breakup is the following: e ≈ −iΔt H mapped onto the AIM Hamiltonian (1). This mapping −i½ðΔtÞ=2H −i½ðΔtÞ=2H 1 l¼2 1 e e e . Next, we split off H defines the bath hybridization function ΔðωÞ describing the and iterate this procedure until we end up with influence of the surrounding electrons. N −1 b Since FTPS provide the Green’s function of the AIM on Y Y Δt −iΔtH −i H free mlσ e ≈ e the real-frequency axis, the self-consistency loop is also mσ l¼1 performed directly for real frequencies. For calculating the bath hybridization, we use retarded Green’s functions with Δt −iΔtH mN σ −i H mlσ b 2 × e e ; ð10Þ a finite broadening η in order to avoid numerical SC l¼N −1 difficulties with the poles of the Green’s function. Throughout this work, we use η ¼ 0.005 eV [59]. SC with N thenumberofbathsitesand H ¼ b mlσ The impurity solver calculates the self-energy ΣðωÞ of ϵ n þ V ðc c þ H:c:Þ. In the above equation, we l mlσ l mlσ m0σ the AIM, given the bath hybridization function ΔðωÞ and neglected the term ϵ n that we add to H . 0 m0σ m1σ the interaction Hamiltonian on the impurity. To this end, Equation (10) is a product of two-site gates (an operator acting our solver performs the following steps, which are nontrivially only on two sites) with one of the two sites always explained in more detail in the text below: being the impurity. This means that those two sites are not (1) Obtain bath parameters ϵ and V by a deterministic l l nearest neighbors in the tensor network. To overcome this approach based on integration of the bath hybridi- problem, we use so-called swap gates [14,55]. The two-site zation function ΔðωÞ. operator (2) Calculate the ground state jψ i and ground-state energy E of the interacting problem. n n i j 0 0 0 S ¼ δ δ ð−1Þ ð11Þ ij s ;s s ;s i j j i (3) Apply impurity creation or annihilation operators, and time evolve these states to determine the swaps the state of site i (s with occupation n ) with the i i interacting Green’s function [Eq. (2)]. n n i j state of sitej (s with occupationn ).The factor ð−1Þ gives j j (4) Fourier transform Eq. (2) to obtain GðωÞ and the correct fermionic sign and is negative if an odd number of calculate the local self-energy, particlesonsiteigetsswappedwithanoddnumberofparticles −1 −1 on site j. To be more precise, the matrix representation of the ΣðωÞ¼ G ðωÞ − GðωÞ : ð13Þ swap gates used in this work is To perform step 1 we use S ¼j00ih00jþj10ih01 þj01ih10j − j11ih11j: ð12Þ ij h i V ¼ − ImΔðωÞ dω; Itturnsoutthateveryswapgatecanbecombinedwithanactual l time-evolution gate without additional computational time. h i 1 1 For example, the first step in this time evolution would be to ϵ ¼ ω − ImΔðωÞ dω; ð14Þ −i½ðΔtÞ=2H V π m1σ l l apply e . Immediately afterwards, even before the 031013-5 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) similar to Refs. [10] (NRG) and [44]. Each interval I parameter sets. The next step of our testing was to include corresponds to a bath site. This discretization can be density-density interactions, one term at a time. For interpreted as representing each interval I as a delta peak example, we only included ðU − JÞn n and compared l 10↑ 30↑ at position ϵ and weight V . Sum rules for such discre- energy and Green’s functions to a standard one-orbital MPS tization parameters can be found analytically [60]. In this solver. Finally, we also compared our method to the MPS two-band solver used in Ref. [35]. Indeed, all tests work, we choose the length of each interval such that the performed produced the correct energies and Green’s area of the bath spectral function −ð1=πÞImΔðωÞ is approximately constant for each interval [61]. For the case functions. at hand, this discretization was found to be numerically more stable than using intervals of constant length. Unless IV. RESULTS stated otherwise, we use N ¼ 109 bath sites per orbital We performed DMFT calculations based on a band and spin. We note that this scheme is not restricted to structure obtained from DFT for the prototypical com- diagonal hybridizations. In the general case of off-diagonal pound SrVO , using the approximation of the Kanamori hybridizations, the hybridization function is a matrix Δ. Hamiltonian [Eq. (1)]. It has a cubic crystal structure with a Therefore, instead of taking the imaginary part, we can use nominal filling of one electron in the V-3d shell [64]. the bath spectral function ½i=ð2πÞðΔ − Δ Þ. Similarly to Because of the crystal symmetry, the five orbitals of the Eq. (14), we represent each interval by one delta peak for V-3d shell split into two e and three t orbitals. The latter each orbital. For instance, fixing ϵ to the center of g 2g form the correlated subspace. We performed the DFT the interval, the hopping parameters V can be found calculation with Wien2k [65] and used 34220 k points systematically from the Cholesky factorization of in the irreducible Brillouin zone in order to reach an ½i=ð2πÞðΔ − Δ Þdω. Most importantly, this scheme energy resolution comparable with the η ¼ 0.005 eV does not involve any fitting procedure on the Matsubara SC broadening. axis. A very similar approach was developed independently The TRIQS/DFTTools package (v1.4) [66–68], which is in Ref. [62]. based on the TRIQS library (v1.4) [69], was used to In step 2, we use a DMRG approach with the following generate the projective Wannier functions and to perform parameters, unless specified otherwise. The truncated the DMFT self-consistency cycle. weight t (sum of all discarded singular values of each −8 Figure 3 shows the main results of this paper: the DMFT SVD) is kept smaller than 10 . When spin-flip and pair- spectral function AðωÞ for SrVO , (i) in the approximation hopping terms are neglected, we use an even smaller cutoff −9 of density-density interactions only and (ii) with full of 10 . Note that, except in the five-band calculation, we rotational invariance, including spin-flip and pair-hopping do not restrict the bond dimensions by some hard cutoff terms. Overall, both cases show the well-known features of (see Appendix A2). the SrVO spectral function [70,71]. We see a hole During time evolution (step 3), we use a truncated weight 3 −8 −8 excitation at around −2 eV, as well as the quasiparticle of t ¼ 2 × 10 ,or 10 with density-density interactions −1 peak at zero energy whose shape and position does not only. We time evolve to t ¼ 16 eV , with a time step of −1 Δt ¼ 0.01 eV . Green’s functions are measured every fifth time step. The time-evolution operator of H is loc applied using the zip-up algorithm [55]. Afterwards, the Green’s functions are extrapolated in time using the linear- −1 prediction method [28,35] up to t ¼ 250 eV . 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 iωt−η jtj FT broadening in the kernel e of η ¼ 0.02 eV to FT avoid cutoff effects remaining after the linear prediction. The influence of the linear prediction on our results is discussed in Appendix A1. We stress that although a calculation with full rotational symmetry is more demand- ing, the computational effort is still very feasible. With the parameters mentioned above, one full DMFT cycle takes FIG. 3. Spectral functions AðωÞ for density-density interactions about five hours on 16 cores. (DD) only (blue line), and with spin-flip and pair-hopping To verify that our implementation of DMRG and time terms included (red line). In both calculations, we used U ¼ evolution produces correct results when used with our 4.0 eV and J ¼ 0.6 eV. Both spectra show a three-peak structure tensor network, we first compared Green’s functions and in the upper Hubbard band and additional features at high ground-state energies for U ¼ J ¼ 0 for several bath energies (around 8 eV). 031013-6 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) depend on the inclusion of full rotational invariance. In the (a) upper Hubbard band, a distinctive three-peak structure can be seen. This structure has not been resolved in other exact methods like CTQMC (because of problems with analytic continuation, see below) or NRG (because of the loga- rithmic discretization problem). In our real-time approach, high energies correspond to short times, where the calcu- lations are particularly precise [72]. Most methods allow us 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 noncrossing approximation (NCA) show atomiclike features, but they have very limited accuracy for the description of the low-energy quasiparticle excitations in the metallic phase [74]. Thus, our FTPS solver combines the best of both worlds, with atomic multiplets at high (b) 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 FIG. 4. (a) We take the bath spectral function ΔðωÞ from the between impurities. In all subsequent plots, we show results DMFT self-consistent solution for N ¼ 109 and represent it from calculations with DD interactions only. b using various numbers of bath sites. It is obvious that N ¼ 9 is too small to represent the bath well. (b) Converged DMFT spectral function using the AIM with different numbers of bath A. Effect of bath size sites. Only the smallest bath shows a noticeable difference. This is mostly due to the fact that, in this case, a higher broadening of In order to achieve a reliable high-resolution spectrum η ¼ 0.1 eV had to be used in the Fourier transform and time on the real-frequency axis, it is imperative to have a good FT −1 evolution was only possible for t ¼ 14 eV . The additional representation of the hybridization function ΔðωÞ in small structure at ω ¼ 0 for N ¼ 59 bath sites is most likely a terms of the bath parameters, for which a sufficient linear-prediction artifact. number of bath sites is needed. Figure 4 shows how well a hybridization function can be represented with our B. Comparison to CTQMC approach [Eq. (14)] using a certain number of bath sites. We see that for N ¼ 9 bath sites (we always denote sites In Fig. 5, we compare the converged spectral function of per orbital), ΔðωÞ can be reconstructed only very our approach (FTPS) with a spectrum obtained from roughly, which in turn gives an incorrect spectral function CTQMC and analytic continuation. In both calculations, (Fig. 4, bottom). To some extent, the difference in the we used the same interaction Hamiltonian with density- spectrum is due to the shorter time evolution and there- density interactions only. The CTQMC calculation was fore a higher broadening, η , which we were forced to performed with the TRIQS CTHYB solver (v1.4) [75,76] FT use. For such a small bath, the finite-size effects from with 3.2 × 10 measurements and at inverse temperature −1 reflections at the bath ends appear much earlier in the β ¼ 200 eV . For the analytic continuation, we applied time evolution. the Ω MaxEnt method [77]. Increasing the number of bath sites to N ¼ 29,we The three-peak structure in the upper Hubbard band is observe that the reconstructed bath spectral function not present in the CTQMC spectrum. We show in an already shows the relevant features of ΔðωÞ. The spectrum example below that even for a Green’s function that does is well converged for the largest bath sizes N ¼ 59 contain these peaks, the analytic continuation does not and N ¼ 109. resolve this structure. 031013-7 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) FIG. 5. DMFT spectral functions AðωÞ from CTQMC þ −1 MaxEnt (blue line) at β ¼ 200 eV , and from FTPS (red line). The FTPS result shows a distinctive three-peak structure in the upper Hubbard band. For another comparison, we consider the imaginary-time FIG. 6. Comparison of imaginary-time Green’s functions GðτÞ Green’s functions GðτÞ in Fig. 6. Apart from the effect of from CTQMC (G , blue line) and FTPS using Eq. (15) CTQMC statistical errors, CTQMC provides an exact self-consistent (G , red squares). The agreement is also equally good at β ¼ AðωÞ solution of DMFT on the imaginary-frequency axis. As −1 −1 100 eV and β ¼ 400 eV (not shown). The difference be- mentioned above, when we use the FTPS solver, we tween the two Green’s functions is shown in the bottom panel. formulate the DMFT self-consistency equations on the Note that on both ends, G is smaller than G . The AðωÞ CTQMC real axis. To obtain an approximate finite-temperature normalization of the spectral function demands that imaginary-time Green’s function from FTPS that we can Gðτ ¼ 0Þþ Gðτ ¼ βÞ¼ −1. The CTQMC data deviate in the −2 compare to the CTQMC result, we need to take the finite order of 10 from this constraint due to statistical noise, while temperature of the CTQMC calculation into account. FTPS gives (by construction) the correct result to a precision of −8 10 . This explains the bigger differences of the Green’s Therefore, we use the FTPS spectrum AðωÞ and assume functions around τ ¼ 0 and τ ¼ β. For better visibility of the that we would obtain the same spectrum for a finite (but −3 τ > 0 data, the value of 9 × 10 at τ ¼ 0 is not shown. high enough) inverse temperature β; we use −ωτ energies. Further benchmarks concerning the linear pre- dω e GðτÞ¼ AðωÞ : ð15Þ −βω diction can be found in Appendix A1. 2π e þ 1 Finally, we show that the ill-posedness of the analytic continuation is the most likely explanation for the missing 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 Green’s function GðτÞ and Fourier transform it to get Gðiω Þ: iω τ Gðiω Þ¼ e GðτÞdτ: Looking at the last few DMFT cycles, we estimate it to be around Aðω ¼ 0Þ¼ −ð1=πÞlim ℑGðiω Þ ≈ iω →0 n −1 0.272 eV , 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 FIG. 7. Spectral functions from analytically continued −1 Aðω ¼ 0Þ¼ 0.28 eV , with fluctuations of about imaginary-time Green’s functions GðτÞ calculated by CTQMC −1 0.01 eV . This agreement is very good considering that (blue line) and by FTPS (red line). Clearly, the analytic continu- linear prediction has its strongest influence at small ation cannot resolve the peak structure in the upper Hubbard band. 031013-8 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) TABLE I. Relevant states of the atomic problem of Hamiltonian (1) without spin-flip and pair-hopping terms. Type States Energy difference to ground state Degeneracy N ¼ 1, ground state j↑; 0; 0ij↓; 0; 0ij0; ↑; 0i 06 N ¼ 0 j0; 0; 0i −ϵ 1 N ¼ 2, same spin j↑; ↑; 0ih↑; 0; ↑ij0; ↑; ↑i U − 3J þ ϵ 6 N ¼ 2, different spin j↑; ↓; 0ij↑; 0; ↓ij↓; ↑; 0i  U − 2J þ ϵ 6 N ¼ 2, double occupation j↑↓; 0; 0ij0; ↑↓; 0ij0; 0; ↑↓i U þ ϵ 3 N ¼ 3, all spins equal j↑; ↑; ↑ih↓; ↓; ↓i 3U − 9J þ 2ϵ 2 N ¼ 3, one spin different j↑; ↑; ↓ij↑; ↓; ↑ij↓; ↑; ↑i  3U − 7J þ 2ϵ 6 N ¼ 3, double occupation j↑↓; ↑; 0ij↑↓; ↓; 0ij↑↓; 0; ↑i 3U − 5J þ 2ϵ 12 peak structure in the upper Hubbard band of the spectral Table II shows how bare atomic parameters change when function obtained from the CTQMC data. To do so, we take adding a bath, and we see that our qualitative arguments the FTPS spectrum AðωÞ, calculate GðτÞ as described give a correct idea of how parameters are rescaled. Further evidence that the observed three-peak structure is above, and perform the same analytic continuation that we did for the GðτÞ from CTQMC. We added noise of the indeed a result of atomic physics can be seen in Fig. 8. This order of the CTQMC error to the FTPS data. The resulting figure shows a close-up of the upper Hubbard band for spectrum is shown in Fig. 7, and indeed, the peak structure three different values of J. The corresponding effective in the upper Hubbard band vanishes. parameters J are shown in Table II. We observe that J is also rescaled slightly, but the rescaling gets smaller for higher J. Furthermore, increasing J also increases the total C. Discussion of peak structure—Effective width of the Hubbard band, which scales mostly linearly atomic physics with J. At the same time, measuring the quasiparticle In order to understand the peak structure observed in the spectral weight as a function of J at constant U shows that it spectral functions, we take a look at the underlying atomic increases with increasing J, also implying an increasing problem, where, for simplicity, we start with density- critical U for the metal-to-insulator transition [40]. density interactions only. We show that the same arguments Upon a careful inspection of the spectral function in hold for full rotationally invariant interactions. Fig. 5, we observe small peaks at energies around 8 eV. A Table I shows the relevant atomic states and their close-up of this energy region for different values of J is corresponding energies. The atomic model has a hole shown in Fig. 8. The energy difference between the peaks is excitation at energy −ϵ and three single electron excita- close to 2J and can, again, be well explained by atomic tions with energies U þ ϵ , U − 2J þ ϵ and U − 3J þ ϵ 0 0 0 physics, namely, excitations into states with three electrons relative to the ground state. If we measure the energy on the impurity (Table I) [79]. These excitations originate differences between the three peaks of the upper Hubbard from small admixtures of N ¼ 2 states to the ground state. band in our results, we find values of 1.27 eV and 0.69 eV, With atomic physics in mind, let us take a look again at which are close to the atomic energy differences of 1.2 eV the spectrum using full rotational symmetry (Fig. 3). The and 0.6 eV (J ¼ 0.6 eV). We also find the hole excitation at spin-flip and pair-hopping terms only contribute if there −2.0 eV. This indicates that we can describe the positions are two or more particles present. Thus, the quasiparticle of the observed peaks approximately by atomic physics peak and the hole excitation do not change. The atomic ¯ ¯ with effective parameters ϵ , U, and J and widened peaks. N ¼ 2 sector does change, however. Diagonalizing the Furthermore, the heights of the peaks roughly correspond to the degeneracy of the states in the atomic model (see Table I). TABLE II. Atomic parameters and their effective values ob- We can determine U ¼ 5.97 eV (where U ¼ 4.00 eV) tained from the spectral functions shown in Figs. 5 and 8. For J, from the energy difference of the peak highest in energy to the values themselves were obtained from the energy difference the hole excitation. This increase of U compared to U is of the highest peak to the lowest peak, whereas the uncertainty is plausible, considering the following. When coupling the estimated from ω þ 2J and ω − J. M M impurity to the bath, particles have the possibility to avoid Parameter Atomic value (eV) Effective value (eV) each other by jumping into unoccupied sites of the bath. ϵ −0.86 −2.00 This results in a decrease of hn n i. To model this situation 0 ↑ ↓ U 4.00 5.97 using atomic physics, one needs to increase the interaction J 0.50 0.59(6) strength. Finally, it is well known that J is much less J 0.60 0.66(3) affected by the surrounding electrons than U since the latter J 0.70 0.72(2) is screened significantly stronger [78]. 031013-9 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) regime. In Fig. 9, we show results with U ¼ 5.5 eV at (a) constant J ¼ 0.6 eV. We see a shift of the upper Hubbard band to higher energies, but little shift of the hole excitation. Also, the central peak is shifted and gets slimmer since more weight is transferred into the Hubbard bands. Most importantly, as we approach the strongly correlated 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 Hubbard bands smears out and even vanishes. The close-up 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 addi- (b) tional 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 models on the Bethe lattice [35] upon increasing U. We leave further investigation of this feature to future work. (a) FIG. 8. (a) Close-up of the three-peak structure for various values of J. Additionally, we show vertical lines for the J ¼ 0.6 eV spectrum at energies ω (position of the middle peak) and at ω þ 2J and ω − J. We see that the width of the upper M M Hubbard band is close to 3J. (b) Close-up of the small spectral peaks at high energies. These correspond to excitations into the N ¼ 3 sector of the atomic model (see Table I). The height of each peak can be estimated by the degeneracy of the atomic states. Effective parameters J are 0.53 eV (J ¼ 0.5 eV), 0.59 eV (J ¼ 0.6 eV), and 0.68 eV (J ¼ 0.7 eV). They are obtained from the difference between the two peaks that are highest in energy. (b) Hamiltonian, we find eigenstates with three different energies and differences of 3J ¼ 1.8 eV and 2J ¼ 1.2 eV, respectively. Measuring the energy differences in ¯ ¯ Fig. 3, we find 3J ¼ 1.75 eV and 2J ¼ 1.32 eV. Estimating U ¼ 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 effec- tive 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. FIG. 9. (a) Increasing U results in a slimmer central peak and a shift of the upper Hubbard band. Also, the three-peak structure D. Beyond atomic physics gets smeared out. (b) Close-up of the upper Hubbard band. As in The previous section showed that at U ¼ 4.0 eV, the Fig. 8, additional vertical lines are plotted at ω (position of the spectral features in the Hubbard bands can be well middle peak) and at ω þ 2J and ω − J as a rough guide to M M described by atomic physics with effective parameters where the atomic peaks would be located. With the help of these and widened peaks. It is not clear whether this picture is lines, one can discern a three-peak structure again, but it is valid for higher interaction strengths U in the metallic extended by a feature at the low energy side of the Hubbard band. 031013-10 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) It might at first seem counterintuitive that increasing U Green’s function took about 190 hours on the processors makes the physics less atomiclike. Indeed, at very high specified in Appendix A2. We want to stress, though, that interaction strengths, in the insulating regime, the spectrum the resulting spectrum (as well as Fig. 10) was already fully −1 must become atomiclike again. Here, however, we identify converged at t ¼ 12 eV (70 hours). We note that even −1 an intermediate regime where additional structures appear with only one CPU hour (t ¼ 6 eV ), the resulting when increasing U since we get closer to the Mott metal-to- spectrum is almost converged and barely distinguishable insulator transition. from the final result. The benchmark therefore shows that with our FTPS approach, a full five-orbital DMFT calcu- lation is well within reach. E. Solution of a five-band AIM In this section, we show that FTPS not only applies to three-band models, but it also works in the case of five V. CONCLUSIONS orbitals. To do so, we use the bath parameters ϵ and V k k We have presented a novel multi-orbital impurity solver from the converged N ¼ 59 calculation for SrVO and b 3 which uses a forklike tensor network whose geometry construct an artificial degenerate five-band AIM. resembles that of the Hamiltonian. The network structure is Interaction parameters are U ¼ 4.0 eV and J ¼ 0.6 eV. simple enough to generate Schmidt decompositions, We decrease the on-site energy ϵ to get a similar allowing us to truncate the tensor network safely and to occupation of each impurity orbital as for SrVO use established methods like DMRG and real-time evolu- ðhn i ≈ Þ. Note that, in doing so, we have a model m;0;σ tion. The solver works on the real frequency axis and hence with, in total, electrons on the impurity. We only use allows us to formulate the full DMFT self-consistency density-density interactions and carry out the time evolu- procedure for real frequencies. Therefore, results are not −1 tion to t ¼ 16 eV . We set the truncated weight to plagued by an ill-conditioned analytic continuation. Our −8 t ¼ 10 , but we restrict the bond dimension of the approach exhibits no sign problem, though it does become impurity-impurity links to χ ¼ 200. max more involved for larger numbers of orbitals. In Fig. 10, we compare the results obtained for this five- We tested the solver within DMFT on a Hamiltonian band model to results from CTQMC, where we used the typically used for the test-bed material SrVO and inves- same discretized bath hybridization as input to CTQMC. tigated the influence of full rotational invariance on the We again see excellent agreement, even on a logarithmic results. We found clear spectral structures, in particular, in scale. The spectrum AðωÞ (not shown) again exhibits strong the upper Hubbard band that have not been accessible by structure in the upper Hubbard band. Of course, the CTQMC, for which the necessary analytic continuation computational complexity is larger than in the three-orbital prohibits the resolution of fine structures in the spectral case, and it grows during time evolution. Calculating the 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 approach is feasible for five-orbital models, by comparing results from the FTPS solver to CTQMC for an artificial five-band model. ACKNOWLEDGMENTS The authors acknowledge financial support by the Austrian Science Fund (FWF) through SFB ViCoM F41 (P04 and P03), through Project No. P26220, and through the START Program No. 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. FIG. 10. Comparison of imaginary-time Green’s functions GðτÞ Held, F. Maislinger, and G. Kraberger. The computational from CTQMC (G , blue line) and FTPS using Eq. (15) CTQMC resources have been provided by the Vienna Scientific (G , red squares). As in Fig. 6, they compare very well. AðωÞ 031013-11 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) Cluster (VSC). All calculations involving tensor networks were performed using the ITensor library [81]. APPENDIX A: INFLUENCE OF COMPUTATIONAL PARAMETERS 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. A1). Then, we show that the results are converged with respect to the usual MPS approximation (Sec. A2). 1. Linear prediction In order to obtain smooth and sharp spectra, we used FIG. 12. Different truncation values t in the time evolution do linear prediction to extrapolate the Green’s function in time not influence the shape of the spectrum AðωÞ. [28,35]. Without going into detail, we state the fact that linear prediction has two parameters, the pseudo-inverse cutoff p and the order N of the linear prediction. Figure 11 parameter is the only approximation in the representation of inv shows that the results are converged in these parameters. a state as a tensor-product state, as we do not impose any We also show a DMFT run without any linear prediction, hard cutoff on the bond dimensions. Figure 12 shows that which is only possible if we increase the broadening the spectrum is well converged with respect to the trunca- parameter of the Fourier transform to η ¼ 0.1 eV since FT tion error during time evolution. otherwise we would get oscillations due to the hard cutoff Finally, we comment on the required computational of the time series. Except for a shift towards the right, effort. In the calculation without full rotational symmetry, omitting the linear prediction only changes the height (and the size of the largest tensor to represent the ground state −9 width) of the peaks but not the overall structure. This is a was 35 × 22 × 9 × 2 (t ¼ 10 ) [82] and at the end of −8 strong indicator of the stability of these features. time evolution 127 × 79 × 30 × 2 (t ¼ 10 ). For a trun- −7 cated weight of t ¼ 10 , calculating the Green’s function 2. Truncation of the tensor network took about 17 minutes on a node with two processors (Intel Xeon E5-2650v2, 2.6 GHz with eight cores, and G and One of the most important parameters in any MPS-like G each calculated on one processor). This time increases calculation is the sum of discarded singular values in each to five hours for the lowest truncated weight of SVD (truncated weight t ). We want to emphasize that this −9 t ¼ 5 × 10 . Using the full rotationally invariant Hamiltonian, the biggest tensor in the ground-state search −8 was 90 × 40 × 10 × 2 (t ¼ 10 ), and at the end of time −8 evolution, 79 × 46 × 21 × 2 (t ¼ 2 × 10 ). The Green’s function takes about three hours, and we need five hours for one full DMFT cycle on the same two processors as above. [1] P. W. Anderson, More Is Different, Science 177, 393 (1972). [2] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical Mean-Field Theory of Strongly Correlated Fermion Systems and the Limit of Infinite Dimensions, Rev. Mod. Phys. 68, 13 (1996). [3] W. Metzner and D. Vollhardt, Correlated Lattice Fermions FIG. 11. Spectrum AðωÞ using different linear prediction in d ¼ ∞ Dimensions, Phys. Rev. Lett. 62, 324 (1989). parameters for a calculation without spin-flip and pair-hopping [4] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. terms. The calculations with linear prediction were performed Troyer, and P. Werner, Continuous-Time Monte Carlo with a broadening of η ¼ 0.02 eV. Except for small changes Methods for Quantum Impurity Models, Rev. Mod. Phys. FT around ω ¼ 0, the effect of the various linear prediction param- 83, 349 (2011). eters is minor. The blue line lies directly below the red and green [5] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. lines. We also show a DMFT calculation without any linear Millis, Continuous-Time Solver for Quantum Impurity prediction. The main features are still present. Models, Phys. Rev. Lett. 97, 076405 (2006). 031013-12 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) [6] M. Caffarel and W. Krauth, Exact Diagonalization [24] M. Karski, C. Raas, and G. S. Uhrig, Single-Particle Approach to Correlated Fermions in Infinite Dimensions: Dynamics in the Vicinity of the Mott-Hubbard Metal-to- Mott Transition and Superconductivity, Phys. Rev. Lett. 72, Insulator Transition, Phys. Rev. B 77, 075116 (2008). 1545 (1994). [25] M. Karski, C. Raas, and G. S. Uhrig, Electron Spectra Close [7] M. Capone, L. de’ Medici, and A. Georges, Solving the to a Metal-to-Insulator Transition, Phys. Rev. B 72, 113110 Dynamical Mean-Field Theory at Very Low Temperatures (2005). Using the Lanczos Exact Diagonalization, Phys. Rev. B 76, [26] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The 245116 (2007). Kernel Polynomial Method, Rev. Mod. Phys. 78, 275 [8] J. Kolorenč, A. B. Shick, and A. I. Lichtenstein, Electronic (2006). Structure and Core-Level Spectra of Light Actinide Diox- [27] A. Holzner, A. Weichselbaum, I. P. McCulloch, U. Schollwöck, and J. von Delft, Chebyshev Matrix Product ides in the Dynamical Mean-Field Theory, Phys. Rev. B 92, 085125 (2015). State Approach for Spectral Functions, Phys. Rev. B 83, [9] K. G. Wilson, The Renormalization Group: Critical Phe- 195115 (2011). nomena and the Kondo Problem, Rev. Mod. Phys. 47, 773 [28] S. R. White and I. Affleck, Spectral Function for the S ¼ 1 (1975). Heisenberg Antiferromagnetic Chain, Phys. Rev. B 77, [10] R. Bulla, T. A. Costi, and T. Pruschke, Numerical Renorm- 134437 (2008). [29] M. Ganahl, P. Thunström, F. Verstraete, K. Held, and H. G. alization Group Method for Quantum Impurity Systems, Rev. Mod. Phys. 80, 395 (2008). Evertz, Chebyshev Expansion for Impurity Models Using [11] Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort, Matrix Product States, Phys. Rev. B 90, 045144 (2014). Efficient Real-Frequency Solver for Dynamical Mean-Field [30] F. A. Wolf, I. P. McCulloch, O. Parcollet, and U. Theory, Phys. Rev. B 90, 085102 (2014). Schollwöck, Chebyshev Matrix Product State Impurity [12] D. Zgid, E. Gull, and G. K. L. Chan, Truncated Configu- Solver for Dynamical Mean-Field Theory, Phys. Rev. B ration Interaction Expansions as Solvers for Correlated 90, 115124 (2014). Quantum Impurity Models and Dynamical Mean-Field [31] G. Vidal, Efficient Classical Simulation of Slightly Entangled Theory, Phys. Rev. B 86, 165128 (2012). Quantum Computations, Phys. Rev. Lett. 91, 147902 (2003). [13] S. R. White, Density Matrix Formulation for Quantum [32] G. Vidal, Efficient Simulation of One-Dimensional Quan- Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992). tum Many-Body Systems, Phys. Rev. Lett. 93, 040502 [14] U. Schollwöck, The Density-Matrix Renormalization (2004). [33] S. R. White and A. E. Feiguin, Real-Time Evolution Using Group in the Age of Matrix Product States, Ann. Phys. (Amsterdam) 326, 96 (2011). the Density Matrix Renormalization Group, Phys. Rev. Lett. [15] K. M. Stadler, Z. P. Yin, J. von Delft, G. Kotliar, and A. 93, 076401 (2004). Weichselbaum, Dynamical Mean-Field Theory Plus [34] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Numerical Renormalization-Group Study of Spin-Orbital Time-Dependent Density-Matrix Renormalization-Group Using Adaptive Effective Hilbert Spaces, J. Stat. Mech. Separation in a Three-Band Hund Metal, Phys. Rev. Lett. 115, 136401 (2015). (2004) P04005. [16] A. Horvat, R. Žitko, and J. Mravlje, Low-Energy Physics of [35] M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Three-Orbital Impurity Model with Kanamori Interaction, Held, and F. Verstraete, Efficient DMFT Impurity Solver Phys. Rev. B 94, 165140 (2016). Using Real-Time Dynamics with Matrix Product States, [17] K. M. Stadler, A. K. Mitchell, J. von Delft, and A. Phys. Rev. B 92, 155132 (2015). [36] F. A. Wolf, A. Go, I. P. McCulloch, A. J. Millis, and U. Weichselbaum, Interleaved Numerical Renormalization Group as an Efficient Multiband Impurity Solver, Phys. Schollwöck, Imaginary-Time Matrix Product State Impurity Rev. B 93, 235101 (2016). Solver for Dynamical Mean-Field Theory, Phys. Rev. X 5, [18] K. A. Hallberg, Density-Matrix Algorithm for the Calcu- 041032 (2015). lation of Dynamical Properties of Low-Dimensional [37] A. Holzner, A. Weichselbaum, and J. von Delft, Matrix Product State Approach for a Two-Lead Multilevel Ander- Systems, Phys. Rev. B 52, R9827 (1995). [19] D. J. García, K. Hallberg, and M. J. Rozenberg, Dynamical son Impurity Model, Phys. Rev. B 81, 125126 (2010). Mean Field Theory with the Density Matrix Renormaliza- [38] P. Werner, E. Gull, M. Troyer, and A. J. Millis, Spin tion Group, Phys. Rev. Lett. 93, 246403 (2004). Freezing Transition and Non-Fermi-Liquid Self-Energy in [20] T. D. Kühner and S. R. White, Dynamical Correlation a Three-Orbital Model, Phys. Rev. Lett. 101, 166405 (2008). Functions Using the Density Matrix Renormalization [39] L. de’ Medici, J. Mravlje, and A. Georges, Janus-Faced Influence of Hund’s Rule Coupling in Strongly Correlated Group, Phys. Rev. B 60, 335 (1999). [21] E. Jeckelmann, Dynamical Density-Matrix Renormaliza- Materials, Phys. Rev. Lett. 107, 256401 (2011). tion-Group Method, Phys. Rev. B 66, 045114 (2002). [40] A. Georges, L. de’ Medici, and J. Mravlje, Strong Corre- [22] S. Nishimoto and E. Jeckelmann, Density-Matrix Renorm- lations from Hund’s Coupling, Annu. Rev. Condens. Matter alization Group Approach to Quantum Impurity Problems, Phys. 4, 137 (2013). J. Phys. Condens. Matter 16, 613 (2004). [41] A. J. Kim, H. O. Jeschke, P. Werner, and R. Valentí, [23] R. Peters, Spectral Functions for Single- and Multi-impurity J-Freezing and Hund’s Rules in Spin-Orbit-Coupled Multi- Models Using Density Matrix Renormalization Group, orbital Hubbard Models, Phys. Rev. Lett. 118, 086401 Phys. Rev. B 84, 075139 (2011). (2017). 031013-13 DANIEL BAUERNFEIND et al. PHYS. REV. X 7, 031013 (2017) [42] F. Verstraete and J. I. Cirac, Matrix Product States [59] For stability reasons, a larger broadening of η ¼ 0.01 eV SC Represent Ground States Faithfully, Phys. Rev. B 73, was used in the first two DMFT cycles. 094423 (2006). [60] E. Koch, G. Sangiovanni, and O. Gunnarsson, Sum Rules [43] V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, and Bath Parametrization for Quantum Cluster Theories, Simulating Strongly Correlated Quantum Systems with Tree Phys. Rev. B 78, 115102 (2008). Tensor Networks, Phys. Rev. B 82, 205105 (2010). [61] I. de Vega, U. Schollwöck, and F. A. Wolf, How to [44] F. A. Wolf, I. P. McCulloch, and U. Schollwöck, Solving Discretize a Quantum Bath for Real-Time Evolution, Phys. Nonequilibrium Dynamical Mean-Field Theory Using Rev. B 92, 155126 (2015). Matrix Product States, Phys. Rev. B 90, 235131 (2014). [62] J. G. Liu, D. Wang, and Q. H. Wang, Quantum Impurities [45] H. Guo and S. R. White, Density Matrix Renormalization in Channel Mixing Baths, Phys. Rev. B 93, 035102 Group Algorithms for Y-Junctions, Phys. Rev. B 74, 060401 (2016). (2006). [63] T. Barthel, Precise Evaluation of Thermal Response [46] Y. Y. Shi, L. M. Duan, and G. Vidal, Classical Simulation of Functions by Optimized Density Matrix Renormalization Quantum Many-Body Systems with a Tree Tensor Network, Group Schemes, New J. Phys. 15, 073010 (2013). [64] Indeed, model calculations performed for fillings of N ¼ 2 Phys. Rev. A 74, 022320 (2006). [47] L. Tagliacozzo, G. Evenbly, and G. Vidal, Simulation of and N ¼ 3 electrons (the latter in the insulating phase) show Two-Dimensional Quantum Systems Using a Tree Tensor that these calculations are of comparable computational Network that Exploits the Entropic Area Law, Phys. Rev. B effort. For any N, we do expect increased numerical effort 80, 235127 (2009). close to the Mott transition though. [48] I. Pižorn, F. Verstraete, and R. M. Konik, Tree Tensor [65] P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Networks and Entanglement Spectra, Phys. Rev. B 88, Luitz, WIEN2k, An Augmented Plane Wave +Local Orbitals 195102 (2013). Program for Calculating Crystal Properties (Techn. [49] P. Schmitteckert, Nonequilibrium Electron Transport Using Universitat Wien, Austria, 2001). the Density Matrix Renormalization Group Method, Phys. [66] M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola, M. Zingl, Rev. B 70, 121302 (2004). O. E. Peil, X. Deng, J. Mravlje, G. J. Kraberger, C. Martins, [50] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. M. Ferrero, and O. Parcollet, TRIQS/DFTTools: A TRIQS Verschelde, and F. Verstraete, Time-Dependent Variational Application for Ab Initio Calculations of Correlated Principle for Quantum Lattices, Phys. Rev. Lett. 107, Materials, Comput. Phys. Commun. 204, 200 (2016). 070601 (2011). [67] M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. [51] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and Parcollet, T. Miyake, A. Georges, and S. Biermann, F. Verstraete, Unifying Time Evolution and Optimization with Dynamical Mean-Field Theory within an Augmented Matrix Product States, Phys. Rev. B 94, 165116 (2016). Plane-Wave Framework: Assessing Electronic Correlations [52] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and in the Iron Pnictide LaFeAsO, Phys. Rev. B 80, 085101 F. Pollmann, Time-Evolving a Matrix Product State with (2009). Long-Ranged Interactions, Phys. Rev. B 91, 165112 [68] M. Aichhorn, L. Pourovskii, and A. Georges, Importance of (2015). Electronic Correlations for Structural and Magnetic Prop- [53] M. Suzuki, Fractal Decomposition of Exponential erties of the Iron Pnictide Superconductor LaFeAsO, Phys. Operators with Applications to Many-Body Theories Rev. B 84, 054529 (2011). and Monte Carlo Simulations, Phys. Lett. A 146, 319 [69] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. (1990). Krivenko, L. Messio, and P. Seth, TRIQS: A Toolbox for [54] Note that we use the term MPO a bit loosely here. What we Research on Interacting Quantum Systems, Comput. Phys. are referring to is an operator factorized in the same forklike Commun. 196, 398 (2015). structure as the state in Fig. 2. [70] A. Liebsch, Surface versus Bulk Coulomb Correlations in [55] E. M. Stoudenmire and S. R. White, Minimally Entangled Photoemission Spectra of SrVO and CaVO , Phys. Rev. 3 3 Typical Thermal State Algorithms, New J. Phys. 12, 055026 Lett. 90, 096401 (2003). (2010). [71] A. Sekiyama, H. Fujiwara, S. Imada, S. Suga, H. Eisaki, S. I. [56] A. Georges, Strongly Correlated Electron Materials: Uchida, K. Takegahara, H. Harima, Y. Saitoh, I. A. Dynamical Mean-Field Theory and Electronic Structure, Nekrasov et al., Mutual Experimental and Theoretical in American Institute of Physics Conference Series, edited Validation of Bulk Photoemission Spectra of by A. Avella and F. Mancini (AIP, New York, 2004), Sr Ca VO , Phys. Rev. Lett. 93, 156402 (2004). 1−x x 3 Vol. 715, pp. 3–74. [72] We note that high-energy peaks already appear in the first [57] F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. DMFT iteration, for which the bath does not have any Posternak, A. Yamasaki, and O. K. Andersen, Dynamical spectral weight at high energies. Mean-Field Theory Using Wannier Functions: A Flexible [73] G. Sangiovanni, A. Toschi, E. Koch, K. Held, M. Capone, Route to Electronic Structure Calculations of Strongly C. Castellani, O. Gunnarsson, S.-K. Mo, J. W. Allen, H.-D. Correlated Materials, Phys. Rev. B 74, 125120 (2006). Kim, A. Sekiyama, A. Yamasaki, S. Suga, and P. Metcalf, [58] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Static versus Dynamical Mean-Field Theory of Mott Parcollet, and C. A. Marianetti, Electronic Structure Antiferromagnets, Phys. Rev. B 73, 205121 (2006). Calculations with Dynamical Mean-Field Theory, Rev. [74] E. Gull, D. R. Reichman, and A. J. Millis, Bold-Line Mod. Phys. 78, 865 (2006). Diagrammatic Monte Carlo Method: General Formulation 031013-14 FORK TENSOR-PRODUCT STATES: EFFICIENT MULTIORBITAL … PHYS. REV. X 7, 031013 (2017) and Application to Expansion around the Noncrossing [78] L. Vaugier, H. Jiang, and S. Biermann, Hubbard U Approximation, Phys. Rev. B 82, 075109 (2010). and Hund Exchange J in Transition Metal Oxides: Screen- [75] P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, TRIQS/ ing versus Localization Trends from Constrained Random CTHYB: A Continuous-Time Quantum Monte Carlo Phase Approximation, Phys. Rev. B 86, 165105 (2012). Hybridisation Expansion Solver for Quantum Impurity [79] Nevertheless, the effective parameters J differ a little from Problems, Comput. Phys. Commun. 200, 274 (2016). those obtained from the main Hubbard band. [76] P. Werner and A. J. Millis, Hybridization Expansion Impu- [80] Note that the peak highest in energy has an atomic energy of rity Solver: General Formulation and Application to Kondo E ¼ U þ 2J þ 2ϵ . Therefore, U can only be determined Lattice and Two-Orbital Models, Phys. Rev. B 74, 155107 after J is found. (2006). [81] ITensor library, http://itensor.org. [77] D. Bergeron and A.-M. S. Tremblay, Algorithms for [82] The numbers 35 and 22 correspond to the impurity links, 9 Optimized Maximum Entropy and Diagnostic Tools for is the bond dimension to the first bath site, and 2 is the Analytic Continuation, Phys. Rev. E 94, 023303 (2016). physical bond dimension. 031013-15

Journal

Physical Review XAmerican Physical Society (APS)

Published: Jul 1, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off