Generalization of the Time-Dependent Numerical Renormalization Group Method to Finite Temperatures and General Pulses
GGeneralization of the Time-Dependent Numerical Renormalization Group Method to FiniteTemperatures and General Pulses
H. T. M. Nghiem and T. A. Costi Peter Gr¨unberg Institut and Institute for Advanced Simulation, Research Centre J¨ulich, 52425 J¨ulich, Germany (Dated: November 11, 2018)The time-dependent numerical renormalization group (TDNRG) method [Anders et al.
Phys. Rev. Lett. ,196801 (2005)] o ff ers the prospect of investigating in a non-perturbative manner the time-dependence of localobservables of interacting quantum impurity models at all time scales following a quantum quench. Here, wepresent a generalization of this method to arbitrary finite temperature by making use of the full density matrixapproach [Weichselbaum et al. Phys. Rev. Lett. , 076402 (2007)]. We show that all terms in the projectedfull density matrix ρ i → f = ρ ++ + ρ −− + ρ + − + ρ − + appearing in the time-evolution of a local observable maybe evaluated in closed form at finite temperature, with ρ + − = ρ − + =
0. The expression for ρ −− is shown tobe finite at finite temperature, becoming negligible only in the limit of vanishing temperatures. We prove thatthis approach recovers the short-time limit for the expectation value of a local observable exactly at arbitrarytemperatures. In contrast, the corresponding long-time limit is recovered exactly only for a continuous bath, i.e.when the logarithmic discretization parameter Λ → + . Since the numerical renormalization group approachbreaks down in this limit, and calculations have to be carried out at Λ >
1, the long-time behavior followingan arbitrary quantum quench has a finite error, which poses an obstacle for the method, e.g., in its applicationto the scattering states numerical renormalization group method for describing steady state non-equilibriumtransport through correlated impurities [Anders, Phys. Rev. Lett. , 066804, (2008)]. We suggest a way toovercome this problem by noting that the time-dependence, in general, and the long-time limit, in particular,becomes increasingly more accurate on reducing the size of the quantum quench. This suggests an improvedgeneralized TDNRG approach in which the system is time-evolved between the initial and final states via asequence of small quantum quenches within a finite time interval instead of by a single large and instantaneousquantum quench. The formalism for this is provided, thus generalizing the TDNRG method to multiple quantumquenches, periodic switching, and general pulses. This formalism, like our finite temperature generalization ofthe single-quench case, rests on no other approximation than the NRG approximation. The results are illustratednumerically by application to the Anderson impurity model.
PACS numbers: 75.20.Hr, 71.27. + a, 72.15.Qm, 73.63.Kv I. INTRODUCTION
The numerical renormalization group (NRG) has provento be one of the most powerful methods for dealing with equi-librium properties of strongly correlated quantum impuritysystems, allowing the calculation at arbitrary temperatures,and in a non-perturbative manner, of thermodynamic, dynamic, and, linear transport properties of such sys-tems. The method is also applicable, within dynamical meanfield theory (DMFT), to correlated lattice models, suchas the Kondo lattice, Anderson lattice or Hubbard and Hubbard-Holstein models. The introduction of the cor-relation self-energy, the reduced density matrix, the com-plete basis set of eliminated states, novel non-Abeliansymmetries, and new discretization schemes havesignificantly improved NRG calculations, particularly fordynamic and transport properties of single and multi-channel models. A further application of the method, based on the com-plete basis set, in combination with the reduced density ma-trix, is to the time-dependent transient response of quantumimpurity systems following a quantum quench. The use ofthe complete basis set is particularly important here, as it re-solves the problem of summing up multiple shell contribu-tions to transient quantities encountered in a previous relatedapproach. This time-dependent numerical renormalization group (TDNRG) method has been used to study the transientdynamics of a number of quantum impurity models, includ-ing the Anderson, Kondo, interacting resonant level and spin-boson models.
In addition, within the scattering statesNRG approach, the TDNRG approach o ff ers the prospectof investigating another important class of problems, namelytruly non-equilibrium steady state transport through correlatednanostructures, such as through correlated quantum dots at fi-nite bias voltage far from the linear response regime. Here,the use of the TDNRG to evolve the density matrix to longtimes, a key ingredient in the approach of Ref. 45, is ham-pered by di ffi culties in obtaining the correct long-time limit ofobservables, including that of the density matrix, e.g., for theAnderson impurity model. In addition, there is signifi-cant noise at intermediate to long times.
It has been ar-gued that the logarithmically discretized Wilson chain, whoseheat capacity is non-extensive for discretization parameters Λ >
1, prohibits thermalization to the correct long-time limitand possibly a ff ects also the short-time limit. This wouldpose a severe limitation on the TDNRG and its application tosteady-state non-equilibrium situations. Reflections, associ-ated with the non-constant hoppings along the Wilson chain,have also been argued to a ff ect the long-time limit. Thishas lead to the formulation of hybrid methods which com-bine the merits of the TDNRG method at short times, with, forexample, methods, such as the time-dependent density matrixrenormalization group method (DMRG) , on tight-binding a r X i v : . [ c ond - m a t . s t r- e l ] F e b chains, to extract the evolution on longer time scales.In this paper, we re-examine the TDNRG method, general-izing it to arbitrary finite temperature, for a general initial stateensemble, within the full density matrix (FDM) approach. We prove a number of exact results, including the exactnessof the short-time limit of local observables following a quan-tum quench, and a trace conserving property of the projectedfull density matrix. These prove useful in testing the finitetemperature formalism. By applying this approach to the pro-totype model of strongly correlated electron systems, the An-derson impurity model, we analyze in detail, and shed fur-ther light on, the origin of noise in the TDNRG method andgive a detailed description of errors in the local occupation anddouble occupancy in the long-time limit, identifying trends inswitching protocols, bath discretization and quench size thatminimize these errors. The results suggest that the short, in-termediate and particularly the long-time behavior can be im-proved by replacing a large quantum quench by a sequence ofsmaller quantum quenches, while retaining the standard Wil-son chain for the bath. With this motivation, the formalism forsuch multiple quenches is derived and shown to rest solely onthe NRG approximation [Eq. (10)], as in our finite tempera-ture generalization of the single quench case. This formalismis important since, (a), it points a way to an overall improvedtime-dependence within a generalized multiple-quench TD-NRG approach, (b), it allows a general pulse (appropriatelydiscretized, if continuous) to be treated, and, (c), it includesperiodic switching as a special case. The last, (c), for squarepulses, has also been considered in Ref. 51, where, however,additional approximations, beyond the NRG approximation,were made, and without use of the FDM (for a discussion ofthis, see the end of Sec. VI).Besides the NRG, a large number of other methods are be-ing used to investigate real-time and non-equilibrium dynam-ics of correlated systems. These include analytic approaches,such as the functional renormalization group, real-timerenormalization group, perturbative scaling approach, Keldysh perturbation theory, real-time and renormalizedperturbation theory, flow-equation, , dual-fermion ,slave-boson, and 1 / N -expansion techniques . Applica-tions of these to a number quantum impurity models have beenmade, including, to the anisotropic Kondo and spin-bosonmodels, to the interacting resonant level model , andto the Anderson impurity model . One advantage ofthese analytic approaches over the TDNRG, and other numer-ical approaches, is their ability to deal with a continuous bath.On the other hand, most of these approaches are perturbative,and they may have di ffi culty accessing the low temperaturestrong coupling limit. Numerical methods include the afore-mentioned DMRG method , also at finite temperature, and several quantum Monte Carlo approaches. The for-mer is limited, as compared to TDNRG, to finite size systemsand is not applicable to exponentially long times, whereasthe latter become computationally expensive for strong in-teractions or low temperatures. Indeed, the main motivationfor developing the TDNRG approach, is, on the one hand,its inherently non-perturbative nature, allowing strongly cor-related systems to be treated accurately at moderate compu- tational cost, with, on the other hand, its ability to describearbitrary low energy scales and temperatures. The latter, inprinciple, has the prospect of addressing accurately the time-dependence of such systems at arbitrarily long times. Theseadvantages, however, come at the expense of using a logarith-mically discretized bath, which incurs errors that we shall ad-dress in this paper. Finally, we mention that most of the abovetechniques, are, in principle, also of use in the description ofnon-equilibrium DMFT for correlated lattice models. The present work is also motivated by the increasing num-ber of experiments probing time-dependent properties of cor-related systems, including time-dependent spectroscopies ofcorrelated electron systems, such as pump-probe investiga-tions of Mott insulators, quasiparticle lifetime e ff ects insurface states, coherent control and relaxation times in solid-state qubits, determination of relaxation rates of excitedspin states of atoms on surfaces via voltage pulses, and non-equilibrium e ff ects in cold atom systems. The paper is organized as follows. In Sec. II, we providerequired background information (Sec. II A-II C), discuss lim-iting cases and exact results (Sec. II D), and state the prob-lem that has to be overcome in generalizing the existing TD-NRG approach of Refs. 14 and 40 to finite temperature withinthe FDM formalism (Sec. II E). This generalization requirescalculating all terms appearing in a projected density matrix,which is accomplished in Sec. III. In Sec. IV we derive re-cursion relations for these terms, which allow them to be cal-culated in a numerically e ffi cient manner. We then apply thisfinite-temperature generalization of the TDNRG approach tothe Anderson impurity model in Sec. V, providing detailednumerical results for local observables on all time and tem-perature scales. In particular we discuss the accuracy of thisapproach for the occupation number and double occupancy atdi ff erent temperature and time scales. In Sec. VI, we gener-alize the TDNRG approach from a single instantaneous quan-tum quench to a sequence of multiple quantum quenches overa finite time interval. The primary motivation for this gener-alization was stated above, namely to obtain an overall im-provement in the accuracy of the TDNRG method at longertimes. However, in the process we generalize the approachto other cases of great interest in their own right, such as toperiodic switching / driving and to general pulses (suitably dis-cretized, if continuous). A summary and outlook on futureprospects and applications of our formalism is presented inSec. VII. The trace conserving property of the projected fulldensity matrix is shown in Appendix A, and Appendix B givesa detailed proof that the short-time limit of observables is ex-act within TDNRG. Appendix C gives explicit expressions forrecursion relations for the projected density matrix for the caseof SU(2) spin symmetry. We used this symmetry in obtain-ing the numerical results for the Anderson impurity model inSec. V. Additional numerical results for the Anderson modelare given in Appendix D, while further details of the TDNRGderivation for multiple quenches are given in Appendix E. II. PRELIMINARIES
The NRG and its extension to time-dependent problemsapply generally to any quantum impurity model. For thepurposes of this section, we outline these, together with thegeneric quantum quench of interest in Sec. II A, postponing adetailed description of the specific model that we shall uselater in the numerical calculations (the Anderson impuritymodel) and the specific quench protocols to Sec. V. Somegeneral aspects of the NRG together with the complete basisformalism , and the FDM are introduced in Sec. II B. Werecapitulate the use of the complete basis set to derive a formalexpression for the time-dependence of a local observable fol-lowing a quantum quench in Sec. II C. Limiting cases and ex-act results that hold generally within the time-dependent NRGare given in Sec. II D. Finally, Sec. II E states the problem tobe overcome in generalizing the approach of Refs. 14 and 40to finite temperature, using the FDM, the main aim of thepresent paper. A. Quantum impurity models and quantum quenches
The Hamiltonian of a quantum impurity system is given by H = H imp + H int + H bath , (1)where H imp represents the impurity, a small quantum mechan-ical system with a few degrees of freedom and local many-body interactions, H bath represents a quasi-continuous bathof non-interacting particles, and H int is the interaction be-tween the two. Typical examples include the Anderson im-purity model (Sec. V A), the anisotropic Kondo model, orthe spin-boson model of two-level systems coupled to anenvironment. We are interested in the dynamics of a lo-cal observable ˆ O following a quantum quench in which one ormore system parameters of H change suddenly at t =
0. Thus,the time dependence of H is described by H ( t ) = θ ( − t ) H i + θ ( t ) H f , with H i and H f being time-independent initial ( t < t >
0) Hamiltonians, respectively. The timeevolution of ˆ O at t > O ( t ) = Tr (cid:104) ρ ( t ) ˆ O (cid:105) where ρ ( t ) = e − iH f t ρ e iH f t is the time-evolved density matrixand ρ = e − β H i / Tr[ ρ ] is the density matrix of the initial state atinverse temperature β . B. NRG, complete basis set, and full density matrix
Since the initial and final state Hamiltonians H i and H f are time-independent they can be iteratively diagonalized inthe usual way within the NRG method. In brief, this ap-proach, described in detail in Refs. 1–4, yields the eigen-states and eigenvalues of a sequence of truncated Hamilto-nians H i / fm , m = , , . . . , which approximate the spectra of H i / f , on successively decreasing energy scales ω m ∼ Λ − m / .The discretization parameter Λ > H i / f , such that an iterative di-agonalization scheme can be successfully carried out. This procedure is performed up to a maximum iteration m = N (“the longest Wilson chain”). At each m , the states generated,denoted | qm (cid:105) , are partitioned into the lowest energy retainedstates, denoted | km (cid:105) , and the high energy eliminated (or dis-carded) states, | lm (cid:105) . In order to avoid an exponential increasein the dimension of the Hilbert space, only the former are usedto set up and diagonalize the Hamiltonian for the next itera-tion m +
1. The eliminated states, while not used in the it-erative NRG procedure, are nevertheless crucial as they areused to set up a complete basis set with which the expres-sions for the time dependent dynamics are evaluated. Thiscomplete basis set is defined by the product states | lem (cid:105) = | lm (cid:105)| e (cid:105) , m = m , . . . , N , where m is the first iteration at whichtruncation occurs, and | e (cid:105) = | α m + (cid:105)| α m + (cid:105) . . . | α N (cid:105) are environ-ment states at iteration m such that the product states | lem (cid:105) ,for each m = m , m + , . . . , N , reside in the same Fock space(that of the largest system diagonalized, m = N ). The α m rep-resent the configurations of site m in a linear chain representa-tion of the quantum impurity system (e.g. the four states 0, ↑ , ↓ and ↑↓ at site m for a single channel Anderson model) and“ e ” in | lem (cid:105) denotes the collection e = { α m + . . . α N } . Thesestates satisfy the completeness relation N (cid:88) m = m (cid:88) le | lem (cid:105)(cid:104) lem | = , (2)where for m = N all states are counted as discarded (i.e. thereare no kept states at iteration m = N ). We shall also use thefollowing representations of this relation, = − m + + m , (3)1 − m = m (cid:88) m (cid:48) = m (cid:88) l (cid:48) e (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105)(cid:104) l (cid:48) e (cid:48) m (cid:48) | (4)1 + m = N (cid:88) m (cid:48) = m + (cid:88) l (cid:48) e (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105)(cid:104) l (cid:48) e (cid:48) m (cid:48) | = (cid:88) ke | kem (cid:105)(cid:104) kem | . (5)By using the complete basis set, we can construct an ini-tial state density matrix which is valid at any temperature, theFDM, ρ = N (cid:88) m = m w m ˜ ρ m , (6)˜ ρ m = (cid:88) le | lem (cid:105) i e − β E ml ˜ Z m i (cid:104) lem | , (7)which takes into account the discarded states of H i from allshells. For later use, we note that, (a), Tr (cid:2) ˜ ρ m (cid:3) = Tr (cid:2) ρ (cid:3) = (cid:80) Nm = m w m =
1, and, (b), Tr (cid:2) ˜ ρ m (cid:3) = = (cid:80) le e − β Eml ˜ Z m = (cid:80) l d N − m e − β Eml ˜ Z m = d N − m Z m ˜ Z m where Z m = (cid:80) l e − β E ml , i.e., ˜ Z m = d N − m Z m , and d is the degeneracy of theWilson site α m . C. Time-dependence
With the above definitions, we can write for the time evolu-tion of a local observable ˆ OO ( t ) = Tr (cid:104) e − iH f t ρ e iH f t ˆ O (cid:105) = N (cid:88) m = m (cid:88) le f (cid:104) lem | e − iH f t ρ e iH f t ˆ O | lem (cid:105) f = N (cid:88) mm (cid:48) = m (cid:88) lel (cid:48) e (cid:48) f (cid:104) lem | e − iH f t ρ e iH f t | l (cid:48) e (cid:48) m (cid:48) (cid:105) f f (cid:104) l (cid:48) e (cid:48) m (cid:48) | ˆ O | lem (cid:105) f . (8)Here, and in the remainder of the paper, A i / f denotesan expression in the basis of initial / final states. Using (cid:80) ke | kem (cid:105)(cid:104) kem | = (cid:80) Nm (cid:48) = m + (cid:80) l (cid:48) e (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105)(cid:104) l (cid:48) e (cid:48) m (cid:48) | [i.e., the sec-ond equality in Eq. (5)], we can express Eq. (8) as O ( t ) = N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) (cid:88) e f (cid:104) sem | e − iH f t ρ e iH f t | rem (cid:105) f f (cid:104) rem | ˆ O | sem (cid:105) f ≈ N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) (cid:88) e f (cid:104) sem | e − iH fm t ρ e iH fm t | rem (cid:105) f f (cid:104) rem | ˆ O | sem (cid:105) f = N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) (cid:16) (cid:88) e f (cid:104) sem | ρ | rem (cid:105) f (cid:17) e − i ( E ms − E mr ) t O mrs = N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) ρ i → fsr ( m ) e − i ( E ms − E mr ) t O mrs , (9)in which r and s may not both be kept states, O mrs = f (cid:104) lem | ˆ O | rem (cid:105) f are the final state matrix elements of ˆ O , whichare independent of e , the NRG approximation H f | rem (cid:105) ≈ H fm | rem (cid:105) = E mr | rem (cid:105) , (10)is adopted [in the second line of Eq. (9)], and ρ i → fsr ( m ) = (cid:80) e f (cid:104) sem | ρ | rem (cid:105) f represents the reduced density matrix of theinitial state projected onto the basis of final states (henceforthsimply called the projected density matrix ). D. Limiting cases and exact results
A number of limiting cases and exact results follow gener-ally from the above formalism. First, it is clear from Eq. (9)that in the absence of a quench, H i = H f , one recovers for O ( t )the time independent equilibrium thermodynamic average O ( t ) = O i = O f . (11)Second, in the special case that ˆ O is the identity operator, wehave from Eq. (9) upon using O mrs = δ rs that1 = N (cid:88) m = m (cid:88) l ρ i → fll ( m ) , (12) i.e., the trace of the projected density matrix over the elim-inated states is conserved (see Appendix A). The above tworesults will serve as useful checks on the correctness of ourfinite temperature generalization of the time dependent NRGas well as on the precision of the numerical calculations inSec. V.Third, it has been noted previously that the TDNRG yieldsvery accurate results for the short-time limit of observables,i.e., to a very good approximation O ( t → + ) ≈ O i where O i = Tr (cid:104) ρ ˆ O (cid:105) is the thermodynamic value in the initial state. In fact, we can show that the short-time limit is exact, i.e., O ( t → + ) = O i . An explicit proof of this, starting from theexpression for O ( t → + ), is given in Appendix B. One canalso see this from the following argument: the time evolutionin Eq. (9) involves approximate NRG eigenvalues via the fac-tor e i ( E mr − E ms ) t . The NRG approximation for the eigenvalues istherefore not operative for t → + and Eq. (9) is as exact asEq. (8) in this limit (the latter yielding O i = Tr[ ρ ˆ O ]), hencewe have the exact result, O ( t → + ) ≡ N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) ρ i → fsr ( m ) O mrs = O i . (13)This will also be verified numerically in Sec. V.In contrast, as soon as t is finite, we expect that the time evo-lution in Eq. (9) will di ff er from that in Eq. (8), and this will,in part, be due to the appearance of approximate NRG eigen-values in the former, resulting in errors and noise in the timeevolution, which we shall analyze numerically in Sec. V C.We thus also expect, and find, that the long-time limit of O ( t ),given by O ( t → ∞ ) = N (cid:88) m = m (cid:88) rs ∈ DD (cid:48) ρ i → fsr ( m ) δ E ms , E mr O mrs , (14) = N (cid:88) m = m (cid:88) l ρ i → fll ( m ) O mll , (15)is not exact. The deviation of this from its expected thermody-namic value in the final state, O f , will be extensively analyzedin Sec. V E. E. Generalizing the time-dependent NRG to finite temperature
We now return to the formal expression for O ( t ) in Eq. (9)and consider its explicit evaluation. The matrix elements O mrs in Eq. (9) may be calculated recursively at each shell, in theusual way within the NRG. Calculating the projected den-sity matrix ρ i → fsr ( m ) requires more e ff ort and is the main prob-lem in generalizing the TDNRG to finite temperatures. In or-der to see this more clearly, we use a modification of Eqs. (3)-(5) for the decomposition of unity as in Ref. 401 = I + m + I − m = (cid:88) qe | qem (cid:105) ii (cid:104) qem | + m − (cid:88) m (cid:48) = m (cid:88) l (cid:48) e (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105) ii (cid:104) l (cid:48) e (cid:48) m (cid:48) | , (16)in which the label q includes both retained and eliminatedstates at iteration m , and break ρ i → fsr ( m ) into four terms ρ i → fsr ( m ) = (cid:88) e f (cid:104) sem | ( I + m + I − m ) ρ ( I + m + I − m ) | rem (cid:105) f = ρ ++ sr ( m ) + ρ + − sr ( m ) + ρ − + sr ( m ) + ρ −− sr ( m ) , (17)with ρ pp (cid:48) sr ( m ) = (cid:80) e f (cid:104) sem |I pm ρ I p (cid:48) m | rem (cid:105) f and p , p (cid:48) = ± . Calcu-lating ρ ++ sr ( m ) is straightforward, since it includes only overlapmatrix elements between initial and final states at the sameshell, i (cid:104) qem | rem (cid:105) f . But for ρ + − sr ( m ), ρ − + sr ( m ), and ρ −− sr ( m ) this ismanifestly not the case, requiring overlap matrix elements be-tween initial and final states at di ff erent shells, i (cid:104) l (cid:48) e (cid:48) m (cid:48) | rem (cid:105) f with m (cid:48) < m .In previous work, the problem of evaluating ρ + − sr ( m ), ρ − + sr ( m ), and ρ −− sr ( m ), for general situations, was avoided bychoosing a special form for the initial state density matrix,namely ρ = (cid:88) l | lN (cid:105) i e − β E Nl Z N i (cid:104) lN | , (18)with Z N = (cid:80) l e − β E Nl , in which only the discarded states ofthe last NRG iteration enter, while discarded states at m < N are neglected. This simplifies the projected density matrix to ρ i → fsr ( m ) = ρ ++ sr ( m ) as ρ + − sr ( m ) = ρ − + sr ( m ) = ρ −− sr ( m ) = m ≤ N for this special choice of ρ . However, this limits thecalculations to T ∼ T N where T N ∼ ω N is the characteris-tic scale of H iN . Calculations at higher temperature T > T N ,within this approach, would require choosing a shorter chainof length M < N such that T ≈ T M , which introduces addi-tional errors due to the shorter chain. Hence, a formulationthat uses the FDM of the initial state [Eq. (6)], which encodesall discarded states from all shells m ≤ N , would be advan-tageous, since it would be valid at all temperatures T ≥ T N automatically, while a fixed chain of length N is used for all T . In the next section, we show that this is possible, hence al-lowing the calculation of the time-dependence at an arbitraryfinite temperature starting from a general initial state. Weshall thus use the FDM in Eq. (6), which takes into accountthe discarded states from all shells, and present the calcula-tion of each of the four terms appearing in ρ i → fsr ( m ), especially ρ + − sr ( m ), ρ − + sr ( m ), and ρ −− sr ( m ). III. FINITE-TEMPERATURE PROJECTED DENSITYMATRIX
In this section we evaluate the four terms contributing to ρ i → fsr ( m ) in Eq. (17) starting from the FDM of the initial state ρ = (cid:80) Nm = m w m ˜ ρ m . Substituting the latter into the first term of ρ i → fsr ( m ) in Eq. (17), we have ρ ++ sr ( m ) = N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) e f (cid:104) sem |I + m ˜ ρ m (cid:48) I + m | rem (cid:105) f = N (cid:88) m (cid:48) = m (cid:88) eqe (cid:48) q (cid:48) e (cid:48)(cid:48) f (cid:104) sem | qe (cid:48) m (cid:105) ii (cid:104) qe (cid:48) m | w m (cid:48) ˜ ρ m (cid:48) | q (cid:48) e (cid:48)(cid:48) m (cid:105) ii (cid:104) q (cid:48) e (cid:48)(cid:48) m | rem (cid:105) f , in which the overlap matrix elements, i (cid:104) qe (cid:48) m | rem (cid:105) f = S mq i r f δ ee (cid:48) , between initial and final states, are diagonal in, andindependent of the environment degrees of freedom. Then ρ ++ sr ( m ) = N (cid:88) m (cid:48) = m (cid:88) qq (cid:48) S ms f q i (cid:88) e i (cid:104) qem | w m (cid:48) ˜ ρ m (cid:48) | q (cid:48) em (cid:105) i S mq i r f = N (cid:88) m (cid:48) = m (cid:88) qq (cid:48) S ms f q i (cid:88) e , le (cid:48) i (cid:104) qem | le (cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) i (cid:104) le (cid:48) m (cid:48) | q (cid:48) em (cid:105) i S mq (cid:48) i r f . Since (cid:80) q | qem (cid:105) = (cid:80) l | lem (cid:105) + (cid:80) k | kem (cid:105) , i (cid:104) l (cid:48) e (cid:48) m (cid:48) | lem (cid:105) i = δ ll (cid:48) δ ee (cid:48) δ mm (cid:48) , and i (cid:104) le (cid:48) m (cid:48) | kem (cid:105) i = m (cid:48) ≤ m , the above equa-tion is equivalent to ρ ++ sr ( m ) = (cid:88) le S ms f l i w m e − β E ml ˜ Z m S ml i r f + N (cid:88) m (cid:48) = m + (cid:88) kk (cid:48) S ms f k i (cid:88) ele (cid:48) i (cid:104) kem | le (cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) i (cid:104) le (cid:48) m (cid:48) | k (cid:48) em (cid:105) i S mk (cid:48) i r f . Using the definition of the full reduced density matrix R m red ( k , k (cid:48) ) = N (cid:88) m (cid:48) = m + (cid:88) ele (cid:48) i (cid:104) kem | le (cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) i (cid:104) le (cid:48) m (cid:48) | k (cid:48) em (cid:105) i , (19)and ˜ Z m = d N − m Z m from Sec. II, we have ρ ++ sr ( m ) = (cid:88) l S ms f l i w m e − β E ml Z m S ml i r f + (cid:88) kk (cid:48) S ms f k i R m red ( k , k (cid:48) ) S mk (cid:48) i r f . (20)In the other terms with I − m , we notice that I − m ˜ ρ m (cid:48) = ˜ ρ m (cid:48) I − m = (cid:40) ˜ ρ m (cid:48) if m (cid:48) < m ;0 otherwise , then ρ + − sr ( m ) = N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) e f (cid:104) sem |I + m ˜ ρ m (cid:48) I − m | rem (cid:105) f = m − (cid:88) m (cid:48) = m w m (cid:48) (cid:88) eqe (cid:48) f (cid:104) sem | qe (cid:48) m (cid:105) ii (cid:104) qe (cid:48) m | ˜ ρ m (cid:48) | rem (cid:105) f = m − (cid:88) m (cid:48) = m (cid:88) eqe (cid:48) le (cid:48)(cid:48) f (cid:104) sem | qe (cid:48) m (cid:105) ii (cid:104) qe (cid:48) m | le (cid:48)(cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) i (cid:104) le (cid:48)(cid:48) m (cid:48) | rem (cid:105) f . This term equals zero since i (cid:104) qe (cid:48) m | le (cid:48)(cid:48) m (cid:48) (cid:105) i = m (cid:48) < m (i.e., discarded states | le (cid:48)(cid:48) m (cid:48) (cid:105) at iteration m (cid:48) < m have no over-lap with states | qe (cid:48) m (cid:105) of later iterations). Similarly, ρ − + sr ( m ) =
0. Finally, the last term contributing to ρ i → fsr ( m ) is given by ρ −− sr ( m ) = N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) e f (cid:104) sem |I − m ˜ ρ m (cid:48) I − m | rem (cid:105) f = m − (cid:88) m (cid:48) = m (cid:88) e , le (cid:48) f (cid:104) sem | le (cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) i (cid:104) le (cid:48) m (cid:48) | rem (cid:105) f . Inserting 1 = + m (cid:48) + − m (cid:48) from Eq. (3) into the overlap matrixelements appearing above, we have that f (cid:104) sem | le (cid:48) m (cid:48) (cid:105) i = f (cid:104) sem | (1 + m (cid:48) + − m (cid:48) ) | le (cid:48) m (cid:48) (cid:105) i = (cid:88) k f (cid:104) sem | ke (cid:48) m (cid:48) (cid:105) f f (cid:104) ke (cid:48) m (cid:48) | le (cid:48) m (cid:48) (cid:105) i + m (cid:48) (cid:88) m (cid:48)(cid:48) = m (cid:88) l (cid:48) e (cid:48)(cid:48) f (cid:104) sem | l (cid:48) e (cid:48)(cid:48) m (cid:48)(cid:48) (cid:105) f f (cid:104) l (cid:48) e (cid:48)(cid:48) m (cid:48)(cid:48) | le (cid:48) m (cid:48) (cid:105) i = (cid:88) k f (cid:104) sem | ke (cid:48) m (cid:48) (cid:105) f S m (cid:48) k f l i (21)since f (cid:104) sem | l (cid:48) e (cid:48)(cid:48) m (cid:48)(cid:48) (cid:105) f = m (cid:48)(cid:48) < m , then ρ −− sr ( m ) = m − (cid:88) m (cid:48) = m (cid:88) ee (cid:48) kk (cid:48) f (cid:104) sem | ke (cid:48) m (cid:48) (cid:105) f × (cid:34) (cid:88) l S m (cid:48) k fli (cid:16) w m (cid:48) e − β E m (cid:48) l ˜ Z m (cid:48) (cid:17) S m (cid:48) l i k (cid:48) f (cid:35) f (cid:104) k (cid:48) e (cid:48) m (cid:48) | rem (cid:105) f . (22)In general, we proved that ρ + − sr ( m ) = ρ − + sr ( m ) = ρ i → fsr ( m ) = ρ ++ sr ( m ) + ρ −− sr ( m ) , (23)without the restriction of using a density matrix defined forthe last iteration, therefore the calculation of O ( t ) is possibleat arbitrary temperatures for a general initial state. Further-more, in Eq. (22), ρ −− sr ( m ) is no longer expressed in terms ofoverlap matrix elements ( S -matrix elements) between initialand final states belonging to di ff erent shells but instead in-volves only shell diagonal matrix elements. This simplifiesthe calculation of ρ −− sr ( m ) making it as simple as the calcula-tion of R m red ( k i , k (cid:48) i ) : notice the structural similarities betweenEq. (19) and Eq. (22). The R m red are obtained e ffi ciently bytracing out shell degrees of freedom backwards starting from N in a recursive manner. Similarly, the ρ −− ( m ) can be ob-tained recursively by a forward iteration procedure startingfrom m = m . This forward recursive procedure, described inthe next section, makes the calculation of ρ −− ( m ) as e ffi cientas that for R m red . IV. RECURSION RELATIONS FOR THE PROJECTEDDENSITY MATRIX
As shown in the last section, the projected density matrixis simplified into two contributions [Eq. (23)]. In this section,we show how to calculate them recursively. For convenience,we rewrite Eq. (23) as ρ i → fsr ( m ) = ˜ ρ ++ sr ( m ) + ρ sr ( m ) + ρ −− sr ( m ) , (24)where ˜ ρ ++ sr ( m ) = (cid:88) kk (cid:48) S ms f k i R m red ( k , k (cid:48) ) S mk (cid:48) i r f , (25) ρ sr ( m ) = (cid:88) l S ms f l i (cid:16) w m e − β E ml Z m (cid:17) S ml i r f . (26) ρ sr ( m ) can be calculated easily at each shell once we haveall the eigenvalues and overlap matrix elements. ˜ ρ ++ sr ( m ) canalso be calculated at each shell, once R m red ( k i , k (cid:48) i ) has been cal-culated recursively. Here, we present the recursive calcu-lation for ρ −− sr ( m ), with the calculation for R m red ( k i , k (cid:48) i ) beingsimilar, Using the above definition of ρ sr ( m ) and Z m (cid:48) = Z m (cid:48) d N − m (cid:48) ,we have that Eq. (22) is equivalent to ρ −− sr ( m ) = m − (cid:88) m (cid:48) = m d N − m (cid:48) (cid:88) ee (cid:48) kk (cid:48) f (cid:104) sem | ke (cid:48) m (cid:48) (cid:105) f ρ kk (cid:48) ( m (cid:48) ) f (cid:104) k (cid:48) e (cid:48) m (cid:48) | rem (cid:105) f . We note that ρ −− sr ( m ) =
0, and ρ −− sr ( m + = d N − m (cid:88) ee (cid:48) kk (cid:48) f (cid:104) se ( m + | ke (cid:48) m (cid:105) f ρ kk (cid:48) ( m ) f (cid:104) k (cid:48) e (cid:48) m | re ( m + (cid:105) f = d (cid:88) α m + kk (cid:48) A α m + † sk ρ kk (cid:48) ( m ) A α m + k (cid:48) r (27)since f (cid:104) k (cid:48) e (cid:48) m | re ( m + (cid:105) f = δ e (cid:48) e A α m + k (cid:48) r with e (cid:48) = { e , α m + } ,and (cid:80) e = d N − ( m + . In the above we used the notation ofRef. 16 for the transformation matrix A α m + k (cid:48) r relating eigen-states | r ( m + (cid:105) of H m + to the product basis | k (cid:48) m (cid:105)| α m + (cid:105) ,i.e., | r ( m + (cid:105) = (cid:88) k (cid:48) α m + A α m + k (cid:48) r | k (cid:48) m (cid:105)| α m + (cid:105) . (28)Similarly, we have that ρ −− sr ( m + = d (cid:88) α m + kk (cid:48) A α m + † sk ρ kk (cid:48) ( m + A α m + k (cid:48) r + d (cid:88) α m + α m + kk (cid:48) k k (cid:48) A α m + † sk A α m + † kk ρ k k (cid:48) ( m ) A α m + k (cid:48) k (cid:48) A α m + k (cid:48) r . Using the definition of ρ −− sr ( m +
1) in Eq. (27), we obtain ρ −− sr ( m + = d (cid:88) α m + kk (cid:48) A α m + † sk [ ρ kk (cid:48) ( m + + ρ −− kk (cid:48) ( m + A α m + k (cid:48) r . Consequently, we have the recursion relation ρ −− sr ( m ) = m = m ;1 d (cid:88) α m kk (cid:48) A α m † sk (cid:2) ρ kk (cid:48) ( m − + ρ −− kk (cid:48) ( m − (cid:3) A α m k (cid:48) r otherwise . (29)In this relation, ρ −− sr ( m ) is no longer expressed in terms of allthe earlier shells as in Eq. (22), but just in terms of one ear-lier shell, ( m − R m red ( k , k (cid:48) ) = m = N ; (cid:88) l α m + A α m + kl (cid:16) w m + e − β E m + l Z m + (cid:17) A α m + † lk (cid:48) + (cid:88) k k (cid:48) α m + A α m + kk R m + ( k , k (cid:48) ) A α m + † k (cid:48) k (cid:48) otherwise , (30)in which R m red ( k , k (cid:48) ) depends only on terms of one later shell,( m + ρ −− sr ( m ) re-cursively in a single forward run, while Eq. (30) may be usedto determine R m red ( k , k (cid:48) ) [and hence ρ ++ sr ( m ) from Eq. (20)] re-cursively in a backward run. While the ρ ++ sr ( m ) is the contri-bution to the projected density matrix at iteration m arisingfrom lower energy states at subsequent iterations m (cid:48) ≥ m , ρ −− sr ( m ) is the contribution to the projected density matrix atiteration m arising from higher energy states at earlier itera-tions m (cid:48) < m [see Eq. (22)]. The former is finite in the equi-librium limit H f = H i , yielding the correct Boltzmann factorsin thermodynamic averages, while the latter is finite only fornon-equilibrium H f (cid:44) H i [as can be seen from Eq. (22)]. Fi-nally, the modifications to the recursive expressions Eq. (29)and Eq. (30) in the case where SU(2) symmetry is used aregiven in Appendix C. V. NUMERICAL RESULTS
In this section, we apply the finite-temperature TDNRG ap-proach to the Anderson impurity model, a prototype modelof strong electronic correlations. Section V A describes themodel and the switching protocols that we shall consider. Sec-tion V B tests the trace conservation property [Eq. (12)] of theprojected density matrix ρ i → f . This is both a useful test ofthe correctness of our expression for ρ i → f as calculated fromthe recursion relations in Eqs. (24)-(30), and also indicates therelative importance of the various contributions to this quan-tity at di ff erent temperatures. The time-dependence of observ-ables is presented in Sec. V C, where we attempt to identifywhich contributions result in noise in the long-time dynamicsand how this noise is a ff ected by averaging over discretiza-tions of the band. Sections V D-V E present results for theshort ( t → + ) and long-time ( t → + ∞ ) limits of local observ-ables, such as the impurity occupation and double occupation,discussing in particular how di ff erent switching protocols, theorder of switchings, and quench strengths a ff ect the accuracyof the long-time limit. A comparison between results withinthe present TDNRG approach and that of Anders and Schillerin Refs. 14 and 40 is also presented [Fig. 3 and Fig. 8]. A. Anderson Impurity Model
We illustrate the finite-temperature TDNRG method by ap-plying it to the Anderson impurity model. This is defined by Eq. (1) with H imp = (cid:88) σ ε d ( t ) n d σ + U ( t ) n d ↑ n d ↓ , (31) H bath = (cid:88) k σ (cid:15) k c † k σ c k σ , (32) H int = (cid:88) k σ V ( t )( c † k σ d σ + d † σ c k σ ) , (33)where n d σ is the number operator for electrons with spin σ in a local level with energy ε d ( t ), U ( t ) is the Coulomb re-pulsion between electrons in this level, V ( t ) is the hybridiza-tion matrix element of the local d state with the conductionstates, and (cid:15) k is the kinetic energy of the conduction elec-trons with wavenumber k . For simplicity, we shall consideronly quenches in which the d -level energy or Coulomb re-pulsion (or both) are changed at t = V ( t ) = V is kept constant. Specifically, we consider, (a), switchingthe d -orbital energy, ε d = θ ( − t ) ε i + θ ( t ) ε f while keeping U ( t ) = U i = U f constant, and, (b), switching the Coulombinteraction U ( t ) = θ ( − t ) U i + θ ( t ) U f with simultaneous changeof ε d ( t ) = θ ( − t ) ε i + θ ( t ) ε f . In the former, the interest is inswitching between di ff erent regimes of the Anderson model,e.g., from the mixed-valence regime with ε i = ε f = − U f /
2, and vice versa (or between twoKondo regimes with di ff ering level energies), while in the lat-ter it is of interest to investigate switching between an uncor-related ( U i =
0) and a correlated system ( U f > ε d , U and the hybridization strength Γ = πρ V , where ρ = / W is the density of states per spin for a flat band ofwidth W = D = D = U (cid:29) π Γ and − ε d (cid:29) Γ it also develops a renormalized low energy scale, the Kondoscale T K = √ U Γ / e πε d ( ε d + U ) / Γ U . Since ε d will often be var-ied below, whenever T K appears in the numerical calcula-tions below, it will denote the symmetric Kondo scale with U = max( U i , U f ). We shall mainly focus on the expectationvalue of the local level occupation (cid:104) n d ( t ) (cid:105) and of the doubleoccupancy (cid:104) n d ↑ ( t ) n d ↓ ( t ) (cid:105) . In the results reported below, weshall express temperature and time dependences in terms of T / T K and t Γ , respectively. Short, intermediate and long timeswill correspond to t Γ (cid:28) t Γ ∼ t Γ (cid:29)
1, respectively.The short-time ( t → + ) and long-time ( t → + ∞ ) limits of thetdNRG will be tested against the exactly known results (ob-tained, for example, from thermodynamic calculations for theinitial / final state Hamiltonians). This allows the accuracy ofthe TDNRG approach to be benchmarked at finite tempera-ture. B. Trace of ρ i → fsr Since the projected density matrix is the most importantquantity in our formulation of the finite-temperature TDNRG,we consider Tr[ ρ i → f ] vs temperature first. Appendix A showsthat Tr[ ρ i → f ] = -4 -3 -2 -1 T r [ ρ ] T/T K ρ i → f ρ ~ ++ ρ ρ -- FIG. 1. (Color online) Tr[ ρ i → f ] and its contributions vs tempera-ture. The system is switched from ε i = ε f = − Γ (symmetric Kondo regime). The other parameters are U i = U f = U = Γ , and Γ = − D . T K ≈ . × − D is the Kondotemperature in the final state. The calculations are for Λ =
2, no z averaging, and keeping 660 states per NRG iteration. case of time evolution in Eq. (9) when ˆ O =
1. This resultis useful as a test of the numerical accuracy of the calcula-tions as well as a way of estimating the relative contributionsof ˜ ρ ++ , ρ and ρ −− to ρ i → f and hence their importance for thetime-dependence. Figure 1 shows that the trace is conservedto within an error less than 10 − at all temperatures. Thesame accuracy holds also for any choice of ε d , U , and Λ . Thetrace of each contribution to ρ i → f , i.e., ˜ ρ ++ , ρ , and ρ −− , isalso shown in Fig. (1). We see that Tr[ ˜ ρ ++ ] has its main con-tribution at low to intermediate temperatures, while the maincontribution of Tr[ ρ ] is at intermediate to high temperatures.The contribution of Tr[ ρ −− ] can be as large as 20% at interme-diate temperatures, but vanishes only at T (cid:28) T K and T (cid:29) T K .These results exhibit the numerical precision of calculation,and also show the correctness of generalizing the TDNRG atfinite temperature, which essentially depends on ρ −− . C. Time-dependence
Figure 2 (a) shows the time dependence of the occupa-tion number n d ( t ) after switching the system from the mixed-valence regime to the symmetric Kondo regime at three rep-resentative temperatures T = . T K , T = T K ≈ Γ and T = T K ≈ W . The results at lower temperatures, T (cid:28) T K ,are qualitatively similar to those in the lowest-temperaturecase ( T = . T K ) shown in Fig. 2, namely, starting fromthe thermodynamic value of the initial state at t = n d ( t )varies smoothly at short ( t Γ (cid:28)
1) and intermediate t Γ ∼ t Γ (cid:38) T = T K ∼ Γ , n d ( t ) evolves from ahigher initial value at t =
0, as expected from the behavior ofthe occupation number with temperature in the mixed-valenceregime, and exhibits less noise at long times than the low -2 -1 n d ( t ) t Γ (b) N z =32 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) N z =1T=10 -1 T K T=10 +2 T K T=10 +5 T K -0.1 0 0.1 -0.1 0 0.1 0.95 1 0.95 1 FIG. 2. (Color online) Time dependence of occupation number n d ( t )at di ff erent temperatures with switching as in Fig. 1 (i.e. from themixed-valence to the symmetric Kondo regime with U = U i = U f = Γ ). (a) Calculation for Λ =
2, no z averaging, and keeping 660states per NRG iteration. The upper and lower insets, respectively,represent the contributions of the kept states and discarded states tothe occupation number in the range 1 ≤ t Γ ≤ . (b) Calculationas in (a), but with z averaging, using N z =
32. Insets as in (a) showcontributions from kept (upper) and discarded (lower) states. temperature curve. Since the final state is particle-hole sym-metric, one expects that n d ( t → + ∞ ) = n d ( t ) deviates from the exact value by 2% − T = T K ≈ W , the occupa-tion number n d ( t ) is almost independent of time and and liesvery close to its expected thermodynamic value. The noise isalso seen to be negligible in this limit. The main panel of Fig-ure 2 (b) shows the e ff ect of z averaging with N z =
32 onthe above results. We see that z averaging reduces the noiseat long times, but does not improve on the expected long-timelimit for n d (i.e., the error for the low temperature results re-mains 2% − (cid:88) rs (cid:60) KK (cid:48) = (cid:88) rs ∈ KD + (cid:88) rs ∈ DK + (cid:88) rs ∈ DD , we decompose O ( t )in Eq. (9) as follows O ( t ) = O KD ( t ) + O DK ( t ) + O DD (cid:48) ( t ). O KD ( t ) + O DK ( t ) represents the contribution from the keptstates, while O DD (cid:48) ( t ) represents the contribution from just thediscarded states. These two contributions to the total time evo-lution without and with the z averaging are respectively shownin the insets to Fig. 2 (a) and (b) in the time range t Γ ≥
1. Theupper inset to Fig. 2 (a) shows the contribution of kept states, n KDd ( t ) + n DKd ( t ) = × (cid:80) Nm = m (cid:80) ∈ KDkl ρ i → fkl ( m ) cos[( E mk − E ml ) t ] n mlk .This contribution, evolving in time as cos[( E mk − E ml ) t ], exhibitssignificant noise at long times. The features in the noise fromthis contribution are reflected in the total. We do not observesuch large noise in the remaining contribution n DD (cid:48) d ( t ), whichevolves in time as cos[( E ml − E ml (cid:48) ) t ] (lower inset). The insetsin Fig. 2 (b) also show the contributions of kept and discardedstates to the total time evolution, calculated with Nz = z averaging also reduces the noise appearing in the keptstate contribution, but has a small e ff ect on the discarded statecontribution. On the other hand, n KDd ( t ) + n DKd ( t ) equilibratesto zero in the long-time limit, whereas n DD (cid:48) d ( t ) makes up allthe contribution in this limit (as expected, since the long-timelimit consists of selecting states with E mr = E ms , implying thatthey are both discarded states). This implies that the noiseis mainly coming from the eigenvalues of the kept states, ormore generally, the NRG approximation, as the low-energykept states at each iteration are the ones most a ff ected by theNRG truncation procedure. On the other hand, the observederror in the long-time limit is likely not due to the NRG ap-proximation, since this approximation generally has a smallere ff ect on the high energy (discarded) states which are also theonly ones contributing to this limit. n d ( t ) t Γ N z =32presentAS 0.9 0.95 1 N z =1 FIG. 3. (Color online) Time dependence of the occupation num-ber n d ( t ) calculated within the present TDNRG (full line) and theAnders-Schiller approach (AS). Switching as in Fig. 1 with the sameparameters (mixed-valence to symmetric Kondo state). The calcula-tion is at T = T K with Λ = z averaging with N z =
32, andkeeps 660 states per NRG iteration. The inset represents the resultswith no z averaging. It is of interest to compare the present TDNRG calculationswith those from previous work . As mentioned in Sec. II,previous TDNRG calculations at finite temperature T werecarried out by choosing a shorter chain of length M < N , suchthat T ≈ T M = ω M where ω M is the characteristic energy scale of the Wilson chain of length M . In this approach, the initialstate density matrix includes only the contribution from all thestates of shell M [i.e., it uses Eq. (18) with N = M ]. Figure 3shows n d ( t ) calculated within the present TDNRG and in theprevious approach of Anders and Schiller (AS) in Refs. 14and 40 at T = T K , with z averaging. The two calculationsgive very similar results at short times, while there is a cleardi ff erence at long times. The long-time limit of the occupationnumber in the present work is closer to the expected value thanin the AS approach. We also see that the time evolution in theAS approach exhibits more noise than in the present one, bothwith and without z averaging (see the inset to Figure 3). Thesedi ff erences show that at finite temperature T ≈ T M , the useof a short chain of length M in the approach of AS does notcapture the time evolution of local observables as accuratelyas within our approach based on using a FDM for the initialstate. On the other hand, at low temperatures ( T < T K ), thechain of length M corresponding to T is su ffi ciently long, thatthe two approaches give very similar results. D. Short-time limit -4 -3 -2 -1 n d ( t → + ) T/T K ε i = -5 Γε i = -4 Γε i = -3 Γε i = -2 Γε i = -1 Γε i = 0 Γ FIG. 4. (Color online) Local level occupation number n d ( t ) in theshort-time limit vs temperature for di ff erent initial ε i keeping the final ε f = − Γ and U i = U f = U = Γ fixed (i.e., switching from theasymmetric model to the symmetric Kondo regime). Γ = − D , and T K ≈ . × − D is the symmetric Kondo temperature of the finalstate. The symbols represent results with the TDNRG; the lines arethe thermodynamic values in the initial state. The calculations arefor Λ =
2, no z averaging, and keeping 660 states per NRG iteration. In the short-time limit, t → + , the NRG approximationis not operative in Eq. 9, and, as indicated in Sec II, O ( t =
0) exactly equals O i , the thermodynamic value in the initialstate, calculated within the FDM approach . In Fig. 4, weshow the short-time limit of the occupation number as wellas the initial state thermodynamic value in a wide range oftemperatures. The system is switched from the asymmetricto the symmetric Kondo model. We can see that n d ( t = − . We have the sameprecision with other sets of parameters, i.e., ε d , U , and Λ , orother local observables, e.g., double occupancy. This providesa check on our generalization to arbitrary temperatures, on theexpression for ρ i → f entering the calculation of O ( t → + ) andon the accuracy of the numerical calculations.We have also tested the short-time limit within the AS ap-proach for the same system and switching parameters as inFig. 4. In this case, the short-time limit for n d ( t ) correspondsexactly to the thermodynamic value in the initial state, cal-culated within the conventional approach to thermodynamicsvia the NRG . As shown in Ref. 110, the FDM and con-ventional approaches to thermodynamics give essentially thesame results, so we conclude that the short-time limit of theAS approach is also numerically exact. E. Long-time limit and thermalization
In the long-time limit, t → + ∞ , we already clarified in theprevious section that the NRG approximation is not the mainsource of error . However, there is another source of error,the logarithmic discretization of the band. In order to estimatethe accuracy of the long-time limit, we define the followingpercentage relative error δ O ( t → + ∞ ) = × O ( t → + ∞ ) − O f O f , (34)in which O f is the thermodynamic value in the final state cal-culated within the FDM approach to thermodynamics. Thecloseness of the long-time limit to O f indicates either the ex-tent to which the system thermalizes at long times, or the ex-tent of the errors arising from the discretization (the two mayalso be related ).Figure 5 shows the relative error in n d ( t → + ∞ ) uponswitching from the asymmetric to the symmetric model[Fig. 5 (a)], and vice versa [Fig. 5 (b)]. We can see two trends.The first one is that the larger the quench, measured by thedi ff erence ∆ ε d between ε i and ε f in this case, the larger theerror in the long-time limit in both Figs. 5 (a) and 5 (b). Theerror in the former lies below 4% in all cases and in the latterbelow 10%. The error vanishes for ∆ ε d → T ≈ T K ≈ Γ = U / T ≈ T K = . Γ . Signatures ofthis, as shoulders, are present also in Fig. 5 (a) for the largestquenches. The second trend is the strong dependence of theerror on the magnitude of the largest incoherent excitation inthe final state, denoted by ε maxinc and defined by ε maxinc = max( | ε f | , | ε f + U | , Γ ) . (35)We see that ε maxinc = Γ for all final states of Fig. 5 (a), butincreases from 7 Γ to 12 Γ for the corresponding final states ofFig. 5 (b), with the error increasing monotonically with in-creasing ε maxinc . -4 -3 -2 -1 δ n d ( t → + ∞ ) T/T K (b) ε i =-6 Γε f = 0 Γε f =-1 Γε f =-2 Γε f =-3 Γε f =-4 Γε f =-5 Γ -3-2-1 0 (a) ε f =-6 Γε i = 0 Γε i =-1 Γε i =-2 Γε i =-3 Γε i =-4 Γε i =-5 Γ FIG. 5. (Color online) The (percentage) relative error of the lo-cal level occupation number in the long-time limit vs temperatures. U i = U f = U = Γ and Γ = − D . (a) Switching is from theasymmetric model to the symmetric Kondo regime, as in Fig. 4. (b)Reverse switching, from the symmetric Kondo regime to the asym-metric model [the ε f are indicated and are the same as the ε i in (a)]. T K ≈ . × − D is the Kondo temperature of the symmetric modeland calculations are for Λ = .
6, no z averaging, and keeping 660states per NRG iteration. We consider a di ff erent switching protocol in order to checkthat the above two trends persist. We switch both the Coulombinteraction U ( t ) = θ ( − t ) U i + θ ( t ) U f , and the level energy ε d ( t )such that ε d ( t ) + U ( t ) / =
0, thereby maintaining particle-holesymmetry in both the initial and final states. By doing so, weare considering switching from a strongly correlated system toa less correlated one (or vice versa). Since ε d ( t ) + U ( t ) / = n d ( t ) = ∆ U = | U f − U i | , thebigger the error, with a small violation of this trend only forthe case of switching to an uncorrelated system with U f = -6-4-2 010 -4 -3 -2 -1 δ 〈 n d ↑ n d ↓ 〉 ( t → + ∞ ) T/T K (b) U i =12 Γ U f = 0 Γ U f = 2 Γ U f = 4 Γ U f = 6 Γ U f = 8 Γ U f =10 Γ f =12 Γ U i = 0 Γ U i = 2 Γ U i = 4 Γ U i = 6 Γ U i = 8 Γ U i =10 Γ FIG. 6. (Color online) The (percentage) relative error of the doubleoccupation number in the long-time limit vs temperature.
Γ = − D ,and switching is applied maintaining the particle-hole symmetry, i.e., ε d ( t ) + U ( t ) / =
0. (a) Switching from a less correlated system toa more correlated one ( U f > U i ). (b) Reverse switching, from amore correlated to a less correlated system ( U f < U i ). The U f areindicated and are the same as the U i in (a). T K ≈ . × − D is theKondo temperature of the most correlated system and calculationsare for Λ = .
6, no z averaging, and keeping 660 states per NRGiteration. most correlated system with U f = Γ . We also note that theextrema in the error for the double occupancy are qualitativelysimilar to those found for the occupation number in Figs. 5 (a)and 5(b) above. Comparing now the same size of quench, inboth Figs. 6 (a) and 5(b), we see that the error is larger inFig. 6 (a) with ε maxinc constant at 6 Γ , than in Fig. 6 (b) with ε maxinc smaller and ranging from Γ to 5 Γ . We note also thatthe errors in the long-time limit of the double occupancy arelarger than the corresponding ones in the occupation discussedabove, reaching values as large as 25% for highly correlatedfinal states. This is due to the fact that the double occupancycan become very small in a correlated system, hence while theabsolute errors in the long-time limit for di ff erent observablesare actually found to be comparable, the relative error for thedouble occupancy can become large for highly correlated finalstates.Combining the observations from Fig. 6 with those of Fig. 5, we conclude that the error in the long-time limit of alocal observable depends strongly on, (a), the size of the quan-tum quench, and, (b), the magnitude of the highest energy in-coherent excitation in the final state. A possible reason for(b), is that the logarithmic discretization scheme used in theTDNRG calculations does not capture accurately enough thehigh energy excitations (see the next paragraph for supportingresults). Further discussion of the long-time limit, supportingthe above conclusions, is presented in Appendix D, where weshow results for additional switching protocols. Σ m m s he ll n d ( m ) m shell Λ =1.6 Λ =2.0 Λ =4.0 Λ =10. FIG. 7. (Color online) Λ dependence of the shell accumulated lo-cal level occupation number in the long-time limit at T = − T K (symbols). Switching is from the mixed-valence ( ε i =
0) to the sym-metric Kondo regime ( ε f = − Γ ) with U = U i = U f = Γ . TheKondo scale in the final state is T K ≈ . × − D . The calcula-tions are done with no z averaging and keeping 660 states per NRGiteration. The solid lines indicate the e ff ect of increasing the Wilsonchains above by six additional sites, maintaining the same T . The to-tal n d ( t → ∞ ) = (cid:80) N (cid:48) = N + m = m n d ( m ) is unchanged, indicating convergedresults (see the text for further discussion). In order to support the point made above concerning thelogarithmic discretization, we consider observables in thelong-time limit, calculated with di ff erent discretization pa-rameters, Λ , and at T = − T K (essentially zero tempera-ture). The system considered is that in Fig. 2 (a), with switch-ing being from the mixed-valence to the symmetric Kondoregime. The expected value of the occupation number in thelong-time limit is 1. Fig. 7 shows the shell accumulated valuefor n d ( t → + ∞ ), i.e. the quantity (cid:80) m shell m = m n d ( m ) where n d ( m )is the contribution to the long-time limit from shell m . Weobserve that the final result for m shell = N , with contributionsfrom all shells included, n d ( t → + ∞ ), is closer to the ex-pected value of 1 when Λ is closer to 1, the continuum limit.A smaller Λ resolves the spectrum better at high energies andcould account for the observation (b) above. The trend withdecreasing Λ suggests that the exact result can be obtained inthe limit Λ → + .It is worth commenting on whether the results for the totaloccupation number n d ( t → ∞ ) = (cid:80) Nm = m n d ( m ) shown in Fig. 7(symbols), are indeed converged with respect to the chain2 -6-5-4-3-2-1 0 (a) ε i = 0 Γ , ε f =-6 Γ presentAS 0 2 4 6 8 1010 -4 -3 -2 -1 δ n d ( t → + ∞ ) T/T K (b) ε i =-6 Γ , ε f = 0 Γ FIG. 8. (Color online) The (percentage) relative error of the locallevel occupation number in the long-time limit calculated with thepresent TDNRG (full line) and the approach of Ref. 14 and 40 (AS) atall temperatures. (a) The system is switched from the mixed-valenceregime with ε i = ε f = − Γ keeping U i = U f = U = Γ fixed and Γ = − D . The calculation isfor Λ = . N z =
32, and keeps 660 states per NRG iteration. (b)As in (a), but with the reverse switching from the symmetric Kondo( ε i = − Γ ) to the mixed-valence ( ε f = Γ ) regime. Other parametersare as in (a) and T K ≈ . × − D denotes the symmetric Kondoscale. length N = N ( Λ ) used in the calculations. For each Λ wechoose N such that Λ − ( N − / ≈ − T K which is su ffi cientlylong to carry out calculations for all temperatures down to thesmallest scale accessible, i.e., down to T ≈ − T K , the tem-perature used in Fig. 7. Nevertheless, in order to show thatthe above results are indeed converged, and that the observedtrend that n d ( t → ∞ ) comes closer to the exact value withdecreasing Λ is correct, we carried out calculations also forlonger Wilson chains and showed that n d ( t → ∞ ) remainedunchanged. This is illustrated by the solid lines in Fig. 7,which show the shell accumulated occupation numbers for aWilson chain of length N (cid:48) = N + T = − T K fixed.The final value of n d ( t → ∞ ) at shell N (cid:48) is unchanged, but theapproach of the shell accumulated occupation to this valueis modified. The upturn in the contribution to an observablefrom shells close to m = N T where Λ − ( N T − / ≈ T is a fea- ture of multiple shell calculations, but does not imply thatthe cumulative sum over all shells is not converged.Finally, we compare the long-time limit of the present TD-NRG with previous work in Refs. 14 and 40. Figure 8 exhibitsthe relative error of the occupation number in the long-timelimit with, (a), switching from the mixed-valence to the sym-metric Kondo regime, and, (b), vice versa. This figure illus-trates further details of the comparison discussed in Sec. V Cthat at low temperatures T < T K the two schemes give essen-tially the same results, while, a significant di ff erence arises inthe temperature range T K ≤ T < T K , with the presentTDNRG calculations showing improved accuracy for bothFig. 8 (a) and (b). VI. MULTIPLE QUENCHES AND GENERAL PULSES
From the results in Sec. V we see that the TDNRG approachbecomes increasingly more accurate at long times with de-creasing size of the quantum quench, while in the short-timelimit t → + it is always exact [see Eq. 13 and Fig. 4]. Wecan also show that the TDNRG calculation remains more ac-curate at short to intermediate times with decreasing quenchsize. To show this, we require a comparison to exact resultsat finite time. These are available for the non-interacting limitof the Anderson model U i = U f = U =
0, i.e. for the reso-nant level model. In Fig. 9 we show the (percentage) relativeerror δ n d ( t ) in n d ( t ) from TDNRG calculations, using the an-alytic result for the resonant level model as reference. Theerror is shown for several quantum quenches correspondingto varying ε i keeping ε f = Γ fixed, i.e. the quench sizesare ∆ ε d = ε f − ε i = Γ , Γ , Γ , Γ . From Fig. 9 we see thatfor su ffi ciently small quantum quenches ∆ ε d ≤ Γ , the relativeerror in n d ( t ) remains below 1% for times exceeding t Γ ∼ t Γ ≈ . n d ( t ) for the largest and smallest quenches. The analytic re-sults in both quenches saturate to the same value, lying ontop of each other at long times (since the final state energy ε f is fixed for all quenches), whereas we see that the TDNRGresult in the case of the smaller quench fits very well the an-alytic one in the whole time interval, while, in the case of thelarger quench, the TDNRG result fits the analytic one only upto t Γ (cid:46)
2, and deviates from this at longer times.In summary, then, the TDNRG is a controlled method forshort times and su ffi ciently small quantum quenches. Thissuggests a way to improve the long-time limit within a TD-NRG approach, by replacing a large quantum quench by a se-quence of smaller quantum quenches over a short time scale˜ τ n . The time evolution from the initial state at t < t > ˜ τ n is achieved by applying the TDNRG tothe sequence, i = , . . . , n , of small quantum quenches in theintervals ˜ τ i ≤ t < ˜ τ i + with ˜ τ = τ n = (cid:80) ni = τ i (seeFig. 10). For su ffi ciently small intervals, this should then al-low the TDNRG to access longer times more accurately thanwould be possible with a single large quantum quench.In this section we show how the above idea of using mul-tiple quenches to improve the long-time limit may be accom-3 -4 0 4 8 12 16 20 0 0.5 1 1.5 2 2.5 3 δ n d ( t ) t Γ ε i =-2 Γε i =-1 Γε i = 0 Γε i = 1 Γ n d ( t ) t Γε i =-2 Γε i = 1 Γ FIG. 9. (Color online) Percentage relative error in the occupationnumber calculated with TDNRG applied to the resonant level model,using the exact analytic results for reference. Calculations are withvarying ε i and fixed ε f = Γ . D = Γ , T = − Γ , Λ =
2, and N z =
32. The inset shows the time evolution of the occupation numbercalculated within TDNRG (dots) and the analytic expressions (solidlines) in the cases that the quenches are the biggest and smallest inthe main figure. plished within the formalism presented in Sec. III-IV. As il-lustrated schematically in Fig. 10, such an approach wouldalso allow the treatment of general continuous pulses actingin a finite time interval ˜ τ n , since a sequence of small quantumquenches acts as a discrete approximation for a continuouspulse. Such pulses, acting over a finite time interval, also cor-respond more closely to the actual situation in experiments.As a special case of the general multiple quench that we shalltreat, we mention the case of periodic switching between twostates (i.e., a train of square pulses). For a system driven through a set of quenches, as in Fig. 10,the time evolved density matrix at a general time in the inter-val ˜ τ p + > t ≥ ˜ τ p is given by ρ ( t ) = e − iH Qp + ( t − ˜ τ p ) e − iH Qp τ p ... e − iH Q τ ρ e iH Q τ ... e iH Qp τ p e iH Qp + ( t − ˜ τ p ) , (36)in which H Q p , p = , . . . , n are intermediate quench Hamilto-nians, acting during time intervals of length τ p , p = , . . . , n ,that determine the time evolution at intermediate times and H Q = H i and H Q n + = H f are the initial and final state Hamil-tonians, respectively (see Fig. 10).The time evolution of a local observable at ˜ τ p + > t ≥ ˜ τ p isexpressed in terms of the complete basis set of H Q p + as O ( t ) = (cid:88) mle Q p + (cid:104) lem | ρ ( t ) ˆ O | lem (cid:105) Q p + = (cid:88) mle Q p + (cid:104) lem | e − iH Qp + ( t − ˜ τ p ) ... e − iH Q τ ρ e iH Q τ ... e iH Qp + ( t − ˜ τ p ) × ˆ O | lem (cid:105) Q p + , (37)in which we used (cid:80) m for (cid:80) Nm = m . Inserting 1 Q p + = (cid:80) mle | lem (cid:105) Q p + Q p + (cid:104) lem | , the identity belonging to H Q p + , in ......... - τ ˜ τ ˜ τ n t H i = H Q H Q H Q H Q n H f = H Q n + τ z}|{ τ z}|{ τ n z}|{ ......... FIG. 10. A system driven from an initial to a final state via a se-quence of quantum quenches at times ˜ τ = , ˜ τ , . . . , ˜ τ n with evo-lution according to H Q p in the time interval ˜ τ p > t ≥ ˜ τ p − . Sucha sequence of multiple quantum quenches could also be used to de-scribe periodic switching or to approximate any general continuouspulse (smooth solid line). front of the operator ˆ O in the above equation, and using (cid:80) ke | kem (cid:105)(cid:104) kem | = (cid:80) Nm (cid:48) = m + (cid:80) l (cid:48) e (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105)(cid:104) l (cid:48) e (cid:48) m (cid:48) | , we have O ( t ) = (cid:60) KK (cid:48) (cid:88) mrse Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ ρ e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + × e − i ( E mr − E ms )( t − ˜ τ p ) O msr = (cid:60) KK (cid:48) (cid:88) mrs ρ i → Q p + rs ( m , ˜ τ p ) e − i ( E mr − E ms )( t − ˜ τ p ) O msr , (38)with ρ i → Q p + rs ( m , ˜ τ p ) = (cid:88) e Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ ρ e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + . Evidently, this projected density matrix can recover itscounterpart for a single quench since ρ i → Q rs ( m , ˜ τ ) = (cid:80) e Q (cid:104) rem | ρ | sem (cid:105) Q = ρ i → Q rs ( m ). Substituting the FDM of theinitial state into the above projected density matrix, we have ρ i → Q p + rs ( m , ˜ τ p ) = (cid:88) m l e e Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | l e m (cid:105) i × w m e − β E m l ˜ Z m i (cid:104) l e m | e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + . (39)We decompose ρ i → Q p + rs ( m , ˜ τ p ) into three terms, ρ i → Q p + rs ( m , ˜ τ p ) = ˜ ρ ++ rs ( m , ˜ τ p ) + ρ rs ( m , ˜ τ p ) + ρ −− rs ( m , ˜ τ p ) , (40)corresponding to the m > m , m = m , and m < m contribu-tions in Eq. (39) [compare with Eq. (24) for the single quench4case]. Explicitly written out, these are given by˜ ρ ++ rs ( m , ˜ τ p ) = N (cid:88) m = m + (cid:88) l e e Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | l e m (cid:105) i × w m e − β E m l ˜ Z m i (cid:104) l e m | e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + = N (cid:88) m = m + (cid:88) l e e Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ (1 i + m + i − m ) | l e m (cid:105) i × w m e − β E m l ˜ Z m i (cid:104) l e m | (1 i + m + i − m ) e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + = (cid:88) kk (cid:48) ee (cid:48) e (cid:48)(cid:48) Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | ke (cid:48) m (cid:105) i × N (cid:88) m = m + (cid:88) l e i (cid:104) ke (cid:48) m | l e m (cid:105) i w m e − β E m l ˜ Z m i (cid:104) l e m | k (cid:48) e (cid:48)(cid:48) m (cid:105) i × i (cid:104) k (cid:48) e (cid:48)(cid:48) m | e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + , (41) ρ rs ( m , ˜ τ p ) = (cid:88) lee Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | le m (cid:105) i w m e − β E ml ˜ Z m × i (cid:104) le m | e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + , (42)and, ρ −− rs ( m , ˜ τ p ) = m − (cid:88) m = m (cid:88) l e e Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | l e m (cid:105) i × w m e − β E m l ˜ Z m i (cid:104) l e m | e iH Q τ ... e iH Qp τ p | sem (cid:105) Q p + = m − (cid:88) m = m (cid:88) l e e Q p + (cid:104) rem | (1 Q p + + m + Q p + − m ) e − iH Qp τ p ... e − iH Q τ | l e m (cid:105) i × w m e − β E m l ˜ Z m i (cid:104) l e m | e iH Q τ ... e iH Qp τ p (1 Q p + + m + Q p + − m ) | sem (cid:105) Q p + = (cid:88) kk (cid:48) ee (cid:48) e (cid:48)(cid:48) m − (cid:88) m = m (cid:88) l e Q p + (cid:104) rem | ke (cid:48) m (cid:105) Q p + × Q p + (cid:104) ke (cid:48) m | e − iH Qp τ p ... e − iH Q τ | l e m (cid:105) i w m e − β E m l ˜ Z m × i (cid:104) l e m | e iH Q τ ... e iH Qp τ p | k (cid:48) e (cid:48)(cid:48) m (cid:105) Q p + Q p + (cid:104) k (cid:48) e (cid:48)(cid:48) m | sem (cid:105) Q p + , (43)in which we used 1 i − m | l e m (cid:105) i = m > m ,1 Q p + − m | sem (cid:105) Q p + = m < m , and 1 + m = (cid:80) ke | kem (cid:105)(cid:104) kem | . Noting that the matrix element, Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | se (cid:48) m (cid:105) i , appears in all threeterms above, we denote this matrix element by S mr Qp + s i ( ee (cid:48) , − ˜ τ p ) = Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | se (cid:48) m (cid:105) i , (44) and call it the generalized overlap matrix element. We canprove, see Appendix E, that these generalized overlap matrixelements are diagonal in the environment variables, S mr Qp + s i ( ee (cid:48) , − ˜ τ p ) = S mr Qp + s i ( − ˜ τ p ) δ ee (cid:48) , (45)with S mr Qp + s i ( − ˜ τ p ) being decomposed into three terms S mr Qp + s i ( − ˜ τ p ) = S m ++ r Qp + s i ( − ˜ τ p ) + S m r Qp + s i ( − ˜ τ p ) + S m −− r Qp + s i ( − ˜ τ p )(46)These terms are determined by the following equations, S m ++ r Qp + s i ( − ˜ τ p ) = (cid:88) k S mr Qp + k Qp e − iE mk τ p S mk Qp s i ( − ˜ τ p − ) (47) S m r Qp + s i ( − ˜ τ p ) = (cid:88) l S mr Qp + l Qp e − iE ml τ p S ml Qp s i ( − ˜ τ p − ) (48) S m −− r Qp + s i ( − ˜ τ p ) = (cid:88) α m (cid:88) kk (cid:48) A α m † rk (cid:104) S ( m − k Qp + k (cid:48) i ( − ˜ τ p ) + S ( m − −− k Qp + k (cid:48) i ( − ˜ τ p ) (cid:105) A α m k (cid:48) s with S m −− r Qp + s i ( − ˜ τ p ) = . (49)In the above, and in subsequent equations, indices such as k Q p + k (cid:48) i appearing in full within summands, are abbreviated to kk (cid:48) in summation subscripts (as in (cid:80) kk (cid:48) above). The recursionrelations given by Eqs. (45)-(49) allow the generalized over-lap matrix elements in Eq. (45) to be calculated at an arbitrarytime step, ˜ τ p > ˜ τ . While S m ++ and S m depend on the gen-eralized overlap matrix element at one earlier time step, S m −− depends on the 0 and −− components of these matrix elementsat the same time step, but at one earlier shell. For a detailedproof of the above equations as well as the closed form of thegeneralized overlap matrix elements at ˜ τ we refer the readerto Appendix E.Returning to the terms in the projected density matrix inEqs. (41)-(43), we substitute the generalized overlap matrixelement, which is diagonal in the environment variable, intoeach term as follows˜ ρ ++ rs ( m , ˜ τ p ) = (cid:88) kk (cid:48) S mr Qp + k i ( − ˜ τ p ) R m red ( k , k (cid:48) ) S mk (cid:48) i s Qp + (˜ τ p ) (50) ρ rs ( m , ˜ τ p ) = (cid:88) l S mr Qp + l i ( − ˜ τ p ) w m e − β E ml Z m S ml i s Qp + (˜ τ p ) (51) ρ −− rs ( m , ˜ τ p ) = (cid:88) kk (cid:48) e m − (cid:88) m = m (cid:88) l e Q p + (cid:104) rem | ke m (cid:105) Q p + S m k Qp + l i ( − ˜ τ p ) × w m e − β E m l ˜ Z m S m l i k (cid:48) Qp + (˜ τ p ) Q p + (cid:104) k (cid:48) e m | sem (cid:105) Q p + = (cid:88) kk (cid:48) e m − (cid:88) m = m (cid:88) e Q p + (cid:104) rem | ke m (cid:105) Q p + d N − m ρ kk (cid:48) ( m , ˜ τ p ) × Q p + (cid:104) k (cid:48) e m | sem (cid:105) Q p + = d (cid:88) kk (cid:48) α m A α m † rk (cid:26) ρ kk (cid:48) ( m − , ˜ τ p ) + ρ −− kk (cid:48) ( m − , ˜ τ p ) (cid:27) A α m k (cid:48) s with ρ −− rs ( m , ˜ τ p ) = , (52)5in which S mr i s Qp + (˜ τ p ) = [ S ms Qp + r i ( − ˜ τ p )] † , R m red ( k , k (cid:48) ) is definedas in Eq. (19), and the recursion relation for ρ −− rs ( m , ˜ τ p ) inEq. (52) is derived as for ρ −− rs ( m ) in section IV. The aboveequations determine the projected density matrix of Eq. (40)at an arbitrary time step ˜ τ i ≥ ˜ τ .One can use the scheme in previous works to calcu-late the projected density matrix for multiple quenches, whichyield a simplified ρ i → Q p + rs ( m , ˜ τ p ) = ˜ ρ ++ rs ( m , ˜ τ p ) + ρ rs ( m , ˜ τ p )since ρ −− rs ( m , ˜ τ p ) = S m −− r i s Qp + (˜ τ p ) of the generalizedoverlap matrix element, presented above, then S mr i s Qp + (˜ τ p ) ≈S m ++ r i s Qp + (˜ τ p ) + S m r i s Qp + (˜ τ p ). Using the definitions in Eqs. (47)-(48), we have ρ i → Q p + rs ( m , ˜ τ p ) ≈ (cid:88) qq (cid:48) S mr Qp + q Qp e − i ( E mq − E mq (cid:48) ) τ p ρ i → Q p qq (cid:48) ( m , ˜ τ p − ) S mq (cid:48) Qp s Qp + . This equation is similar to Eqs. (56)-(57) in the recent workof Ref. 51, which investigates time evolution due to periodicswitching. However, S m −− r i s Qp + (˜ τ p ) equals zero only at m = m .Omitting this term leads to an accumulated error after eachtime step due to the causality in Eqs. (45)-(49), where gener-alized overlap matrix elements at a given time step are shownto depend on those at one earlier time step. This omission maybe equivalent to the approximation in the mentioned work,where an accumulated error is observed in the time evolution.In general, by following the procedure for calculating theprojected density matrix for the single quench case in Sec. III,and especially the term ρ −− rs ( m ) in Sec. IV, we explicitlyformulated all the terms of the projected density matrix inEqs. (50)-(52) for the multiple quench case. Substitutingthis projected density matrix into Eq. (38), one can estimatethe time evolution of a local observable at an arbitrary time, t ≥ ˜ τ . For 0 = ˜ τ ≤ t < ˜ τ , one recovers the TDNRGwith a single quench as in Secs II-IV. The entire procedure ofgeneralizing the TDNRG to an arbitrary number of quenches,presented here, adopts no further approximation than does theTDNRG with a single quench. Therefore, by ensuring su ffi -ciently small quantum quenches, one can expect, by the ar-guments presented in the introduction to this section, to ob-tain high accuracy for the time evolution even in the limit oflong times. The numerical implementation of this lies out-side the scope of the present paper and will be reported else-where, however, a few remarks about the computational costof such an approach are in order. In comparison to the singlequench case, the main increase in computational cost comesfrom, (a), the requirement to calculate the projected densitymatrix ρ i → Q p + sr ( m , ˜ τ p ) for p = , . . . , n , i.e. for n quenchesinstead of a single quench, and, (b), the requirement to cal-culate the generalized overlap matrix elements, S mr Qp + s i ( − ˜ τ p ),also for n quenches, using an algorithm which resembles thatused for the calculation of reduced density matrices. Sincethe latter are implemented highly e ffi ciently, we estimate thecomputational cost to be approximately linear in the numberof quenches and therefore feasible. VII. CONCLUSIONS AND OUTLOOK
In this paper, we generalized the TDNRG method for de-scribing the time evolution of a local observable following asingle quantum quench to arbitrary finite temperature, startingfrom the full density matrix of the initial state. In this gener-alization, we made no further approximation than the originalNRG approximation [Eq. (10)]. The generalization relies ondetermining the projected density matrix, which is made up ofthree parts ˜ ρ ++ , ρ −− and ρ , for which explicit expressions, andrecursion relations for their evaluation, were obtained. Thetrace conserving property of the projected density matrix wasproven and was used to check the numerical precision of thecalculations. We find that the contribution from ρ −− , absentin previous approaches, is significant at finite temperatureand cannot be neglected. The short-time limit of local observ-ables was shown to be exactly equal to the thermodynamicvalue in the initial state, both analytically and numerically(Sec. V D). Furthermore, we clarified that the NRG approxi-mation is the main source of noise in local observables at longtimes (Sec. V C), whereas it appears not responsible for theerror in the long-time limit. The latter is found to be reducedby either reducing Λ (i.e. approaching the continuum limit Λ → + ), as noted previously , or by reducing the size of thequantum quench (Sec. V E). Also, the short-time behavior isimproved with decreasing quench size, see Fig. 9.The above results suggested, since the limit Λ → + isimpractical in NRG calculations, a formulation in which alarge quench is replaced by multiple quenches over a finitetime interval as a means to improve the TDNRG at all timescales, and particularly for longer times. We provided thismultiple quench formulation (Sec. VI) and showed that, asfor the single quench case, it rests solely on a single approx-imation (the NRG approximation). The structure of this the-ory, while formally similar to the single quench case, bringsout more clearly the importance of including the “ −− ” terms,i.e., ρ −− and S m −− . Namely, neglecting S m −− in the multiplequench case leads to an accumulated error after each time step˜ τ p , p = , . . . , n due to the causality of Eqs. (45)-(49), hence,it will be particularly important not to approximate this in afuture numerical implementation of the multiple-quench TD-NRG. The latter, as discussed at the end of the previous sec-tion is feasible, scaling at most linearly with the number ofquenches, but lies beyond the scope of the present paper andwill be presented elsewhere. Our formalism for the multiplequench case also generalizes the TDNRG to general pulsesand periodic switching / driving. The latter has been consideredalso in Ref. 51, within a hybrid TDNRG DMRG approach,although several approximations in addition to the NRG ap-proximation were made (see the discussion in Sec. VI).We also compared our single-quench TDNRG approach tothe previous scheme, finding that the present one gives re-sults with less noise at long times, with the results being closerto the expected long-time limit, particularly at finite tempera-tures. This is due to the contribution of the full Wilson chainin our TDNRG using the full density matrix of the initial state,which at finite temperature always results in a finite contribu-tion from ρ −− , absent in the scheme of Refs. 14 and 40 which6uses a truncated Wilson chain for finite temperatures.In future, it would be of interest, especially in the lightof time-resolved experiments, e.g., time-resolved photoemis-sion and scanning tunneling microscope spectroscopies, to extend the formalism in this paper to the time-evolutionof dynamical quantities following a quench / pulse, e.g.,to single-particle spectral functions A ( ω, t ), dynami-cal spin susceptibilities, or optical conductivities. Thiswould allow basic questions about the time-dependence ofthe Kondo e ff ect in quantum dots to be addressed, e.g. thetime-evolution of the Kondo resonance upon instantaneouslyswitching from the mixed-valence to the Kondo regime, or,the transient response of the current following a bias voltagepulse of finite duration. Given the common matrix productstate structure of the NRG and DMRG eigenstates, someof the formalism developed in this paper for general mul-tiple quenches might find application also within the time-dependent DMRG. Finally, improved TDNRG approacheswould be of interest in the development of methods for study-ing steady-state non-equilibrium transport through correlatedsystems, such as quantum dots. ACKNOWLEDGMENTS
We thank D. P. DiVincenzo, S. Andergassen, V. Meden, A.Weichselbaum and F. B. Anders for useful discussions and en-couraging remarks, and acknowledge supercomputer supportby the John von Neumann institute for Computing (J¨ulich).
Appendix A: Trace of the Projected Density Matrix
Starting from the exact Eq. (8) and setting ˆ O = = Tr[ ρ ], a tautology, since ρ is normalized. However,the result we are interested in starts from Eq. (9), which ingeneral will be approximate due to the use of the NRG ap-proximation. Inserting ˆ O = O mrs = δ rs results in Eq. (12)which states that the trace of the projected density matrix ispreserved. The result is exact, for the same reason that theshort-time limit is exact (as discussed in Sec. II D), namelyany error that would arise from approximate NRG eigenval-ues in Eq. (9), does not arise for ˆ O = ρ = (cid:80) Nm = m w m ˜ ρ m , from Sec. II B and the completeness re-lation (cid:80) Nm = m (cid:80) le | lem (cid:105)(cid:104) lem | =
1, we have for the RHS of Eq. (12), N (cid:88) m = m (cid:88) l ρ i → fll ( m ) ≡ N (cid:88) m = m (cid:88) l , e f (cid:104) lem | ρ | lem (cid:105) f = N (cid:88) m = m N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) l , e (cid:88) l (cid:48) , e (cid:48) f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) i (cid:104) l (cid:48) e (cid:48) m (cid:48) | lem (cid:105) f = N (cid:88) m = m N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) l , e (cid:88) l (cid:48) , e (cid:48) i (cid:104) l (cid:48) e (cid:48) m (cid:48) | lem (cid:105) f f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) = N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) l (cid:48) , e (cid:48) i (cid:104) l (cid:48) e (cid:48) m (cid:48) | l (cid:48) e (cid:48) m (cid:48) (cid:105) i e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) = N (cid:88) m (cid:48) = m w m (cid:48) (cid:88) l (cid:48) e − β E m (cid:48) l (cid:48) Z m (cid:48) = , (A1)where, in the last line, (cid:80) Nm (cid:48) = m w m (cid:48) ≡ Tr[ ρ ] = ρ i → f ( m ) in Eq. (24)have been evaluated correctly, and, (b), the numerical calcula-tion of ρ i → f ( m ) is su ffi cient. This was the case, since we foundan error of less than 10 − for all Λ and all system parametersused in the Anderson model ( ε d , U ); see Fig. 1. Appendix B: Short-Time Limit
In this appendix we show explicitly that the TDNRG ex-pression for the expectation value of a local observable in theshort-time limit ( t → + ) recovers the correct thermodynamicvalue in the initial state, i.e. that O ( t → + ) ≡ N (cid:88) m = m (cid:88) rs (cid:60) KK (cid:48) ρ i → fsr ( m ) O mrs = O i . (B1)where O i = Tr (cid:104) ρ ˆ O (cid:105) .Similar to the calculation of the trace of the projected den-sity matrix in Eq. (A1), we substitute ρ = (cid:80) Nm = m w m ˜ ρ m intoEq. (B1), obtaining O ( t → + ) = N (cid:88) m = m (cid:44) KK (cid:48) (cid:88) rse f (cid:104) rem | ρ | sem (cid:105) f f (cid:104) sem | ˆ O | rem (cid:105) f = N (cid:88) m = m (cid:44) KK (cid:48) (cid:88) rse N (cid:88) m (cid:48) = m (cid:88) l (cid:48) e (cid:48) f (cid:104) rem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i w m (cid:48) e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) × i (cid:104) l (cid:48) e (cid:48) m (cid:48) | sem (cid:105) f f (cid:104) sem | ˆ O | rem (cid:105) f = N (cid:88) m (cid:48) = m (cid:88) l (cid:48) e (cid:48) w m (cid:48) e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) (cid:110) N (cid:88) m = m (cid:44) KK (cid:48) (cid:88) rse i (cid:104) l (cid:48) e (cid:48) m (cid:48) | sem (cid:105) f × f (cid:104) sem | ˆ O | rem (cid:105) f f (cid:104) rem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i (cid:111) . (B2)Decomposing | rem (cid:105) f into kept, | kem (cid:105) f , and discarded, | lem (cid:105) f ,7states, we obtain for the part inside the curly brackets N (cid:88) m = m (cid:88) lse i (cid:104) l (cid:48) e (cid:48) m (cid:48) | sem (cid:105) f f (cid:104) sem | ˆ O | lem (cid:105) f f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i + N (cid:88) m = m (cid:88) kle i (cid:104) l (cid:48) e (cid:48) m (cid:48) | lem (cid:105) f f (cid:104) lem | ˆ O | kem (cid:105) f f (cid:104) kem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i = N (cid:88) m = m N (cid:88) m = m (cid:88) le (cid:88) l e i (cid:104) l (cid:48) e (cid:48) m (cid:48) | l e m (cid:105) f f (cid:104) l e m | ˆ O | lem (cid:105) f × f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i + N (cid:88) m = m N (cid:88) m = m + (cid:88) l e (cid:88) le i (cid:104) l (cid:48) e (cid:48) m (cid:48) | lem (cid:105) f f (cid:104) lem | ˆ O | l e m (cid:105) f × f (cid:104) l e m | l (cid:48) e (cid:48) m (cid:48) (cid:105) i . (B3)Replacing (cid:80) Nm = m (cid:80) Nm = m → (cid:80) Nm = m (cid:80) m m = m in the first term,and interchanging | lem (cid:105) ↔ | l e m (cid:105) in the second one, wehave Eq. (B3) equal to N (cid:88) m = m m (cid:88) m = m (cid:88) le (cid:88) l e i (cid:104) l (cid:48) e (cid:48) m (cid:48) | l e m (cid:105) f f (cid:104) l e m | ˆ O | lem (cid:105) f × f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i + N (cid:88) m = m N (cid:88) m = m + (cid:88) le (cid:88) l e i (cid:104) l (cid:48) e (cid:48) m (cid:48) | l e m (cid:105) f f (cid:104) l e m | ˆ O | lem (cid:105) f × f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i = N (cid:88) m = m N (cid:88) m = m (cid:88) le (cid:88) l e i (cid:104) l (cid:48) e (cid:48) m (cid:48) | l e m (cid:105) f f (cid:104) l e m | ˆ O | lem (cid:105) f × f (cid:104) lem | l (cid:48) e (cid:48) m (cid:48) (cid:105) i = i (cid:104) l (cid:48) e (cid:48) m (cid:48) | ˆ O | l (cid:48) e (cid:48) m (cid:48) (cid:105) i . (B4)Substituting Eq. (B4) into Eq. (B2), we have O ( t → + ) = N (cid:88) m (cid:48) = m (cid:88) l (cid:48) e (cid:48) w m (cid:48) e − β E m (cid:48) l (cid:48) ˜ Z m (cid:48) i (cid:104) l (cid:48) e (cid:48) m (cid:48) | ˆ O | l (cid:48) e (cid:48) m (cid:48) (cid:105) i = N (cid:88) m (cid:48) = m (cid:88) l (cid:48) w m (cid:48) e − β E m (cid:48) l (cid:48) Z m (cid:48) O l (cid:48) l (cid:48) , (B5)which is the thermodynamic average of the observable in theinitial state, as calculated within the FDM approach. Thisproves that O ( t → + ) = O i is always true. We also verifiedthe latter numerically, for the Anderson impurity model, find-ing that the calculated value of O ( t → + ) lies within 10 − of O i regardless of the parameter sets used and regardless of thediscretization parameter Λ . Appendix C: Recursion Relations for SU(2) Symmetry
In the actual calculations for the Anderson impurity model,presented in this paper, SU(2) spin symmetry was used in the numerical diagonalizations. Accordingly, the states were clas-sified according to their values of total charge Q , total spin S ,and z component of total spin S z . Moreover, as the eigenval-ues are independent of S z , we need only write the recursionrelations in the ( Q , S ) subspace. The transformation matrix A α m kr in Eq. (28) is now denoted by U m ( ki , r ) as in the originalwork in Ref. 2, with the correspondence i ↔ α m . In detail, ρ −− sr ( m ) is expressed in terms of the following four terms, cor-responding to i taking one of four values ( d = ρ −− sr ( Q , S , m ) = (cid:88) k , k (cid:48) U m † Q , S ( k , s ) ˜ ρ k , k (cid:48) ( Q + , S , m − U mQ , S ( k (cid:48) , r ) + (cid:88) k , k (cid:48) U m † Q , S ( k , s ) ˜ ρ k , k (cid:48) ( Q , S − , m − U mQ , S ( k (cid:48) , r ) + (cid:88) k , k (cid:48) U m † Q , S ( k , s ) ˜ ρ k , k (cid:48) ( Q , S + , m − U mQ , S ( k (cid:48) , r ) + (cid:88) k , k (cid:48) U m † Q , S ( k , s ) ˜ ρ k , k (cid:48) ( Q − , S , m − U mQ , S ( k (cid:48) , r ) , (C1)with ˜ ρ k , k (cid:48) ( Q , S , m ) = ρ k , k (cid:48) ( Q , S , m ) + ρ −− k , k (cid:48) ( Q , S , m ).Similarly, the recursion relation for the reduced full densitymatrices takes the form[ R m red ( Q , S )] k , k (cid:48) = (cid:88) q , q (cid:48) U m + Q − , S ( k , q )[ ˜ R m + ( Q − , S )] q , q (cid:48) U m + † Q − , S ( k (cid:48) , q (cid:48) ) + S + S + (cid:88) q , q (cid:48) U m + Q , S + ( k , q )[ ˜ R m + ( Q , S +
12 )] q , q (cid:48) U m + † Q , S + ( k (cid:48) , q (cid:48) ) + S S + (cid:88) q , q (cid:48) U m + Q , S − ( k , q )[ ˜ R m + ( Q , S −
12 )] q , q (cid:48) U m + † Q , S − ( k (cid:48) , q (cid:48) ) + (cid:88) q , q (cid:48) U m + Q + , S ( k , q )[ ˜ R m + ( Q + , S )] q , q (cid:48) U m + † Q + , S ( k (cid:48) , q (cid:48) ) . (C2)with (cid:104) ˜ R m red ( Q , S ) (cid:105) = (cid:104) R m red ( Q , S ) (cid:105) k , k (cid:48) · · · w m e − β EmQ , S , l Z m · · · ... ... . . . ... · · · w m e − β EmQ , S , l Z m ,and noting that (cid:80) q = (cid:80) k + (cid:80) l and | qm (cid:105) = | km (cid:105) + | lm (cid:105) . Appendix D: Long-Time Limit: Additional Switching Protocols
In this appendix, we present additional results for the longtime limit of observables with di ff erent switching protocolsto support the conclusions in Sec. V E. In general, the resultsshow us the trend of quench size dependence as observed inFig. 5. Depending on the local observable, i.e., local leveloccupation number or double occupancy, the results dependon di ff erent parameters of the final state.8
1. Double occupancy: switching ε d with U constant -4 -3 -2 -1 δ 〈 n d ↑ n d ↓ 〉 ( t → + ∞ ) T/T K (b) ε i =-6 Γε f = 0 Γε f =-1 Γε f =-2 Γε f =-3 Γε f =-4 Γε f =-5 Γ -5 0 5 10 (a) ε f =-6 Γε i = 0 Γε i =-1 Γε i =-2 Γε i =-3 Γε i =-4 Γε i =-5 Γ FIG. 11. (Color online) The (percentage) relative error of thedouble occupation number in the long time limit vs temperatures. U i = U f = U = Γ and Γ = − D . (a) Switching from the asym-metric to the symmetric Kondo system. (b) Reverse switching, fromthe symmetric to the asymmetric Kondo system [the ε f are indicatedand are the same as the ε i in (a)]. T K ≈ . × − D is the Kondotemperature of the symmetric model and calculations are for Λ = . z averaging. In Fig. 11, we present the error in the double occupancyupon switching the system between the asymmetric case andthe symmetric Kondo regime while keeping the Coulomb in-teraction constant, i.e. ε d ( t ) = ε i θ ( − t ) + ε f θ ( t ) with U = U i = U f . The analogous results for the occupation number werepresented in Fig. 5. As for the error in the occupation numberin Fig. 5, we observe for that the error in the double occupancyexhibits extrema at approximately the same temperatures,namely at T ≈ T K ≈ U / T ≈ T K = . Γ . The errorincreases with increasing quench size for both the switchingfrom the asymmetric model to the symmetric Kondo regime[Fig. 11 (a)] and vice versa [Fig. 11 (b)], supporting the firsttrend discussed in Sec. V E. The highest incoherent excitationin Fig. 11 (a), defined by Eq. (35), is 6 Γ , while in Fig. 11 (b),it varies from 7 Γ up to 12 Γ . For the same sized quench inFigs. 11 (a) and 11(b), we see that the error is always larger for Fig. 11 (b), thus supporting the conclusion that this errordepends on the magnitude of the highest incoherent excitationin the final state.
2. Occupation number and double occupancy: switching U with ε d constant -15-10-5 010 -4 -3 -2 -1 δ n d ( t → + ∞ ) T/T K (b) U i =12 Γ U f = 0 Γ U f = 2 Γ U f = 4 Γ U f = 6 Γ U f = 8 Γ U f =10 Γ f =12 Γ U i = 0 Γ U i = 2 Γ U i = 4 Γ U i = 6 Γ U i = 8 Γ U i =10 Γ FIG. 12. (Color online) The (percentage) relative error of the lo-cal level occupation number in the long time limit vs temperature. ε i = ε f = − Γ and Γ = − D . (a) Switching from the asymmetricmodel to the symmetric Kondo regime. (b) Reverse switching, fromthe symmetric Kondo regime to the asymmetric model [the ε f areindicated and are the same as the ε i in (a)]. T K ≈ . × − D is theKondo temperature of the symmetric model and calculations are for Λ = . z averaging. In Sec. V E we discussed the error in the long time limit ofthe double occupancy for the case where both U and ε d wereswitched at t =
0, while maintaining particle-hole symme-try, i.e. for ε d ( t ) + U ( t ) / =
0. Here, we consider relaxing thelast condition, and consider changing the Coulomb interaction U ( t ) = U i θ ( − t ) + U f θ ( t ) while keeping ε d = ε i = ε f = − Γ constant in both initial and final states. We either choose U i = Γ and switch U f , or vice versa. Figure 12 showsthe error in the long time limit of the occupation number andFig. 13 the corresponding error for the double occupancy. Forboth quantities, the error increases with increasing size of thequantum quench and the positions of the extrema in the er-9 -30-20-10 010 -4 -3 -2 -1 δ 〈 n d ↑ n d ↓ 〉 ( t → + ∞ ) T/T K (b) U i =12 Γ U f = 0 Γ U f = 2 Γ U f = 4 Γ U f = 6 Γ U f = 8 Γ U f =10 Γ f =12 Γ U i = 0 Γ U i = 2 Γ U i = 4 Γ U i = 6 Γ U i = 8 Γ U i =10 Γ FIG. 13. (Color online) The (percentage) relative error of the dou-ble occupancy in the long time limit vs temperatures. Parametersand switching in (a) and (b) are respectively the same as those inFig. 12 (a) and (b). ror are consistent with those seen in previous quenches. Thetrend that the error increases with increasing value of the high-est incoherent excitation in the final state [Eq. (35)] can notbe applied in this case as the value of this excitation is thesame (6 Γ ) for both switchings and all cases. However, Fig. 12agrees with the trend in Fig. 5 that switching from the sym-metric to the asymmetric model results in a larger error thanvice versa, and Fig. 13 agrees with the trend in Fig. 6 that alarger error is expected upon switching from a less correlatedto a more correlated system than vice versa. Appendix E: Generalized overlap matrix elements
In this appendix, we prove that the generalized overlap ma-trix elements at an arbitrary time step are diagonal in the envi-ronment variables, and that they can be calculated recursively.We start with the simplest case for the first time step ˜ τ S mr Q s i ( ee (cid:48) , − ˜ τ ) = Q (cid:104) rem | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) m l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ | se (cid:48) m (cid:105) i . Decomposing the matrix element into three terms, corre- sponding to m > m , m = m , and m < m , contributions,respectively, we have S mr Q s i ( ee (cid:48) , − ˜ τ ) = S m ++ r Q s i ( ee (cid:48) , − ˜ τ ) + S m r Q s i ( ee (cid:48) , − ˜ τ ) + S m −− r Q s i ( ee (cid:48) , − ˜ τ ) , with each term shown to be diagonal in the environment vari-ables as follows S m ++ r Q s i ( ee (cid:48) , − ˜ τ ) = N (cid:88) m = m + (cid:88) l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) k Q (cid:104) rem | kem (cid:105) Q Q (cid:104) kem | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) k S mr Q k Q e − iE mk τ S mk Q s i δ ee (cid:48) = S m ++ r Q s i ( − ˜ τ ) δ ee (cid:48) (E1) S m r Q s i ( ee (cid:48) , − ˜ τ ) = (cid:88) l Q (cid:104) rem | lem (cid:105) Q Q (cid:104) lem | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) l S mr Q l Q e − iE ml τ S ml Q s i δ ee (cid:48) = S m r Q s i ( − ˜ τ ) δ ee (cid:48) (E2) S m −− r Q s i ( ee (cid:48) , − ˜ τ ) = m − (cid:88) m = m (cid:88) l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) l e Q (cid:104) rem | (1 Q + m + Q − m ) | l e m (cid:105) Q × Q (cid:104) l e m | e − iH Q τ (1 i + m + i − m ) | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) kk (cid:48) e e (cid:48) Q (cid:104) rem | ke m (cid:105) Q (cid:88) l Q (cid:104) ke m | l e m (cid:105) Q × Q (cid:104) l e m | e − iH Q τ | k (cid:48) e (cid:48) m (cid:105) ii (cid:104) k (cid:48) e (cid:48) m | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) kk (cid:48) e Q (cid:104) rem | ke m (cid:105) Q × (cid:34) (cid:88) l S m k Q l Q e − iE m l τ S m l Q k (cid:48) i (cid:35) i (cid:104) k (cid:48) e m | se (cid:48) m (cid:105) i . Using the definition of S m r Q s i ( − ˜ τ ) in Eq. (E2), we have S m −− r Q s i ( ee (cid:48) , − ˜ τ ) = m − (cid:88) m = m (cid:88) kk (cid:48) e Q (cid:104) rem | ke m (cid:105) Q (cid:104) S m k Q k (cid:48) i ( − ˜ τ ) (cid:105) i (cid:104) k (cid:48) e m | se (cid:48) m (cid:105) i . Noticing that the structure of the above equation is simi-lar to that of Eq. (22), we derive the recursion relation of S m −− r Q s i ( ee (cid:48) , − ˜ τ ) following the procedure for ρ −− in sec. IV,0thus we have S m −− r Q s i ( ee (cid:48) , − ˜ τ ) = (cid:88) α m (cid:88) kk (cid:48) A α m † rk (cid:104) S ( m − k Q k (cid:48) i ( − ˜ τ ) + S ( m − −− k Q k (cid:48) i ( − ˜ τ ) (cid:105) A α m k (cid:48) s δ ee (cid:48) = S m −− r Q s i ( − ˜ τ ) δ ee (cid:48) w ith S m −− r Q s i ( − ˜ τ ) = . (E3)Since each term of the generalized overlap matrix element isdiagonal in the environment variables, we define S mr Q s i ( ee (cid:48) , − ˜ τ ) = S mr Q s i ( − ˜ τ ) δ ee (cid:48) (E4)with S mr Q s i ( − ˜ τ ) = S m ++ r Q s i ( − ˜ τ ) + S m r Q s i ( − ˜ τ ) + S m −− r Q s i ( − ˜ τ ) . In a similar way, we consider the generalized overlap matrixelement at the next time step ˜ τ , S mr Q s i ( ee (cid:48) , − ˜ τ ) = Q (cid:104) rem | e − iH Q τ e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) m l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ e − iH Q τ | se (cid:48) m (cid:105) i = S m ++ r Q s i ( ee (cid:48) , − ˜ τ ) + S m r Q s i ( ee (cid:48) , − ˜ τ ) + S m −− r Q s i ( ee (cid:48) , − ˜ τ ) . (E5)The three terms above correspond to m > m , m = m , and m < m , respectively, and are shown to be diagonal in theenvironment variables as follows, S m ++ r Q s i ( ee (cid:48) , − ˜ τ ) = N (cid:88) m = m + (cid:88) l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) k Q (cid:104) rem | kem (cid:105) Q Q (cid:104) kem | e − iH Q τ e − iH Q τ | sem (cid:105) i = (cid:88) k Q (cid:104) rem | kem (cid:105) Q e − iE mk τ Q (cid:104) kem | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) k S mr Q k Q e − iE mk τ S mk Q s i ( − ˜ τ ) δ ee (cid:48) = S m ++ r Q s i ( − ˜ τ ) δ ee (cid:48) (E6) S m r Q s i ( ee (cid:48) , − ˜ τ ) = (cid:88) l Q (cid:104) rem | lem (cid:105) Q Q (cid:104) lem | e − iH Q τ e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) l Q (cid:104) rem | lem (cid:105) Q e − iE ml τ Q (cid:104) lem | e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) l S mr Q l Q e − iE ml τ S ml Q s i ( − ˜ τ ) δ ee (cid:48) = S m r Q s i ( − ˜ τ ) δ ee (cid:48) (E7) S m −− r Q s i ( ee (cid:48) , − ˜ τ ) = m − (cid:88) m = m (cid:88) l e Q (cid:104) rem | l e m (cid:105) Q Q (cid:104) l e m | e − iH Q τ e − iH Q τ | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) kk (cid:48) e e (cid:48) Q (cid:104) rem | ke m (cid:105) Q (cid:88) l Q (cid:104) ke m | l e m (cid:105) Q × Q (cid:104) l e m | e − iH Q τ e − iH Q τ | k (cid:48) e (cid:48) m (cid:105) ii (cid:104) k (cid:48) e (cid:48) m | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) kk (cid:48) e Q (cid:104) rem | ke m (cid:105) Q (cid:34) (cid:88) l S m k Q l Q e − iE m l τ S l Q k (cid:48) i ( m , − ˜ τ ) (cid:35) × i (cid:104) k (cid:48) e m | se (cid:48) m (cid:105) i = m − (cid:88) m = m (cid:88) kk (cid:48) e Q (cid:104) rem | ke m (cid:105) Q (cid:104) S m k Q k (cid:48) i ( − ˜ τ ) (cid:105) i (cid:104) k (cid:48) e m | se (cid:48) m (cid:105) i . Again, noticing that the structure of the last expression aboveis similar to that in Eq. (22), we derive the following recursionrelation S m −− r Q s i ( ee (cid:48) , − ˜ τ ) = (cid:88) α m (cid:88) kk (cid:48) A α m † rk (cid:104) S ( m − k Q k (cid:48) i ( − ˜ τ ) + S ( m − −− k Q k (cid:48) i ( − ˜ τ ) (cid:105) A α m rs δ ee (cid:48) = S m −− r Q s i ( − ˜ τ ) δ ee (cid:48) with S m −− r Q s i ( − ˜ τ ) = . (E8)Since each term of the overlap matrix element at ˜ τ is alsodiagonal in the environment variables, we have S mr Q s i ( ee (cid:48) , − ˜ τ ) = S mr Q s i ( − ˜ τ ) δ ee (cid:48) (E9)with S mr Q s i ( − ˜ τ ) = S m ++ r Q s i ( − ˜ τ ) + S m r Q s i ( − ˜ τ ) + S m −− r Q s i ( − ˜ τ ) . Consequently, the generalized overlap matrix elements at anarbitrary time step ˜ τ p + > ˜ τ can decomposed as follows S mr Qp + s i ( ee (cid:48) , − ˜ τ p ) = Q p + (cid:104) rem | e − iH Qp τ p ... e − iH Q τ | se (cid:48) m (cid:105) i = (cid:88) m p l p e p Q p + (cid:104) rem | l p e p m p (cid:105) Q p Q p (cid:104) l p e p m p | e − iH Qp τ ... e − iH Q τ | se (cid:48) m (cid:105) i = S m ++ r Qp + s i ( ee (cid:48) , − ˜ τ p ) + S m r Qp + s i ( ee (cid:48) , − ˜ τ p ) + S m −− r Qp + s i ( ee (cid:48) , − ˜ τ p ) . Again, the three terms above correspond to the m > m p , m = m p , and m < m p contributions, and satisfy the general rule S m ++ r Qp + s i ( ee (cid:48) , − ˜ τ p ) = (cid:88) k S mr Qp + k Qp e − iE mk τ p S mk Qp s i ( − ˜ τ p − ) δ ee (cid:48) = S m ++ r Qp + s i ( − ˜ τ p ) δ ee (cid:48) (E10) S m r Qp + s i ( ee (cid:48) , − ˜ τ p ) = (cid:88) l S mr Qp + l Qp e − iE ml τ p S ml Qp s i ( − ˜ τ p − ) δ ee (cid:48) = S m r Qp + s i ( − ˜ τ p ) δ ee (cid:48) (E11) S m −− r Qp + s i ( ee (cid:48) , − ˜ τ p ) = (cid:88) α m (cid:88) kk (cid:48) A α m † rk (cid:104) S ( m − k Qp + k (cid:48) i ( − ˜ τ p ) + S ( m − −− k Qp + k (cid:48) i ( − ˜ τ p ) (cid:105) A α m k (cid:48) s δ ee (cid:48) = S m −− r Qp + s i ( − ˜ τ p ) δ ee (cid:48) with S m −− r Qp + s i ( − ˜ τ p ) = . (E12)1Since each term of S mr f s i ( ee (cid:48) , − ˜ τ p ) is diagonal in the environ- ment variables, we have S mr Qp + s i ( ee (cid:48) , − ˜ τ p ) = S mr Qp + s i ( − ˜ τ p ) δ ee (cid:48) (E13)with S mr Qp + s i ( − ˜ τ p ) = S m ++ r Qp + s i ( − ˜ τ p ) + S m r Qp + s i ( − ˜ τ p ) + S m −− r Qp + s i ( − ˜ τ p ) . In general, we showed that the generalized overlap matrixelements at an arbitrary time step are diagonal in the envi-ronment variables. Moreover, one can use Eqs. (E1)-(E3) tocalculate each term of the generalized overlap matrix elementat ˜ τ , and Eqs. (E10)-(E12) to calculate those at ˜ τ p > ˜ τ . K. G. Wilson, Rev. Mod. Phys. , 773 (1975). H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, Phys.Rev. B , 1003 (1980). H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, Phys.Rev. B , 1044 (1980). R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. , 395(2008). A. C. Hewson,
The Kondo Problem to Heavy Fermions (Cam-bridge University Press, Cambridge, 1997). L. N. Oliveira and J. W. Wilkins, Phys. Rev. Lett. , 1553(1981). C. Gonzalez-Buxton and K. Ingersent, Phys. Rev. B , 14254(1998). L. Merker and T. A. Costi, Phys. Rev. B , 075150 (2012). H. O. Frota and L. N. Oliveira, Phys. Rev. B , 7871 (1986). O. Sakai, Y. Shimizu, and T. Kasuya, Journal of the PhysicalSociety of Japan , 3666 (1989). T. A. Costi and A. C. Hewson, Philosophical Magazine Part B , 1165 (1992). R. Bulla, A. C. Hewson, and T. Pruschke, Journal of Physics:Condensed Matter , 8365 (1998). W. Hofstetter, Phys. Rev. Lett. , 1508 (2000). F. B. Anders and A. Schiller, Phys. Rev. Lett. , 196801 (2005). R. Peters, T. Pruschke, and F. B. Anders, Phys. Rev. B ,245114 (2006). A. Weichselbaum and J. von Delft, Phys. Rev. Lett. , 076402(2007). T. A. Costi, A. C. Hewson, and V. Zlati´c, J. Phys.: Condens.Matter , 2519 (1994). M. Misiorny, I. Weymann, and J. Barna´s, Phys. Rev. B ,245415 (2012). M. Hanl, A. Weichselbaum, T. A. Costi, F. Mallet, L. Samina-dayar, C. B¨auerle, and J. von Delft, Phys. Rev. B , 075146(2013). W. Metzner and D. Vollhardt, Phys. Rev. Lett. , 324 (1989). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev.Mod. Phys. , 13 (1996). G. Kotliar and D. Vollhardt, Physics Today , 53 (2004). D. Vollhardt, Annalen der Physik , 1 (2012). T. A. Costi and N. Manini, Journal of Low Temperature Physics , 835 (2002). O. Bodensiek, R. ˇZitko, M. Vojta, M. Jarrell, and T. Pruschke,Phys. Rev. Lett. , 146406 (2013). T. Pruschke, R. Bulla, and M. Jarrell, Phys. Rev. B , 12799(2000). R. Bulla, Phys. Rev. Lett. , 136 (1999). U. Schneider, L. Hackerm¨uller, S. Will, T. Best, I. Bloch, T. A.Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science , 1520 (2008). J. Bauer and A. C. Hewson, Phys. Rev. B , 235113 (2010). C. P. Moca, A. Alex, J. von Delft, and G. Zar´and, Phys. Rev. B , 195128 (2012). A. Weichselbaum, Annals of Physics , 2972 (2012). W. C. Oliveira and L. N. Oliveira, Phys. Rev. B , 11986 (1994). V. L. Campo and L. N. Oliveira, Phys. Rev. B , 104432 (2005). R. ˇZitko and T. Pruschke, Phys. Rev. B , 085106 (2009). A. K. Mitchell, M. R. Galpin, S. Wilson-Fletcher, D. E. Logan,and R. Bulla, ArXiv e-prints (2013), arXiv:1308.1903 [cond-mat.str-el]. A. I. T´oth, C. P. Moca, O. Legeza, and G. Zar´and, Phys. Rev. B , 245109 (2008). M. Yoshida, A. C. Seridonio, and L. N. Oliveira, Phys. Rev. B , 235317 (2009). T. A. Costi and V. Zlati´c, Phys. Rev. B , 235127 (2010). T. A. Costi, Phys. Rev. B , 3003 (1997). F. B. Anders and A. Schiller, Phys. Rev. B , 245113 (2006). F. B. Anders, New Journal of Physics , 115007 (2008). S. Tornow, R. Bulla, F. B. Anders, and A. Nitzan, Phys. Rev. B , 035434 (2008). A. Hackl, D. Roosen, S. Kehrein, and W. Hofstetter, Phys. Rev.Lett. , 196601 (2009). P. P. Orth, D. Roosen, W. Hofstetter, and K. Le Hur, Phys. Rev.B , 144423 (2010). F. B. Anders, Phys. Rev. Lett. , 066804 (2008). S. Schmitt and F. B. Anders, Phys. Rev. Lett. , 056801 (2011). F. B. Anders, Journal of Physics-Condensed Matter , 195216(2008). A. Rosch, The European Physical Journal B-Condensed Matterand Complex Systems , 1 (2012). In fact, we show in this paper that the short-time limit is alwaysrecovered exactly by the td-NRG. P. Schmitteckert, Journal of Physics Conference Series ,012022 (2010). E. Eidelstein, A. Schiller, F. G¨uttge, and F. B. Anders, Phys. Rev.B , 075118 (2012). F. G¨uttge, F. B. Anders, U. Schollw¨ock, E. Eidelstein, andA. Schiller, Phys. Rev. B , 115115 (2013). M. A. Cazalilla and J. B. Marston, Phys. Rev. Lett. , 256403(2002). H. G. Luo, T. Xiang, and X. Q. Wang, Phys. Rev. Lett. ,049701 (2003). A. J. Daley, C. Kollath, U. Schollw¨ock, and G. Vidal, Journalof Statistical Mechanics: Theory and Experiment , P04005(2004). S. R. White and A. E. Feiguin, Phys. Rev. Lett. , 076401(2004). W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, andK. Sch¨onhammer, Rev. Mod. Phys. , 299 (2012). H. Schoeller, The European Physical Journal Special Topics ,179 (2009). A. Rosch, J. Paaske, J. Kroha, and P. W¨olfle, Phys. Rev. Lett. ,076804 (2003). A. Kamenev,
Field Theory of Non-Equilibrium Systems (Cam-bridge University Press, Cambridge, 2011). R. B. Saptsov and M. R. Wegewijs, Phys. Rev. B , 235432(2012). A. C. Hewson, J. Bauer, and A. Oguri, Journal of Physics Con-densed Matter , 5413 (2005). D. Lobaskin and S. Kehrein, Phys. Rev. B , 193303 (2005). M. Moeckel and S. Kehrein, Phys. Rev. Lett. , 175702 (2008). P. Wang and S. Kehrein, Phys. Rev. B , 125124 (2010). C. Jung, A. Lieder, S. Brener, H. Hafermann, B. Baxevanis,A. Chudnovskiy, A. Rubtsov, M. Katsnelson, and A. Lichten-stein, Annalen der Physik , 49 (2012). D. C. Langreth and P. Nordlander, Phys. Rev. B , 2541 (1991). J. Merino and J. B. Marston, Phys. Rev. B , 6982 (1998). Z. Ratiani and A. Mitra, Phys. Rev. B , 245111 (2009). M. Pletyukhov, D. Schuricht, and H. Schoeller, Phys. Rev. Lett. , 106801 (2010). D. M. Kennes, O. Kashuba, M. Pletyukhov, H. Schoeller, andV. Meden, Phys. Rev. Lett. , 100405 (2013). S. Andergassen, M. Pletyukhov, D. Schuricht, H. Schoeller, andL. Borda, Phys. Rev. B , 205103 (2011). C. Karrasch, S. Andergassen, M. Pletyukhov, D. Schuricht,L. Borda, V. Meden, and H. Schoeller, EPL (Europhysics Let-ters) , 30003 (2010). D. M. Kennes, S. G. Jakobs, C. Karrasch, and V. Meden, Phys.Rev. B , 085113 (2012). D. M. Kennes and V. Meden, Phys. Rev. B , 245101 (2012). Y. Meir, N. S. Wingreen, and P. A. Lee, Phys. Rev. Lett. , 2601(1993). H. Shao, D. C. Langreth, and P. Nordlander, Phys. Rev. B ,13929 (1994). H. Shao, D. C. Langreth, and P. Nordlander, Phys. Rev. B ,13948 (1994). A. Goker, B. A. Friedman, and P. Nordlander, Journal of PhysicsCondensed Matter , 376206 (2007). A. Oguri, Phys. Rev. B , 153305 (2001). E. Mu˜noz, C. J. Bolech, and S. Kirchner, Phys. Rev. Lett. ,016601 (2013). R. B. Saptsov and M. R. Wegewijs, ArXiv e-prints (2013),arXiv:1311.1368 [cond-mat.str-el]. A. E. Feiguin and S. R. White, Phys. Rev. B , 220401 (2005). T. Barthel, U. Schollw¨ock, and S. R. White, Phys. Rev. B ,245101 (2009). C. Karrasch, J. H. Bardarson, and J. E. Moore, New Journal ofPhysics , 083031 (2013). G. Cohen, E. Gull, D. R. Reichman, A. J. Millis, and E. Rabani,Phys. Rev. B , 195108 (2013). G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, ArXive-prints (2013), arXiv:1310.4151 [cond-mat.str-el]. E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer,and P. Werner, Rev. Mod. Phys. , 349 (2011). S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B ,195316 (2008). L. M¨uhlbacher and E. Rabani, Phys. Rev. Lett. , 176403(2008). The Kondo model, for example, exhibits interesting physics atexponentially long times, t (cid:29) (cid:126) / T K , where T K ∼ √ JN F e − / JN F is the exponentially small Kondo scale, with J the exchange con- stant and N F the density of states. J. K. Freericks, V. M. Turkowski, and V. Zlati´c, Phys. Rev. Lett. , 266408 (2006). H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner,ArXiv e-prints (2013), arXiv:1310.5329 [cond-mat.str-el]. E. Arrigoni, M. Knap, and W. von der Linden, Phys. Rev. Lett. , 086403 (2013). L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen,H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, andM. Wolf, Phys. Rev. Lett. , 067402 (2006). L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen,M. Wolf, H. Berger, S. Biermann, and A. Georges, New Jour-nal of Physics , 053019 (2008). P. A. Loukakos, M. Lisowski, G. Bihlmayer, S. Bl¨ugel, M. Wolf,and U. Bovensiepen, Phys. Rev. Lett. , 097401 (2007). J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby,M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard,Science , 2180 (2005). F. H. L. Koppens, C. Buizert, K. J. Tielrooij, I. T. Vink, K. C.Nowack, T. Meunier, L. P. Kouwenhoven, and L. M. K. Vander-sypen, Nature (London) , 766 (2006).
S. Loth, M. Etzkorn, C. P. Lutz, D. M. Eigler, and A. J. Heinrich,Science , 1628 (2010).
M. Greiner, O. Mandel, T. W. H¨ansch, and I. Bloch, Nature(London) , 51 (2002).
S. Will, T. Best, U. Schneider, L. Hackermueller, D.-S.Luehmann, and I. Bloch, Nature (London) , 197 (2010).
S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Scholl-woeck, J. Eisert, and I. Bloch, Nature Physics , 325 (2012). U. Schneider, L. Hackerm¨uller, J. P. Ronzheimer, S. Will,S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, andA. Rosch, Nature Physics , 213 (2012). A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher,A. Garg, and W. Zwerger, Rev. Mod. Phys. , 1 (1987). U. Weiss,
Quantum dissipative systems , Vol. 13 (World ScientificPub Co Inc, 2008).
K. L. Hur, Annals of Physics , 2208 (2008).
A. Weichselbaum, Phys. Rev. B , 245124 (2012). By exact we mean the equality stated in Eq. (13). The actual valueevaluated with NRG may deviate by a small amount from exactanalytic results, e.g. from Bethe ansatz results, however, this dif-ference has been shown to be negligible . L. Merker, A. Weichselbaum, and T. A. Costi, Phys. Rev. B ,075153 (2012). The factor e − iH Qp τ p . . . e − iH Q τ is the time evolution op-erator at time ˜ τ p = (cid:80) pi = τ i following p intermedi-ate quantum quenches described by H Q , . . . , H Q p . Hence, Q p + (cid:104) rem | e − iH Qp τ p . . . e − iH Q τ | se (cid:48) m (cid:105) i is the matrix element of thisoperator between the initial states of H i and the states of thequench Hamiltonian H Q p + . For times t > ˜ τ n . At shorter times, ˜ τ p + > t ≥ ˜ τ p , with p < n , oneonly needs to calculate the first p projected density matrices, sothe computational cost becomes less for shorter times. E. Iyoda and S. Ishihara, ArXiv e-prints (2013),arXiv:1312.1077 [cond-mat.str-el].
J. K. Freericks, H. R. Krishnamurthy, and T. Pruschke, Phys.Rev. Lett. , 136401 (2009).
M. Medvedyeva, A. Ho ff mann, and S. Kehrein, Phys. Rev. B ,094306 (2013). M. Eckstein and M. Kollar, Phys. Rev. B , 205119 (2008). P. Nordlander, M. Pustilnik, Y. Meir, N. S. Wingreen, and D. C.Langreth, Phys. Rev. Lett. , 808 (1999). M. Plihal, D. C. Langreth, and P. Nordlander, Phys. Rev. B , R13341 (2000).
F. Verstraete, V. Murg, and J. I. Cirac, Advances in Physics57