Reducing the numerical effort of finite-temperature density matrix renormalization group transport calculations

Finite-temperature transport properties of one-dimensional systems can be studied using the time dependent density matrix renormalization group via the introduction of auxiliary degrees of freedom which purify the thermal statistical operator. We demonstrate how the numerical effort of such calculations is reduced when the physical time evolution is augmented by an additional time evolution within the auxiliary Hilbert space. Specifically, we explore a variety of integrable and non-integrable, gapless and gapped models at temperatures ranging from T=infty down to T/bandwidth=0.05 and study both (i) linear response where (heat and charge) transport coefficients are determined by the current-current correlation function and (ii) non-equilibrium driven by arbitrary large temperature gradients. The modified DMRG algorithm removes an 'artificial' build-up of entanglement between the auxiliary and physical degrees of freedom. Thus, longer time scales can be reached.


Introduction
A physical system is usually characterized by its response to perturbations. In transport setups, one studies charge or energy currents driven by voltage or temperature gradients. From the theoretical point of view this is complicated -computing the quantum mechanical time evolution of a system in non-equilibrium is one of the most active areas of research in condensed matter physics. When the external perturbations are small, one can simplify the problem by resorting to the Kubo formalism. E.g., the charge conductivity σ(ω) describes the linear response current J induced by a small electric field (we will be more specific below): where the dynamical correlation function J(t)J is calculated in thermal equilibrium: with H being the Hamiltonian of the system and T denoting the temperature. Unfortunately, computing transport coefficients such as σ(ω) is generally still difficult: Even if one knows the exact thermal density matrix ρ T (or the exact ground state), extracting correlation functions, which couple all excitations, remains a formidable task.
The key question posed and addressed in this work is: How can linear-response and non-equilibrium transport properties of one-dimensional (1d) systems at finite temperature be calculated efficiently using the density matrix renormalization group (DMRG) [1,2]? DMRG was originally devised [3,4] as tool to accurately determine ground states of 1d Hamiltonians. The reason for its success became understandable when it was formulated using matrix product states (MPS) [5,6,7,8,9]: The area law [10] stipulates that the ground states of 1d systems governed by local Hamiltonians are only entangled locally; this implies that they can be expressed efficiently using a MPS with a small bond dimension χ (which encodes the amount of entanglement). The MPS describing a given system can be determined variationally -this is the very core of a ground state DMRG calculation [3,4]. One way to compute correlation functions A(t)B at T = 0 is to directly simulate the time evolution (for other approaches see Refs. [1,11,12,13,14,15,16]) using a time-dependent DMRG framework [17,18,19,20,21,22,23,24,25]. The corresponding algorithm can again be formulated elegantly using matrix product states. The physical growth of entanglement implies that the bond dimension needed to approximate A(t)B to a certain accuracy grows with time. This limits the accessible time scales. Standard DMRG methods allow computing the time evolution of a pure state and are thus not directly applicable at T > 0. Various approaches for simulating finite-temperature dynamics using DMRG have been explored within the literature [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]. These include probabilistic sampling over an appropriately chosen set of pure states [30], schemes which time-evolve operators instead of states [31,32], transfer-matrix DMRG [33,34,35,36], and exact representations through purification [37,38,39]. Purification expresses the thermal statistical operator ρ T as a partial trace over a pure state |Ψ T living in an enlarged Hilbert space where auxiliary degrees of freedom Q encode the thermal bath [41]: One of the main advantages of this approach is that all the standard methods for time evolving quantum states within DMRG are directly applicable. In this work, we employ the time evolving block decimation [17,22]. One of the first applications of finite-temperature 'purification DMRG' to dynamical problems was the calculation of spin-spin correlation functions of integrable spin-1/2 Heisenberg chains [38]. While DMRG yields data which are 'numerically exact' (this is verified by comparing with analytic results available for 'non-interacting' models [33,34,38]), the time scales accessible at finite temperatures are considerably smaller than at T = 0. This observation, which might be one reason why studies of finite-T dynamics [31,32,33,34,35,36,38,39,40,42,43,44] are rarer than their T = 0 counterparts, can be understood as follows. Assume we want to compute a ground state correlation function, i.e. evaluate Eq. (3). Under the time evolution, the entanglement grows locally around the region on which B acted. In contrast, at T > 0 where one needs to calculate (see below for details) the entanglement generally increases homogeneously throughout the whole system. This holds true even for B = 1, i.e. when the state |Ψ T is exposed to a supposedly trivial time evolution [40]. The reason for this is that purification is not unique and various representations of the same ρ T differ in their degree of entanglement between the auxiliary and physical degrees of freedom. Even when one starts out with a purification that minimizes this entanglement, it rapidly grows under the DMRG time evolution. This observation naturally leads to the question: Can that very same 'freedom' be used to our advantage to undo the non-physical growth of entanglement? In other words, is there a unitary transformation U Q (t) acting on the auxiliary Hilbert space Q which removes the artificial entanglement? In Ref. [40] we proposed to time-evolve Q backwards in time with the physical Hamiltonian acting on the auxiliary degrees of freedom: This renders the time evolution of |Ψ T trivial, and the evaluation of Eq. (5) is therefore eventually only plagued by an entanglement building up around the region where B acts in complete analogy to the ground state calculation (the physical reason being quasilocality [31,32]). This generically leads to a slower increase of the bond dimension χ, and thus longer time scales can be reached. In Ref. [40], the potential of the modified DMRG algorithm was demonstrated for the spin-spin correlation function S z n (t)S z m of the XXZ chain at zero anisotropy ∆ = 0, which maps to free fermions and thus allows for an exact (benchmark) solution. A more thorough comparison between the numerical effort of the standard [38] and modified [40] algorithms in calculating S z n (t)S z m for arbitrary ∆ can be found in Ref. [32], where further optimizations using operator-space DMRG were explored. One of the ideas of Ref. [32] is to use time translation invariance to rewrite Eq. (2) as which allows accessing times twice as large without additional effort. Refs. [31,32] also show that for certain scenarios it can be more efficient to redistribute the evaluation of exp(−iHt) and exp(−βH) over two DMRG simulations. From the point of view of physical applications, the modified DMRG schemes were used to investigate current correlation functions and the Drude weight of the integrable XXZ chain at intermediate to high temperatures [40,44], scaling properties of the DC conductivity in presence of non-integrable perturbations [42], non-equilibrium induced by temperature gradients [43], and spectral functions of hardcore bosons [31]. For reasons of completeness, we mention other approaches to linear response transport properties of the XXZ model (and related ones). These include exact Bethe ansatz calculations [47,48,49,50], integrability arguments [51,52,53,54,55,56], field theories [35,36,57,58,59], quantum Monte Carlo [60,61,62], exact diagonalization [63,64,65,66,67], transfer matrix DMRG [35,36], and dynamical DMRG [14]. For a non-exhaustive list of prior works on non-equilibrium thermal transport see Refs. [23,68,69,70,71,72,73,74,75,76,77,78].
The first and foremost goal of this paper is to present extensive quantitative data on how the build-up of entanglement is reduced if the auxiliary Hilbert space Q is time-evolved with the physical Hamiltonian but reversed time. We focus on transport properties in linear response [40,42,44] and in thermal non-equilibrium [43]. Specifically, we study (i) homogeneous gapless and gapped spin-1/2 XXZ chains, also in presence of various perturbations which break integrability, (ii) the quantum Ising model, and (iii) impurity setups of a quantum dot connected to non-interacting leads. For temperatures from T = ∞ down to T /bandwidth = 0.05, the modified DMRG algorithm leads to a slower increase of the MPS dimension χ. Only at low T does the standard approach U Q = 1 become more efficient (see Ref. [32] for details). For most (but not all; we will be more specific below) problems studied in this work, T /bandwidth = 0.05 is a low enough temperature to correspond to the T = 0 limit, which one can establish, e.g., by comparing with field theory results [77]. We show the typical behavior of χ on the relevant physical time scales for each problem at hand (e.g., on the time scale at which non-equilibrium currents generically reach their steady state). We reiterate how Eq. (7) can be implemented within the purification approach and illustrate (following Refs. [31,32]) how it allows accessing times twice as large.
As the second purpose of this paper, we present technical details of the implementation of the algorithm. E.g., we show how to time-evolve next-nearest neighbors (closely following Ref. [2]), which is necessary because physical degrees of freedom are separated by an auxiliary site due to the purification. We discuss the numerical accuracy of our data, compare with exact results, investigate to what extend the choice of U Q (t) = exp(iHt) is optimal [31,32], and provide further evidence for the reliability of linear prediction extrapolation schemes [38,42,45,46] for the spin-spin correlation function of the isotropic Heisenberg chain.

Models
As a prototypical model, we consider a chain of interacting spin-1/2 degrees of freedom S x,y,z n governed by local Hamiltonians or equivalently spinless fermions through a Jordan-Wigner transformation. By choosing the couplings J n , ∆ n , and b n appropriately: we can study systems which are gapless or gapped, and investigate the role of integrability. For λ = 1 and b = 0, Eq. (8) can be diagonalized via Bethe ansatz [79]; the model is non-integrable otherwise. The energy spectrum is gapless for |∆| ≤ 1 and gapped for ∆ > 1. A gap opens for λ < λ [81,80,42]. In addition, we study the quantum Ising model  where the corresponding current operators J = n j n are defined through a continuity equation [82]: One usually decomposes the real part of σ(ω) as where the so-called Drude weight D is the prefactor of the singular contribution. D can be determined from the long-time asymptote of the current-current correlation function As already stated above, a time-dependent equilibrium correlation function is defined as with Z T = Tr exp(−H/T ) denoting the partition function.

Non-equilibrium.
In addition to linear response, we study two thermal nonequilibrium setups. The first [labeled 'non-equilibrium A' and depicted in Figure 1(A)] is introduced via the following protocol: We initially consider two separate chains, each being in thermal (grand-canonical) equilibrium at temperatures T L and T R . The corresponding density matrix factorizes, At time t = 0, the chains are coupled through h L/2 , and the time evolution of any observable A is computed using H = H 0 + h L/2 : In the second setup ['non-equilibrium B'; see Figure 1(B)] we investigate two noninteracting chains ∆ = b = 0, λ = 1 of length L/2, at different temperatures T L and T R . At t = 0, they are coupled via an interacting resonant level model: The site n = L/2 + 1 is initially in an equal superposition of up and down states (i.e., formally at infinite temperature).

DMRG
In this Section, we give a brief overview of the DMRG method [1,2,3,4]. More details can be found in the Appendix. In order to evaluate Eqs. (15) or (18) by a standard DMRG algorithm (which time evolves wave functions) one first needs to purify the thermal density matrix ρ T by introducing an auxiliary Hilbert space Q such that ρ T = Tr Q |Ψ T Ψ T |. This is analytically possible only at T = ∞ where ρ T factorizes. However, |Ψ T can be obtained from |Ψ ∞ by applying an imaginary time evolution, [26,37,38]. The correlation function of Eq. (15) is exactly recast as and this object is directly accessible in standard time-dependent DMRG frameworks [17,18,19,20,21,22]. It is convenient to first express the initial state |Ψ ∞ in terms of a matrix product state [5,6,7,8], where Here and in the following σ i is a short hand for either a physical or auxiliary degrees of freedom: σ i ∈ {σ n , σ n Q }. The initial matrices associated with |Ψ ∞ read as well as Λ i = 1. After factorizing the evolution operators exp(−λH) using a second or fourth order Trotter decomposition, they can be successively applied to Eq. (23). At each time step ∆λ, two singular value decompositions are carried out to update three consecutive matrices. The matrix dimension χ is dynamically increased such that at each time step the sum of all squared discarded singular values is kept below a threshold value ǫ.
An exact modification to the finite-temperature DMRG algorithm (which allows accessing longer time scales) was recently introduced in Ref. [40]: One has the analytic freedom to apply any time-dependent unitary transformation to the in principle inert auxiliary sites σ n Q -physical quantities are determined by the trace over Q and are thus not affected by U Q (t): For example, Eq. (21) can be rewritten as where Put differently, purification is not unique. It turned out [31,32,40,42,43] that choosing i.e., time-evolving the auxiliary sites with the physical Hamiltonian but reversed time, leads to a slower increase of χ and thus longer time scales can be reached. An intuitive way of understanding this will be given below. Only at low temperatures does using U Q = 1 become more efficient [31,32]. In practice it is most convenient to implement

Results
In this section we present our main results. We first quantify how the entanglement growth is reduced by time-evolving the auxiliaries with the physical Hamiltonian but reversed time. In an extension of earlier studies (see Refs. [31,32,40] and the discussion in Sec. 1), we study gapless and gapped XXZ chains in presence of non-integrable perturbations and focus on transport problems both in linear response and out of equilibrium. We investigate to what extend U aux (t) = exp(itH) is optimal within our approach. In passing, we provide another example of the stability of linear prediction (see also Refs. [2,38,42,45,46]) and recast (one of) the ideas of Refs. [31,32] in the language of purification.  Figure 2 shows the evolution of the bond dimension χ during the calculation of e −iHt U Q (t)|Ψ T . For the standard approach U Q (t) = 1, χ increases with time: the state e −iHt |Ψ T which purifies the density matrix becomes successively more entangled.

It is instructive to first consider the trivial case of
Choosing U Q (t) = e iHt removes this artifact. This immediately indicates why longer time scales can be reached in DMRG evaluations of, e.g., current correlation functions e −iHt U Q (t)j n |Ψ T . The entanglement only builds up locally around site n (the physical reason being quasi-locality [31,32]), and likewise for our non-equilibrium setup it grows locally around the interface region over which the initial temperature gradient falls off.
The artificial global build-up of entanglement is removed if the auxiliaries are evolved in time. Only at low T choosing U Q = 1 becomes more efficient [31,32] since the time evolution of the ground state is trivial (the latter is indicated in Figure 2: χ grows more slowly at low T = 0.2 than it does at high T ). An intuitive understanding of why the particular choice of U Q (t) = e +iHt removes the 'artifical' entanglement growth was recently provided in Refs. [31,32]: In an operator-space formulation, the modified DMRG algorithm corresponds to a Heisenberg time evolution of the matrix product operators representing A and B in Eq. (27). If A (and likewise for B) is local, most terms in the first Trotter step e iH∆t Ae −iH∆t cancel out. As time evolves, . . . e iH∆t e iH∆t Ae −iH∆t e −iH∆t . . . , entanglement can only build up gradually around the region on which A acted due to the so-called light-cone effect in nonrelativistic systems: Lieb-Robinson bounds [85] generically stipulate that correlation functions A(t)B of local operators A and B fall featuring an initial sharp temperature gradient T L = T R . The system is gapless for zanisotropies |∆| ≤ 1 and gapped otherwise; it becomes non-integrable in presence of a finite dimerization λ < 1. Note that for ∆ = 0.5 and λ = 1 asymptotic low-temperature behavior described by a field theory [77] sets in for T 0.2. (a) Energy current at the interface. Choosing a discarded weight control parameter around ǫ ∼ 10 −6 at a Trotter stepsize of ∆t = 0.2 is sufficient to accurately obtain its steady-state value. Inset: The actual discarded weights during a whole Trotter step and during one substep for ∆χ = 20. The difference to ǫ is explained in the main text. (b) Evolution of the dimension of the corresponding matrix product state. If the auxiliary degrees of freedom which purify the thermal density matrix are time-evolved with the physical Hamiltonian but reversed time (which is an exact modification to the standard algorithm), the build-up of entanglement is reduced. The computational cost of a singular value decomposition (which dominates the DMRG algorithm) scales as χ 3 .
off exponentially for x > vt, with x denoting the spatial distance between the regions on which A and B act, and v being some velocity with which excitations propagate.

Reduction of the growth of entanglement: Non-equilibrium
In this section we illustrate quantitatively the effects of time-evolving the auxiliaries for the non-equilibrium setups sketched in Figure 1. We start by studying an interacting XXZ chain (∆ = 0) with two additional perturbations (dimerization λ < 1 and a staggered field b > 0) that both render the system non-integrable [see Eq. (8) for a precise definition]. At time t = 0, two 'semi-infinite' chains (typical lengths being L/2 ∼ 100 − 200) prepared in thermal equilibrium at temperatures T L,R are coupled by h L/2 to an overall 'translationally-invariant' chain. The temperature gradient drives an energy current j E,n (t) whose typical behavior is shown in Figure 3(a). For b = 0 the current at the interface n = L/2 saturates to a finite value on a scale t ∼ 1 − 10 irrespective of whether the system is gapless or gapped (indicating that either the nonintegrable dimerized chain supports dissipationless transport or that its current decays on a hidden large scale; see Ref. [43] for further details). The steady-state current of the XXZ chain is of the simple 'black-body' form [43,76,77,86] lim implying that it is determined solely by the linear conductance G(T ) ∼ ∂ T f (T ) for arbitrary large T L − T R : This agrees with a field theory result [77] at low temperatures; asymptotic low-T behavior sets in around T 0.2. The time evolution of the corresponding bond dimension χ is shown in Figure  3(b) for the integrable XXZ chain both for parameters where it is gapless (∆ = 0.5) and gapped (∆ = 2). Figure 4(a) shows the same in presence of non-integrable perturbations as well as for the quantum Ising model of Eq. (10) at criticality g = 1. In all cases and for all temperatures from T L,R = ∞ down to T L,R = 0.1 (which for this problem corresponds to zero temperature), time evolution of the auxiliaries leads to a slower increase of χ; note that the computational cost of a singular value decomposition (which dominates the whole DMRG algorithm) scales as χ 3 .
The same reduction of χ holds for the non-equilibrium impurity problem depicted in Figure 1(b) where an interacting resonant level model defined in Eq. (20) is coupled to two non-interacting leads (∆ = b = 0, λ = 1) at different temperatures. An exemplary evolution of the MPS dimension is shown in Figure 4(b); a discussion of the (transport) physics at small interactions or zero temperature can be found in Refs. [87]. Note that due to the appearance of a Kondo-like scale, T /bandwidth = 0.05 does not necessarily correspond to zero temperature for this problem [87]. The discarded weight ǫ = 10 −6 in Figure 3(b) is chosen small enough to 'accurately' determine the current and in particular its steady-state value. This is illustrated in Figure 3(a) where ǫ is successively lowered from ǫ = 10 −4 to ǫ = 10 −7 . Moreover, the non-interacting case ∆ = b = 0, λ = 1 allows for an exact evaluation of the steadystate current [76], which can be used to benchmark the accuracy of the DMRG data (a comparison can be found in Ref. [43]). The inset to Figure 3(a) shows a generic example for the actual discarded weight (for the precise definition of ǫ see Appendix A.2.4). It is always bound by 2ǫ; the total discarded weight during a complete Trotter step e −iH∆t e iH∆t is typically bound by 10ǫ. It is also instructive to carry out a calculation at a fixed MPS dimension χ rather than a fixed discarded weight. The result is shown in Figure 3(a).

Reduction of the growth of entanglement: Equilibrium
We now turn to current correlation functions in thermal equilibrium. Following Eqs. (12) and (A.8), we need to evaluate The AC conductivity σ(ω) is determined by the Fourier transform of J(t)J via Eq. (11). The integrable gapless XXZ chain, on which we focus in this section, supports dissipationless transport at finite temperature, i.e. its Drude weight D in Eq. (13) is finite. D can be obtained from the long-time limit of the current-current correlation function via Eq. (14). The relevant time scale in this problem is thus the scale on which J(t)J saturate to their asymptote. DMRG data for the charge current correlation function at ∆ = 0.71 and for various discarded weights is shown in Figure 5 n e ikn S z n+L/2 (t)S z L/2 and J C (t)J C of the XXZ chain at ∆ = 0. DMRG results are compared with the exact solution obtained by mapping the model to free fermions. Since J C is conserved, J C (t)J C is constant (we dropped some DMRG data points to improve visibility). However, every single term j C,n (t)j C,m that contributes to the sum is time-dependent. (b) Demonstration of the stability of linear prediction (see also Ref. [38] and Refs. [2,31,42,45,46]) for the isotropic chain ∆ = 1. Dashed lines were obtained by fitting the DMRG data for times t fit ∈ [2. 5,5] and subsequent extrapolation to t > 5 using linear prediction. The extrapolated data almost coincides with the DMRG data for t > 5 (solid lines).
T 0.5, we can access time scales on which J(t)J saturates [we stop our simulation once the MPS dimension has reached values around χ ∼ 1500; see the numbers given in Figure 5(a)]. More details can be found in Refs. [40,44]. A fairly large ǫ ∼ 10 −5 is sufficient to obtain D quantitatively. This is supported by the comparison with the exact solution for ∆ = 0 shown in the inset to Figure 6(a). Note that J C (t)J C is constant for ∆ = 0 since J C is conserved; however, every single term j C,n (t)j C,m that contributes to the sum is time-dependent. It is again instructive to carry out a calculation using a constant χ rather than a fixed discarded weight; the result is shown in the upper panel of Figure 5(a).
For all temperatures from T = ∞ down to T = 0.125 time evolution of the auxiliaries leads to a slower increase of the dimension of the matrix product state. This is depicted in Figure 5(b). One might expect low-T behavior of D(T ) to set in for T = 0.125 (see Refs. [35,36,44]), but the current correlation function decays so slowly that we cannot access its asymptotics without extrapolation. Exemplary data for the gapped phase is shown in the lower inset to Figure 5(b).

Spin correlation function
We briefly present data for the spin-spin correlation function S zz (n, t) = S z n+L/2 (t)S z L/2 , S zz (k, t) = n e ikn S zz (n, t) , which provides a standard testing ground for DMRG approaches to time-dependent correlation functions at finite temperature [33,34,38,40]. The XX chain ∆ = 0 again allows for an exact solution by mapping the model to free fermions. Figure 6(a) shows a comparison of DMRG data with the exact result for T = 1 and T = 0.1 (which for this situation corresponds to low temperature). A large discarded weight ǫ ∼ 10 −5 is sufficient to reproduce the analytic result up to times t ∼ 30. If the auxiliaries are time-evolved, the MPS dimension increases only moderately to χ ∼ 100 − 200 in contrast to the standard approach (see Ref. [40] and also Refs. [33,34] for transfermatrix DMRG data). For finite ∆, the choice of U Q (t) = e iHt still outperforms U Q = 1 (expect at low temperatures), but the difference becomes successively less pronounced [32]. For the isotropic chain ∆ = 1 and all T 0.1, the accessible time scales are about a factor 1.5 − 2 larger if the auxiliaries are evolved in time. In Ref. [38] (see also Refs. [2,31,42,45,46]), linear prediction was first established as an accurate tool to extrapolate correlation functions by considering exactly-solvable models and subsequently used to compute the Fourier transform of S zz (k, t) of the isotropic chain. For reasons of completeness, we revisit the isotropic chain and compare the result of linear prediction with DMRG data for larger times. Figure 6(b) illustrates that the agreement is good.

Optimization of U Q
Finally, we shortly investigate to what extend our choice of e iHt for the time evolution of the auxiliaries is optimal (more details on the optimization of U Q can be found in Refs. [31,32]). To this end, we consider with N begin an arbitrary Hermitian matrix. We compute the entanglement entropy (b) By rewriting S zz (n = 0, t) = S z L/2 (t/2)S z L/2 (−t/2) and time-evolving S z L/2 (±t/2) independently the time scale t can be accessed by a MPS of a smaller dimension χ (this idea was introduced in Refs. [31,32]).
between the left and right halves of the time-evolved state e −iHt U η Q (t)|Ψ T which purifies the density matrix. We have verified that for arbitrary N coupling nearest neighbors, S ent (t, η) is at least quadratic in η. Figure 7 illustrates this for the two choices as well as As mentioned above, Refs. [31,32] present thorough details on how to further optimize finite-temperature dynamical DMRG. We here reiterate one of the main ideas using the language of purification. In thermal equilibrium, one can recast any correlation function exploiting time translation invariance: Thus, one only needs to carry out time evolutions up to times |t/2| in order to compute A(t)B . If A = A † = B (e.g., for current-current correlation functions), it is sufficient to perform a single calculation This implies that one can access a time scale twice as large without much additional effort. This is illustrated in Figure 8(b) for the spin-autocorrelation function of the isotropic Heisenberg chain.

Summary
We presented extensive quantitative data for how the growth of entanglement in finite-temperature dynamical DMRG calculations can be reduced by time-evolving the auxiliary degrees of freedom which purify the thermal statistical operator. The time evolution of the auxiliaries is an exact modification to the standard algorithm. Our particular focus was on energy and charge transport properties both in linear response (described by current-current correlation functions) and in non-equilibrium driven by temperature gradients. We studied a variety of gapless and gapped integrable and non-integrable spin-1/2 chains (i.e., interacting spinless fermions), the quantum Ising model at criticality, and impurity (quantum dot) setups. For all temperatures from T = ∞ down to T /bandwidth = 0.05, which for most problems investigated in this work corresponds to low temperatures, the build-up of entanglement is reduced. This speeds up numerics and allows to eventually access longer time scales.
Acknowledgments -We are grateful to Frank Verstraete for very useful suggestions and thank Thomas Barthel, Fabian Heidrich-Meisner, and Steve White for discussions. Support by the Deutsche Forschungsgemeinschaft via KA3360-1/1 (CK), by the DARPA TI program of UCLA (JHB) as well as by the Nanostructured Thermoelectrics program of LBNL (CK and JEM) is acknowledged.
Appendix A. Technical details.
One way to evaluate Eqs. (15) and (18) within the density matrix renormalization group is to purify [41] the thermal density matrix by introducing an auxiliary Hilbert space Q: At infinite temperature, where ρ T is proportional to the unit operator, one can analytically write down the purification: In order to exploit conservation laws within the DMRG algorithm, it is convenient to choose where σ n =↑, ↓ denotes the eigenbasis of S z n , and Q is spanned by One readily verifies that indeed The finite-temperature purified state |Ψ T is obtained from |Ψ ∞ by applying an imaginary time evolution [2], where we formally replaced The correlation function of Eq. (15) is exactly recast as and for two time arguments one similarly finds In order to evaluate Eq. (7) for A = A † = B, it is sufficient to compute The expectation value of any observable A with respect to the time-dependent (nonequilibrium) density matrix of Eq. (18) is given by Appendix A.2. DMRG algorithm Appendix A.2.1. Initial state. It is convenient to first express the initial state |Ψ ∞ in terms of a matrix product state [5,6,7,8], Here and in the following σ i is a short hand for either a physical or auxiliary degrees of freedom: σ i ∈ {σ n , σ n Q }. The initial matrices corresponding to Eq. (A.3) read as well as Λ i = 1. Normalization generally stipulates as well as In order to reduce the error, one can employ a fourth order decomposition [84], and We typically employ a second order Trotter decomposition with a time step of ∆β = 0.005 during the imaginary time evolution from β = 0 down to the physical temperature β = 1/T and a fourth order decomposition with a time step of ∆t = 0.2 during the real time evolution. This is sufficient for all problems studied in this paper (which we verify by comparing against data obtained for smaller ∆β = 0.0025 and ∆t = 0.1). The Trotter decomposition is typically not a significant source of error within time-dependent DMRG calculations [68].
Appendix A.2.3. DMRG step. The physical Hamiltonians of Eqs. (8) and (10) couple matrices with indices σ n and σ n+1 which in the matrix product state of Eq. (A.12) are next-nearest neighbors -they are separated by a matrix associated with an auxiliary degree of freedom σ n Q . Through Eq. (A.7), the local Hamiltonians of Eqs. (8) and (10) are formally replaced by Thus, after each application of exp(−∆λh n ), two singular value decompositions are carried out to update three consecutive matrices; to keep the notation transparent, let us denote those matrices as: and furthermore abbreviate h = h n . This 'DMRG update' can be achieved by a simple and straightforward generalization of the protocol for nearest neighbors outlined in Ref. [2]. First, one forms the three-site tensor, which is then acted on by exp(−∆λh): Now, two singular value decompositions (SVDs) are applied to appropriately reshaped tensors in order to obtain the new matricesλ i andΓ σ i :  16) are preserved automatically (up to errors associated with a potential truncation of the matrices; see below). For a non-unitary exp(−∆βh), however, the matrix product state needs to be brought back to its canonical form in order to simplify the evaluation of expectation values, and, more importantly, because otherwise the truncation below is not optimal. Normalization can be achieved approximately by successively acting with unit operators ∆λ = 0 [22] or exactly using the procedure outlined in Ref.
[83]. χ (a more thorough discussion can be found in Ref. [2]). The computational effort of the DMRG algorithm is dominated by the cost of the singular value decompositions, which scales as χ 3 . Increasing χ → 2χ during each sub-step of a Trotter step is numerically unfeasible (and fortunately generically unnecessary). Thus, we truncate the matrices to a dimension χ ′ < 2χ by dropping the indices a i which correspond to the smallest singular valuesΛ a i . For a normalized state [i.e., if Eqs. (A.15) and (A. 16) hold], this is the best approximation with respect to the Hilbert space norm. The error (i.e., the norm-distance of the states before and after truncation) is given by the sum of all squared singular values which are discarded during all sub-steps of one Trotter step. This socalled discarded weight is the key numerical control parameter; the DMRG algorithm becomes exact as it is sent to zero. One can strictly enforce a fixed discarded weight by analyzing the singular values during all SVDs carried out within one Trotter substep and then retrospectively choosing χ ′ appropriately (or by simply repeating the substep if χ ′ > χ). More pragmatically, one can fix χ ′ = χ and increase χ by ∆χ after one substep if the total discarded weight exceeds a pre-defined threshold ǫ. While this simple approach certainly becomes inapplicable if one aims at strictly maintaining a discarded weight that is very small, it works for all problems studied in this paper, with ∆χ ∼ 10 − 20 being a reasonable choice. We observe that the actual discarded weight is smaller than 2ǫ for all data presented in this paper [it is imperative to always monitor the real discarded weight which determines the actual error; this is illustrated in the inset to Figure 3(a)]. In the following we neglect this difference and simply refer to our control parameter ǫ as 'discarded weight'. We successively lower ǫ from 10 −8 − 10 −13 during the imaginary time evolution and from 10 −3 − 10 −7 during the real time evolution, which turns out to be sufficient to obtain 'accurate' results for all physical quantities of interest. 'Accurate' means 'sufficiently converged with respect to lowering ǫ', or 'in reasonable agreement with exact data' in case that the problem allows for an analytic solution; all this is