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

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 (MPSs) [2–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 [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 MPSs. 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 [10, 14]. 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.

Let us consider a quantum system of n sites with a Hilbert space $\mathcal{H}$. The expectation value of an operator $a:\mathcal{H}\to \mathcal{H}$ in a mixed state described by a density matrix ρ, is given by $\langle a\rangle ={\rm tr}\left( \rho a \right)$ with ${\rm tr}\left( \rho \right)=1$. Like the operator a, the density matrix is also a linear map over $\mathcal{H}$ and these maps form the operator space algebra $\mathcal{K}$ with an inner product $\left( a|b \right)={{\left( {\rm dim}\mathcal{H} \right)}^{-1}}{\rm tr}\left( {{a}^{\dagger }}b \right)$. The expectation value of a can thus be given as

$$\begin{eqnarray}&&\langle a\rangle =\left( \rho |a \right)/\left( \rho |e \right),\end{eqnarray} \tag{ 1 }$$

where the bracket notation ${\rm | }a{\rm )}$ is used when referring to $a\in \mathcal{K}$ and ${\rm | }e{\rm )}$ corresponds to the fully mixed state (the identity). In thermal equilibrium for inverse temperature $\beta ={{\left( {{k}_{{\rm B}}}T \right)}^{-1}}$, the expectation values are given by using ${\rm | }\rho \left( \beta \right){\rm )}={\rm | }{{e}^{-\beta H}}{\rm )}$ as the density matrix whereas the time evolution of equation (1) is obtained by replacing a by the Heisenberg operator $a\left( t \right)={{e}^{itH}}a{{e}^{-itH}}$. Clearly, the simulations for $\rho \left( \beta \right)$ and a(t) can be done independently.

The density matrices $\rho \left( \beta \right)$ are obtained by the imaginary time evolution, starting from the infinite temperature equilibrium state ${\rm | }e{\rm )}$. By defining a map over the operator space, $\hat{\chi }:\mathcal{K}\to \mathcal{K}$, which left-multiplies an operator by H, i.e., $\hat{\chi }{\rm | }a{\rm )}\equiv {\rm | }Ha{\rm )}$, the density matrix ${\rm | }\rho \left( \beta \right){\rm )}={\rm | }{{e}^{-\beta H}}{\rm )}$ becomes

$$\begin{eqnarray}&&{\rm | }\rho \left( \beta \right){\rm )}={{e}^{-\beta \hat{\chi }}}{\rm | }e{\rm )}.\end{eqnarray} \tag{ 2 }$$

Similarly, the Heisenberg evolution $a\left( t \right)={{e}^{itH}}a{{e}^{-itH}}$ can be simulated by introducing a linear map defined as $\hat{H}{\rm | }x{\rm )}\equiv {\rm | }\left[ x,H \right]{\rm )}$ for $x\in \mathcal{K}$. It was shown [19] that the Heisenberg evolution is hence transformed to the Schrödinger evolution in the operator space

$$\begin{eqnarray}&&{\rm | }a\left( t \right){\rm )}={{e}^{-it\hat{H}}}{\rm | }a\left( 0 \right){\rm )}.\end{eqnarray} \tag{ 3 }$$

Both $\hat{\chi }$ and $\hat{H}$ have the same form and the same range of interactions as H, if the basis for $\mathcal{K}$ is local (see appendix A).

Let us now define a suitable basis for $\mathcal{K}$. While any basis could in principle be used, it is advantageous to choose a basis given by direct products of local basis operators $\left\{ p_{{{\nu }_{j}}}^{\left[ j \right]} \right\}$ pertaining to the site j, as ${{P}_{\mathop{\nu }\limits_{\_}}}=p_{{{\nu }_{1}}}^{\left[ 1 \right]}\otimes p_{{{\nu }_{2}}}^{\left[ 2 \right]}\otimes \cdots \otimes p_{{{\nu }_{n}}}^{\left[ n \right]}$. These operators are chosen such that they form an orthonormal set. A suitable choice for $\left\{ p_{{{\nu }_{j}}}^{\left[ j \right]} \right\}$ 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 $\hat{\chi }$, represented by ${{\left[ {\hat{\chi }} \right]}_{\mu ,\nu }}={{\left( {\rm dim}\mathcal{H} \right)}^{-1}}{\rm tr}\left( P_{\mu }^{\dagger }\hat{\chi }{{P}_{\nu }} \right)$, is real if $\left\{ {{P}_{\nu }} \right\}$ are real, which in turn allows us to simulate the thermal density matrices (2) with real arithmetics.

For the real time evolution (3), where $\hat{H}$ is represented by a matrix ${{\left[ {\hat{H}} \right]}_{\mu ,\nu }}={{\left( {\rm dim}\mathcal{H} \right)}^{-1}}{\rm tr}\left( \left[ P_{\mu }^{\dagger },{{P}_{\nu }} \right]H \right)$, a different choice for $\left\{ {{P}_{\nu }} \right\}$ is made to maintain the formulation in real arithmetics. If $P_{\mu }^{\dagger }={{P}_{\mu }}$, then ${{\left[ {\hat{H}} \right]}_{\mu ,\nu }}=-{{\left[ {\hat{H}} \right]}_{\nu ,\mu }}$ and $\hat{H}$ 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_{{{\mu }_{j}}}^{\left[ j \right]}=\sigma _{j}^{{{\mu }_{j}}}$ for ${{\mu }_{j}}\in \left\{ 0,1,2,3 \right\}$. This basis is Hermitian and thus suitable for the real time evolution, whereas the thermal evolution is done using a real local basis $\left\{ {{\sigma }^{0}},{{\sigma }^{1}},-i{{\sigma }^{2}},{{\sigma }^{3}} \right\}$. The local basis for spin-1 systems is provided by Gell–Mann matrices, whereas for spin-$3/2$ it can be chosen as $\left\{ {{S}^{\left( a,b \right)}}={{\sigma }^{a}}\otimes {{\sigma }^{b}} \right\}$. 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}_{{{\nu }_{1}},\ldots ,{{\nu }_{n}}}}=\sigma _{1}^{{{\nu }_{1}}}\cdots \sigma _{n}^{{{\nu }_{n}}}$ form the basis set ${\rm | }\mathop{\nu }\limits_{\_}{\rm )}\equiv {\rm | }{{P}_{\mathop{\nu }\limits_{\_}}}{\rm )}$ and $a\in \mathcal{K}$ is given in terms of an MPS [14]

$$\begin{eqnarray*}&&{\rm | }a{\rm )}=\mathop{\sum }\limits_{{{\nu }_{1}},\ldots ,{{\nu }_{n}}}{\rm tr}\left[ {{{\bf A}}^{\left[ 1 \right]{{\nu }_{1}}}}\cdots {{{\bf A}}^{\left[ n \right]{{\nu }_{n}}}} \right]{\rm | }{{\nu }_{1}},\ldots ,{{\nu }_{n}}{\rm )},\end{eqnarray*}$$

with real matrices ${{{\bf A}}^{\left[ j \right]{{\nu }_{j}}}}$ of dimension at most D×D (see e.g., [4]) where D is the bond dimension. The evolutions (2) and (3) are simulated using the Suzuki–Trotter decomposition and local updates of the MPS, known as the time dependent DMRG algorithm [5–8], which is particularly efficient when H contains at most nearest-neighbor interactions.

The thermal expectation value ${{\langle A\left( t \right)\rangle }_{\beta }}$ given by (1) is obtained by combining the independent simulations (2) and (3). Using the same principles, we can simulate ${{\left\langle a\left( t \right)b \right\rangle }_{\beta }}$ where a(t) and b are MPS, by constructing a matrix product operator $\hat{B}:x\mapsto bx$. The result is given by ${{\left\langle b\;a\left( t \right) \right\rangle }_{\beta }}={\rm (}\rho \left( \beta \right){\rm | }\hat{B}{\rm | }a\left( t \right){\rm )}/\left( \rho \left( \beta \right)|e \right)$ which in the MPS language corresponds to calculating an expectation value. Similarly ${{\left\langle a\left( t \right)b \right\rangle }_{\beta }}$ is obtained.

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}^{\sharp }}\left( x \right)=-{\rm tr}\left[ \hat{R}{{{\rm log} }_{2}}\hat{R} \right]$ for $\hat{R}={\rm t}{{{\rm r}}_{E}}\left( {\rm | }x{\rm )(}x{\rm | } \right)$. 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 ${\rm exp} \left( {{S}^{\sharp }}\left( x \right) \right)$ 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

$$\begin{eqnarray}&&H=\sum \limits_{l=1}^{n-1}\left( \sigma _{l}^{x}\sigma _{l+1}^{x}+\sigma _{l}^{y}\sigma _{l+1}^{y}+\Delta \sigma _{l}^{z}\sigma _{l+1}^{z} \right),\end{eqnarray} \tag{ 4 }$$

where $\Delta$ is the anisotropy parameter. The OSEE of the initial state $\rho \left( \beta =0 \right)$ is zero (the fully mixed state is separable). As shown in figure 1 (left), we observe logarithmic growth of the OSEE in the inverse temperature [17]. Furthermore, we simulate the Heisenberg evolution (3) of the spin current flowing between sites m and $m+1$, ${{j}_{m}}=i\left( \sigma _{m}^{+}\sigma _{m+1}^{-}-\sigma _{m}^{-}\sigma _{m+1}^{+} \right)$, with $\sigma _{m}^{\pm }=\left( \sigma _{m}^{x}\pm i\sigma _{m}^{y} \right)/2$ for $m=n/2$. The results in figure 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 $\Delta$ 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, as in other similar methods.

Zoom In Zoom Out Reset image size

The simulation of the density matrices and the local operators can now be used to calculate e.g., the spin current correlation function ${{\left\langle j{{j}_{n/2}}\left( t \right) \right\rangle }_{\beta }}$. The results are shown in figure 2, with an inverse temperature $\beta =0.25$. The infinite-time limit of these current-correlations 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 [13, 23, 24], however, a detailed discussion is beyond the scope of this manuscript and will be presented elsewhere. When compared to [13] (with our time scale multiplied by a factor of four due to the spin operators) we observe, on one hand, that the reachable times are approximately twice as large in the gapped regime ($\Delta =2$ in figure 2) with a maximal bond dimension D = 1000. On the other hand, the computation in the gapless regime, $\Delta <1$, is less demanding and the curves obtained with a smaller bond dimension D = 800 are essentially indistinguishable in figure 2. This is in complete agreement with the picture obtained from the OSEE scaling. In turn, for small $\Delta$ the simulation can be pushed beyond the time scale shown in figure 2, similarly to [13].

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&H=\sum \limits_{j=-\infty }^{\infty }{{\tau }_{j}}\left( \sigma _{j}^{+}\sigma _{j+1}^{-}+\sigma _{j}^{-}\sigma _{j+1}^{+} \right)+U{{n}_{0}}{{n}_{1}}+{{\epsilon }_{f}}\left( {{n}_{0}}+{{n}_{1}} \right)\end{eqnarray} \tag{ 5 }$$

with ${{n}_{j}}=\sigma _{j}^{-}\sigma _{j}^{+}$ and the interaction strength U. The hopping parameters ${{\tau }_{-1}}={{\tau }_{1}}$ are given by the hybridization whereas the rest of ${{\tau }_{j}}$ are determined from the discretization of the bath (and ${{\tau }_{0}}=0$), see [26] for details. The interaction term is only present between sites 0 and 1 (corresponding to spin-$\uparrow$ and spin-$\downarrow$ of the impurity) which allows longer acessible times and lower temperatures. We simulate the single particle Greenʼs function $G\left( t \right)=-i{{\left\langle \left\{ f_{\uparrow }^{\dagger },{{f}_{\uparrow }}\left( t \right) \right\} \right\rangle }_{\beta }}$ where ${{f}_{\uparrow }}$ is the annihilation operator at the impurity for spin-$\uparrow$, written as ${{f}_{\uparrow }}=\left( {{\prod }_{j<0}}\sigma _{j}^{z} \right)\sigma _{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}_{\uparrow }}+f_{\uparrow }^{\dagger }$ and $w^{\prime} =i\left( {{f}_{\uparrow }}-f_{\uparrow }^{\dagger } \right)$ which are Hermitian and allow the simulation in the real arithmetics. The results for the imaginary part of the Greenʼs function, shown in figure 3, can be used to calculate the spectral functions (see e.g., [27–29]) by means of a Fourier transform, which will be presented in detail elsewhere.

Zoom In Zoom Out Reset image size

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 of temperature and 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 the current–voltage characteristics of strongly correlated nanostructures in the nonlinear 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.

We thank F Verstraete, Lei Wang and C Karrasch for valuable discussions and F Heidrich-Meisner, T Prosen and M Žnidarič for their comments on the manuscript. IP acknowledges the collaboration on related projects with T Prosen and M Žnidarič. SA and VE thank M Rams and B Terhal for discussions. This work was supported by the Swiss National Science Foundation through the NCCR QSIT, by the European Research Council through ERC grants QUERG and SIMCOFE and by SFB ViCoM (FWF project ID F4104). Simulations were run on the Brutus cluster at ETH Zurich and on the Vienna Scientific Cluster.

In table A.1 we present an algorithm to simulate the real time evolution at finite temperatures. In particular, we explain how to generate the super operators $\hat{\chi }$ and $\hat{H}$ for a one-dimensional system of n sites with a local dimension d and a Hamiltonian operator H. Let us first decompose the Hamiltonian into a sum of local product operators

$$\begin{eqnarray}&&H=\sum \limits_{\mu }H_{1}^{\left( \mu \right)}+\sum \limits_{\nu }H_{2}^{\left( \nu \right)},\end{eqnarray} \tag{ A.1 }$$

where $H_{1}^{\left( \mu \right)}$ are local operators of dimension $d\times d$ and $H_{2}^{\left( \nu \right)}$ are two-site operators of dimension $\left( d\times d \right)\times \left( d\times d \right)$. The local Hilbert space for one site is generated by a set of ${{d}^{2}}$ Gell–Mann matrices $\left\{ {{\sigma }^{j}},\;j=0,1,\ldots ,{{d}^{2}}-1 \right\}$ of dimension $d\times d$ where ${{\sigma }^{0}}={\bf 1}$ is the identity. The super operators $\hat{\chi }$ and $\hat{H}$ are defined as $\hat{\chi }:a\mapsto Ha$ and $\hat{H}:a\mapsto \left[ a,H \right]$.

Table A.1.  Practical algorithm to simulate a time correlation function ${{\left\langle ba\left( t \right) \right\rangle }_{\beta }}$ at an inverse temperature β.

Practical algorithm for the calculation of ${{\left\langle ba\left( t \right) \right\rangle }_{\beta }}$
(i) Assume a one-dimensional system of n sites with a local dimension d and Hamiltonian H of the form (A.1)
(ii) Set up Gell–Mann matrices $\left\{ {{\sigma }^{j}},\;j=0,1,\ldots ,{{d}^{2}}-1 \right\}$
(iii) Construct super-maps $\hat{\chi }$ and $\hat{H}$
(iv) Construct an MPS representation of the the unit ${\rm | }e{\rm )}$
(v) Construct an MPS representation for the initial operator ${\rm | }a{\rm )}$
(vi) Simulate ${\rm | }\rho \left( \beta \right){\rm )}={{{\rm e}}^{-\beta \hat{\chi }}}{\rm | }e{\rm )}$ using tDMRG and save ${\rm | }\rho \left( \beta \right){\rm )}$ for various inverse temperatures β
(vii) Simulate ${\rm | }a\left( t \right){\rm )}={{{\rm e}}^{-it\hat{H}}}{\rm | }a{\rm )}$ using tDMRG and save ${\rm | }a\left( t \right){\rm )}$ for various times t
(viii) Construct $\hat{b}$ from b such that $\hat{b}:x\mapsto bx$ (in the same way as $\hat{\chi }$ from H)
(ix) For all t and β, calculate ${\rm (}\rho \left( \beta \right){\rm | }\hat{b}{\rm | }a\left( t \right){\rm )}$ and $\left( \rho \left( \beta \right)|e \right)$

The thermal super-map $\hat{\chi }$ thus contains the same number of product operators as the Hamiltonian

$$\begin{eqnarray*}&&\hat{\chi }=\sum \limits_{\mu }\hat{\chi }_{1}^{\left( \mu \right)}+\sum \limits_{\nu }\hat{\chi }_{2}^{\left( \nu \right)},\end{eqnarray*}$$

where single-site terms are transformed to ${{d}^{2}}\times {{d}^{2}}$ matrices

$$\begin{eqnarray*}&&{{\left[ \hat{\chi }_{1}^{\left( \mu \right)} \right]}_{j,l}}={{d}^{-1}}{\rm tr}\left[ {{\left( {{\sigma }^{j}} \right)}^{\dagger }}H_{1}^{\left( \mu \right)}{{\sigma }^{l}} \right],j,l=0,1,\ldots ,{{d}^{2}}-1,\end{eqnarray*}$$

and the two-site terms to ${{d}^{2}}\times {{d}^{2}}\times {{d}^{2}}\times {{d}^{2}}$ tensors

$$\begin{eqnarray*}&&{{\left[ \hat{\chi }_{2}^{\left( \nu \right)} \right]}_{{{j}_{1}},{{j}_{2}},{{l}_{1}},{{l}_{2}}}}={{d}^{-2}}{\rm tr}\left[ {{\left( {{\sigma }^{{{j}_{1}}}}\otimes {{\sigma }^{{{j}_{2}}}} \right)}^{\dagger }}H_{2}^{\left( \nu \right)}\left( {{\sigma }^{{{l}_{1}}}}\otimes {{\sigma }^{{{l}_{2}}}} \right) \right].\end{eqnarray*}$$

The super-map $\hat{H}$ contains twice as many local terms as the Hamiltonian operator (due to the commutator), namely

$$\begin{eqnarray*}&&\hat{H}=\sum \limits_{\mu }\left( \hat{H}_{1}^{\left( 2\mu -1 \right)}+\hat{H}_{1}^{\left( 2\mu \right)} \right)+\sum \limits_{\nu }\left( \hat{H}_{2}^{\left( 2\nu -1 \right)}+\hat{H}_{2}^{\left( 2\nu \right)} \right).\end{eqnarray*}$$

The single-site terms are obtained as

$$\begin{eqnarray*}&&{{\left[ \hat{H}_{1}^{\left( 2\mu -1 \right)} \right]}_{j,l}}={{d}^{-1}}{\rm tr}\left[ {{\left( {{\sigma }^{j}} \right)}^{\dagger }}{{\sigma }^{l}}{{H}^{\left( \mu \right)}} \right]\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{{\left[ \hat{H}_{1}^{\left( 2\mu \right)} \right]}_{j,l}}={{d}^{-1}}{\rm tr}\left[ {{\left( {{\sigma }^{j}} \right)}^{\dagger }}{{H}^{\left( \mu \right)}}{{\sigma }^{l}} \right].\end{eqnarray*}$$

and similarly for the two-site terms.

Note that we can and we should (to work in real arithmetics) use a different basis set $\left\{ {{\sigma }_{\nu }} \right\}$ for $\hat{\chi }$ and $\hat{H}$, real matrices for the former and Hermitian for the latter.

The transformation from the Hamiltonian to $\hat{\chi }$ and $\hat{H}$ 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].

Here we briefly sketch 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., [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\times n$ rectangular lattice. Any operator can be written as a superposition of basis operators

$$\begin{eqnarray*}&&a=\sum \limits_{{{\nu }_{i,j}}}{{a}_{\mathop{\nu }\limits_{\_}}}\;\sigma _{1,1}^{{{\nu }_{1,1}}}\cdots \sigma _{1,n}^{{{\nu }_{1,n}}}\cdots \sigma _{m,n}^{{{\nu }_{m,n}}}\end{eqnarray*}$$

for ${{\nu }_{i,j}}\in \left\{ 0,x,y,z \right\}$ and the corresponding elements of the operator space ${\rm | }a{\rm )}$ can be represented by a PEPS-Ansatz,

$$\begin{eqnarray}&&{\rm | }a{\rm )}=\sum \limits_{{{\nu }_{1,1}}\cdots {{\nu }_{m,n}}}{\rm tr}\left[ {{{\bf A}}^{\left[ 1,1 \right]{{\nu }_{1,1}}}}\cdots {{{\bf A}}^{\left[ m,n \right]{{\nu }_{m,n}}}} \right]{\rm | }{{\nu }_{1,1}},\ldots ,{{\nu }_{m,n}}{\rm )}.\end{eqnarray} \tag{ B.1 }$$

The trace in (B.1) must be understood as a tensor trace, for details on PEPS see [31]. The operator a is now represented in terms of rank-4 tensors ${{{\bf A}}^{\left[ i,j \right]{{\nu }_{i,j}}}}\in {{\mathbb{R}}^{{{D}_{{\rm l}}}\times {{D}_{{\rm u}}}\times {{D}_{{\rm d}}}\times {{D}_{{\rm r}}}}}$ where $\left( {{D}_{{\rm l}}},{{D}_{{\rm u}}},{{D}_{{\rm d}}},{{D}_{{\rm r}}} \right)$, are the dimensions of the bonds connecting the site $\left[ i,j \right]$ 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-maps $\hat{\chi }$ and $\hat{H}$ are identical to the one-dimensional case and generate the evolution in the Heisenberg picture

$$\begin{eqnarray*}&&{\rm | }a\left( t \right){\rm )}={{{\rm e}}^{-it\hat{a}}}{\rm | }a{\rm )},\end{eqnarray*}$$

and the thermal evolution of the density matrix

$$\begin{eqnarray*}&&{\rm | }\rho \left( \beta \right){\rm )}={{{\rm e}}^{-\beta \hat{\chi }}}{\rm | }e{\rm )}.\end{eqnarray*}$$

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\left( {{D}^{12}} \right)$ (as opposed to $O\left( {{D}^{3}} \right)$ for MPS). It is not clear at present whether the computational costs still scale polynomially in time like in the one-dimensional case, however, results from the ground state simulations are consistent with this conjecture.