Real time evolution at finite temperatures with operator space matrix product states

We propose a method to simulate the real time evolution of one dimensional quantum many-body systems at finite temperature by expressing both the density matrices and the observables as matrix product states. This allows the calculation of expectation values and correlation functions as scalar products in operator space. The simulations of density matrices in inverse temperature and the local operators in the Heisenberg picture are independent and result in a grid of expectation values for all intermediate temperatures and times. Simulations can be performed using real arithmetics with only polynomial growth of computational resources in inverse temperature and time for integrable systems. The method is illustrated for the XXZ model and the single impurity Anderson model.


Introduction
The simulation of strongly interacting quantum many-body systems is, in general, a computationally difficult task. Even in the case of one-dimensional quantum systems with finite range interactions, where ground state properties have become readily accessible using the density matrix renormalization group (DMRG) algorithm [1], accessing the dynamics remains hard. While, the matrix product states (MPS) [2,3,4] that are at the core of DMRG manage to describe the ground states due to their weak entanglement [5,6], the entanglement increases in the real time evolution [3,5,7,8] making the study of dynamics computationally challenging even for pure states.
Even harder is dynamics at finite temperatures. MPS based methods operate on the level of wave functions. Systems at thermal equilibrium, described by a mixture of quantum states, thus need to be purified by being described as a pure state in an enlarged system, consisting of the system and an auxiliary environment. The simulation is then performed on this purified system [4,9] and by tracing over the degrees of freedom of the environment, the properties of the original system are obtained. This idea was first implemented in Ref. [4,10] for density matrices of thermal quantum systems and systems far from equilibrium. The same concept was used in the simulation of dynamical properties of thermal systems where the real time evolution was computed on the level of the purified system and the density operator for the original system was obtained by tracing over the environment [11,12].
A disadvantage of these purification based approaches is that they require additional information on the environment, and especially on its dynamics. In general, these methods are accompanied by an exponential growth of computational resources with time and elaborate procedures are required to control the fast growth of entanglement. Recent advances in the application to the XXZ model have shown that a smart choice of the environmental dynamics can significantly extend the accessible simulation times [13]. However, the question of the optimal environment and its dynamics is still open. Even finding a good choice of the environmental dynamics requires deep knowledge of the system and is far from being straightforward.
Here we propose an alternative strategy, based on operator space concepts and matrix product states. Similarly to the purification approach, we describe the density matrices as many-body states; however, our states are not superpositions of quantum states but superpositions of operators [14,10]. In the same way we also describe the observables and their dynamics, known as the time dependent DMRG in the Heisenberg picture [14,15]. An operator space approach was used in [16], however, still retaining the exponential growth of computational costs. The key feature of our approach is the combination of both concepts where the expectation value of an observable at a given time and temperature is calculated as a scalar product of two vectors in the operator space, one representing the density matrix and the other the real time evolution of an operator. This way, the simulation of density matrices and the real time evolution of operators are decoupled and can be performed independently such that only two simulations are sufficient to calculate a grid of expectation values at various intermediate (inverse) temperatures and times.
The computational costs to simulate the density matrices grow polynomially in time due to a low degree of correlation at thermal equilibrium [17], although such description has mostly been used to simulate systems far from the equilibrium [18]. Similarly, the simulation of local operators also exhibits polynomial growth of resources [14] for integrable systems. This corresponds to a logarithmic growth of the operator space equivalent of the entanglement entropy (OSEE) [19], which is reminiscent of the situation in a pure-state local quantum quench [20,21]. Our method can therefore be used to simulate the real time evolution of local operators in integrable quantum systems in polynomial time.

Methods
Let us consider a quantum system of n sites with a Hilbert space H. The expectation value of an operator a : H → H in a mixed state described by a density matrix ρ, is given by a = tr(ρa) with tr(ρ) = 1. Like the operator a, the density matrix is also a linear map over H and these maps form the operator space algebra K with an inner product (a|b) = (dimH) −1 tr(a † b). The expectation value of a can thus be given as where the bracket notation |a) is used when referring to a ∈ K and |e) corresponds to the fully mixed state (the identity). In thermal equilibrium for inverse temperature β = (k B T ) −1 , the expectation values are given by using |ρ(β)) = e −βH as the density matrix whereas the time evolution of Eq. (1) is obtained by replacing a by the Heisenberg operator a(t) = e itH ae −itH . Clearly, the simulations for ρ(β) and a(t) can be done independently.
The density matrices ρ(β) are obtained by the imaginary time evolution, starting from the infinite temperature equilibrium state |e). By defining a map over the operator space,χ : K → K, which left-multiplies an operator by H, i.e.χ|a) ≡ |Ha), the density matrix |ρ(β)) = e −βH becomes |ρ(β)) = e −βχ |e). ( Similarly, the Heisenberg evolution a(t) = e itH ae −itH can be simulated by introducing a linear map defined asĤ|x) ≡ |[x, H]) for x ∈ K. It was shown [19] that the Heisenberg evolution is hence transformed to the Schrödinger evolution in the operator space Bothχ andĤ have the same form and the same range of interactions as H, if the basis for K is local (see Appendix A).
Let us now define a suitable basis for K. While any basis could in principle be used, it is advantageous to choose a basis given by direct products of local basis operators {p [j] ν j } pertaining to the site j, as P ν = p νn . These operators are chosen such that they form an orthonormal set. A suitable choice for {p [j] ν j } depends on the type of simulation and a different choice can be made for the thermal (2) and the real time evolution (3). For the former, the mapχ, represented by [χ] µ,ν = (dimH) −1 tr(P † µχ P ν ), is real if {P ν } are real which in turn allows us to simulate the thermal density matrices (2) with real arithmetics.
For the real time evolution ν,µ andĤ is represented by a hermitian skew-symmetric matrix. Together with the imaginary unit in the exponent of (3), the Heisenberg evolution is generated by a real valued map. Finally, due to different operator space bases used in the thermal and the Heisenberg evolution, the calculation of expectation values must be preceded by a basis transformation, which in our case is a product operator and has no significant effects on the computational costs.
For reference we list a few typical operator space bases. In the case of spin-1/2 models, spinless fermions or hard core bosons the basis is given by the products of Pauli matrices p [j] µ j = σ µ j j for µ j ∈ {0, 1, 2, 3}. This basis is hermitian and thus suitable for the real time evolution whereas the thermal evolution is done using a real local basis {σ 0 , σ 1 , −iσ 2 , σ 3 }. The local basis for spin-1 systems is provided by Gell-Mann matrices whereas for spin-3/2 it can be chosen as {S (a,b) = σ a ⊗ σ b }. In all cases, the local basis is hermitian and its real representation can be found in a straightforward way.
Let us now focus on the MPS ansatz for the operators. For simplicity we shall assume a spin-1/2 model although the results are readily generalized to any many-body quantum system with a local dimension d. The operators P ν 1 ,...,νn = σ ν 1 1 · · · σ νn n form the basis set |ν) ≡ |P ν ) and a ∈ K is given in terms of a MPS [14] |a) = ν 1 ,...,νn (2) and (3) are simulated using the Suzuki-Trotter decomposition and local updates of the MPS, known as the time dependent DMRG algorithm [5,6,7,8] which is particularly efficient when H contains at most nearestneighbor interactions.
The thermal expectation value A(t) β given by (1) is obtained by combining the independent simulations (2) and (3). Using the same principles, we can simulate a(t)b β where a(t) and b are MPS, by constructing a matrix product operatorB : x → bx. The result is given by b a(t) β = (ρ(β)|B|a(t))/(ρ(β)|e) which in the MPS language corresponds to calculating an expectation value. Similarly a(t) b β is obtained.

Results
The method proposed here can be used to simulate the real time evolution of a wide range of models at finite temperature, for local operators such as the spin projection at a given site, a local current, or correlations of local operators (e.g. Green's functions). We start by analyzing the complexity of both constituents of the method, the evolution of the density matrices and the time evolution of operators in the Heisenberg picture. We quantify the complexity in terms of the OSEE [19] which is an operator space analogue of the entanglement entropy for quantum states, defined as S (x) = −tr[R log 2R ] for R = tr E (|x)(x|). The bipartition is chosen symmetrically and thus the environment E in the partial trace includes half of the system. The effective bond dimension then scales as exp(S (x)) with the OSEE. For density matrices, it was also observed [17] that the OSEE is related to the quantum mutual information which quantifies the total correlations between two parts of the system.
We first consider a spin-1/2 XXZ model with open boundary conditions described by the Hamiltonian where ∆ is the anisotropy parameter. The OSEE of the initial state ρ(β = 0) is zero (the fully mixed state is separable). As shown in Fig. 1 (left), we observe logarithmic growth of the OSEE in the inverse temperature [17]. Furthermore, we simulate the , with σ ± m = (σ x m ± iσ y m )/2 for m = n/2. The results in Fig. 1 (right) show the logarithmic growth of OSEE which is an extension of the results in [14]. Similar to the evolution of density matrices, we observe an increase of OSEE for larger anisotropies ∆ where the simulation is more expensive. The logarithmic growth of OSEE corresponds to polynomial growth of the bond dimension [14,19] and the simulation can be performed efficiently. We note that the computational cost is expected to grow exponentially [14] for non-integrable systems, like in other similar methods.
The simulation of the density matrices and the local operators can now be used to calculate e.g. the spin current correlation function jj n/2 (t) β . The results are shown in Fig. 2, with an inverse temperature β = 0.25. The infinite-time limit of these currentcorrelations yields the Drude weight which can be used to signal ballistic transport in various integrable models, see [22]. For the XXZ chain similar results were obtained in Refs. [23,13,24], a detailed discussion is however beyond the scope of this manuscript and will be presented elsewhere. When compared to Ref. [13] (with our time scale multiplied by a factor of 4 due to the spin operators) we observe, on one hand, that the reachable times are approximately twice as large in the gapped regime (∆ = 2 in Fig. 2) with a maximal bond dimension D = 1000. On the other hand, the computation in the gapless regime, ∆ < 1, is less demanding and the curves obtained with a smaller bond dimension D = 800 are essentially indistinguishable in Fig. 2. This is in complete agreement with the picture obtained from the OSEE scaling. In turn, for small ∆ the simulation can be pushed beyond the time scale shown in Fig. 2, similarly to Ref. [13].
As another illustration of the method we calculate the Green's functions of the  single impurity Anderson model which can be mapped [25,26] to an XXZ-like model with n j = σ − j σ + j and the interaction strength U . The hopping parameters τ −1 = τ 1 are given by the hybridization whereas the rest of τ j are determined from the discretization of the bath (and τ 0 = 0), see [26] for details. The interaction term is only present between sites 0 and 1 (corresponding to spin-↑ and spin-↓ of the impurity) which allows longer acessible times and lower temperatures. We simulate the single particle Green's function G(t) = −i {f † ↑ , f ↑ (t)} β where f ↑ is the annihilation operator at the impurity for spin-↑, written as f ↑ = ( j<0 σ z j )σ + 0 due to Jordan-Wigner transformation. Despite the nonlocality, the Heisenberg evolution of this operator can be simulated efficiently [14,19]. We actually simulate Majorana operators w = f ↑ + f † ↑ and w = i(f ↑ − f † ↑ ) which are hermitian and allow the simulation in the real arithmetics. The results for the imaginary part of the Green's function, shown in Fig. 3, can be used to calculate the spectral functions (see e.g. Ref. [27,28,29]) by means of a Fourier transform which will be in detail presented elsewhere.

Conclusions
We have proposed a tensor network method to simulate the real time dynamics of one dimensional quantum many-body systems at finite temperatures. The main advantage of the method is that it decouples the real time evolution from the imaginary evolution in the inverse temperature such that two independent simulations give the results for all combinations temperature-time. Furthermore, the simulations are done in real arithmetics. The method is most useful for integrable systems  where the costs grow polynomially, and systems close to the integrability where the asymptotic exponential growth only appears at later times. The method can also be used to study nonequilibrium impurity physics [30], in particular current-voltage characteristics of strongly correlated nanostructures in the non-linear response regime. The proposed framework can be extended to two dimensional systems described in terms of projected entangled pair states [31] (see Appendix B); albeit with significantly higher computational costs. (v) Construct a MPS representation for the initial operator |a) (vi) Simulate |ρ(β)) = e −βχ |e) using tDMRG and save |ρ(β)) for various inverse temperatures β (vii) Simulate |a(t)) = e −itĤ |a) using tDMRG and save |a(t)) for various times t (viii) Constructb from b such thatb : x → bx (in the same way asχ from H).
(ix) For all t and β, calculate (ρ(β)|b|a(t)) and (ρ(β)|e). where H The thermal super-mapχ thus contains the same number of product operators as the hamiltonian where single-site terms are transformed to j, l = 0, 1, . . . , d 2 − 1 and the two-site terms to ]. The super-mapĤ contains twice as many local terms as the Hamiltonian operator (due to the commutator), namelŷ ).
The single-site terms are obtained as Real time evolution at finite temperatures with operator space matrix product states 9 and similarly for the two-site terms.
Note that we can and we should (to work in real arithmetics) use a different basis set {σ ν } forχ andĤ,real matrices for the former and Hermitian for the latter.
The transformation from the Hamiltonian toχ andĤ can thus be treated as a black-box routine, providing the "Hamiltonian operators" which can be used in an imaginary (thermal) and a real (Heisenberg picture) time evolution -using an existing implementation of the time dependent DMRG algorithm (tDMRG) [5].

Appendix B. Two-dimensional systems with PEPS
Here we briefly sketch the way how the method can be extended to two dimensional systems. In principle, the method can already be used for small two-dimensional systems as-is, represented by a matrix product state in a snake-like structure [1,5] which gives remarkably accurate results for the ground states, see e.g. Ref. [32] for a recent study. However, for more extended two-dimensional systems the preferred description is given in terms of Projected Entangled Pair States (PEPS) [31]. Here we show how our method can be used with PEPS. For simplicity we assume a spin-1/2 quantum system on a m×n rectangular lattice. Any operator can be written as a superposition of basis operators a = ν i,j a ν σ ν 1,1 1,1 · · · σ ν 1,n 1,n · · · σ νm,n m,n for ν i,j ∈ {0, x, y, z} and the corresponding elements of the operator space |a) can be represented by a PEPS-Ansatz, |a) = ν 1,1 ···νm,n tr[A [1,1]ν 1,1 · · · A [m,n]νm,n ]|ν 1,1 , . . . , ν m,n ). (B.1) The trace in (B.1) must be understood as a tensor-trace, for details on PEPS see Ref. [31]. The operator a is now represented in terms of rank-4 tensors A [i,j]ν i,j ∈ R D l ×Du×D d ×Dr where (D l , D u , D d , D r ), are the dimensions of the bonds connecting the site [i, j] to the neighboring sites (left,up,down,right), respectively. The method proposed in the main text can now be used exactly in the same way as in the one-dimensional case, by operating on PEPS instead of MPS. The super-mapŝ χ andĤ are identical to the one-dimensional case and generate the evolution in the Heisenberg picture |a(t)) = e −itĤ |a) and the thermal evolution of the density matrix |ρ(β)) = e −βχ |e).
The only difference to the one-dimensional case is that the computation with PEPS is significantly more expensive and an accurate description might require a considerably increased computational effort. In particular, the complexity of the real time evolution scales as O(D 12 ) (as opposed to O(D 3 ) for matrix product states). It is not clear at present whether the computational costs still scale polynomially in time like in the onedimensional case, however, results from the ground state simulations are consistent with this conjecture.