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

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

$$\begin{equation} h_{n} = J_n \big(S^x_{n}S^x_{n+1} + S^y_{n}S^y_{n+1} + \Delta_n S^z_{n}S^z_{n+1}\big) + b_n (S_n^z-S_{n+1}^z) , \end{equation} \tag{ 8 }$$

or equivalently spinless fermions through a Jordan–Wigner transformation. By choosing the couplings J_n, Δ_n and b_n appropriately:

$$\begin{equation} J_n = \left\{ \begin{array}{@{}ll@{}} 1, & n\,{\rm odd}, \\ \lambda, & n\,{\rm even}, \end{array} \right.\quad \Delta_n = \Delta,\quad b_n= \frac{(-1)^n b}{2}, \end{equation} \tag{ 9 }$$

we can study systems which are gapless or gapped, and investigate the role of integrability. For λ = 1 and b = 0, equation (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 λ < λ_c or b > b_c, where λ_c < 1 and b_c > 0 if $-1<\Delta <-1/\sqrt {2}$ [42, 80, 81]. In addition, we study the quantum Ising model

$$\begin{equation} h_n = - 4 S^z_{n}S^z_{n+1} - g( S_n^x+S_{n+1}^x ). \end{equation} \tag{ 10 }$$

2.2.1. Linear response

We first consider a homogeneous system of size L governed by $H = \sum _{n=1}^{L-1} h_n$ within linear response. The optical charge (C) and energy (E) conductivities can be computed via the Kubo formula

$$\begin{equation} \sigma(\omega) = \frac{1}{\omega L}\int_0^\infty {\rm e}^{{\rm i}\omega t} \langle [J(t),J]\rangle \mathrm{d}t = \frac{1-\mathrm{e}^{-\omega/T}}{\omega L} \int_0^\infty {\rm e}^{{\rm i}\omega t} \langle J(t)J\rangle\mathrm{ d}t, \end{equation} \tag{ 11 }$$

where the corresponding current operators $J\!=\!\sum _n j_n$ are defined through a continuity equation⁴:

$$\begin{eqnarray} \partial_t h_n &=& j_{{E},n}-j_{{ E},n+1} ~\Rightarrow~ J_{E} = \mathrm{i} \sum_{n=2}^{L-1} [h_{n-1},h_n],\nonumber\\ \partial_t S_n^z & = & j_{{C},n}-j_{{C},n+1} ~\Rightarrow~ J_{C} = \mathrm{i} \sum_{n=2}^{L-1} [h_{n-1},S_n^z]. \end{eqnarray} \tag{ 12 }$$

One usually decomposes the real part of σ(ω) as

$$\begin{equation} {\rm Re }\sigma(\omega) = 2\pi D \delta(\omega) + \sigma_{{\rm reg}}(\omega), \end{equation} \tag{ 13 }$$

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

$$\begin{equation} D = \lim_{t\to\infty}\lim_{{\rm L}\to\infty} \frac {{\rm Re } \langle J(t)J\rangle}{2LT}. \end{equation} \tag{ 14 }$$

As already stated above, a time-dependent equilibrium correlation function is defined as

$$\begin{equation} \langle A(t) B \rangle = {\rm Tr} \left( \rho_T\, {\rm e}^{{\rm i}Ht}A\,{\rm e}^{-{\rm i}Ht}B \right),~\rho_T = \frac{\mathrm{e}^{-H/T}}{Z_T} \end{equation} \tag{ 15 }$$

with $Z_T={\mathrm { Tr}}\exp (-H/T)$ denoting the partition function.

2.2.2. Non-equilibrium

In addition to linear response, we study two thermal non-equilibrium 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

$$\begin{equation} H_0 = H_{\rm L} + H_{\rm R} = \sum_{n=1}^{L/2-1} h_n + \sum_{n=L/2+1}^{L-1} h_n , \end{equation} \tag{ 16 }$$

each being in thermal (grand-canonical) equilibrium at temperatures T_L and T_R. The corresponding density matrix factorizes

$$\begin{equation} \rho_0=\rho_{\rm L}\otimes\rho_{\rm R},\quad \rho_i=\frac{\exp(-H_i/T_i)}{{\rm Tr}\exp(-H_i/T_i)},\quad i=\mathrm{L,R}. \end{equation} \tag{ 17 }$$

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₀ + h_L/2:

$$\begin{equation} \langle A(t)\rangle = {\rm Tr}\, \rho(t) A,~ \rho(t) = {\rm e}^{{\rm i}Ht}\rho_0\,{\rm e}^{-{\rm i}Ht}. \end{equation} \tag{ 18 }$$

In the second setup ('non-equilibrium B'; see figure 1(B)) we investigate two non-interacting chains Δ = b = 0, λ = 1 of length L/2,

$$\begin{equation} H_0 = H_{\rm L} + H_{\rm R} = \sum_{n=1}^{L/2-1} h_n + \sum_{n=L/2+2}^{L} h_n, \end{equation} \tag{ 19 }$$

at different temperatures T_L and T_R. At t = 0, they are coupled via an interacting resonant level model:

$$\begin{eqnarray}\fl h_{\rm IRLM} &=& t' \big(S^x_{{L}/2}S^x_{{L}/2+1} + S^y_{{L}/2}S^y_{{L}/2+1} + U S^z_{{L}/2}S^z_{{L}/2+1}\big) \nonumber\\ &+& t' \big(S^x_{{L}/2+1}S^x_{{L}/2+2} + S^y_{{L}/2+1}S^y_{{L}/2+2} + U S^z_{{L}/2+1}S^z_{{L}/2+2}\big). \end{eqnarray} \tag{ 20 }$$

The site n = L/2 + 1 is initially in an equal superposition of up and down states (i.e. formally at infinite temperature).

Zoom In Zoom Out Reset image size

In this section, we give a brief overview of the DMRG method [1–4]. More details can be found in the appendix. In order to evaluate equations (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, |Ψ_T 〉 = e^−H/(2T)|Ψ_∞〉 [26, 37, 38]. The correlation function of equation (15) is exactly recast as

$$\begin{equation} \langle A(t)B\rangle = \frac{\langle \Psi_\infty|\mathrm{ e}^{-\frac{H}{2T}}\,{\rm e}^{{\rm i}Ht} A \,{\rm e}^{-{\rm i}Ht} B \,\mathrm{e}^{-\frac{H}{2T}}|\Psi_\infty\rangle}{\langle\Psi_\infty|\mathrm{ e}^{-\frac{H}{T}}|\Psi_\infty\rangle}, \end{equation} \tag{ 21 }$$

and this object is directly accessible in standard time-dependent DMRG frameworks [17–22]. It is convenient to first express the initial state |Ψ_∞〉 in terms of a matrix product state [5–8],

$$\begin{equation} |\Psi_\infty\rangle = \sum_{\sigma_n,\sigma_{n_Q}} A^{\sigma_1}A^{\sigma_{1_Q}}\cdots A^{\sigma_{{L}}}A^{\sigma_{{L}_Q}}|\sigma_1\sigma_{1_Q}\ldots\sigma_{{L}}\sigma_{{L}_Q}\rangle, \end{equation} \tag{ 22 }$$

where

$$\begin{equation} A^{\sigma_i}_{a_ia_{i+1}} = \Lambda_{a_i}^i\Gamma_{a_ia_{i+1}}^{\sigma_i}. \end{equation} \tag{ 23 }$$

Here and in the following σ_i is a short hand for either a physical or auxiliary degrees of freedom: σ_i∈{σ_n,σ_nQ}. The initial matrices associated with |Ψ_∞〉 read

$$\begin{eqnarray} \Gamma^{\sigma_{n_{\phantom{Q}}}=\uparrow}& = (1~0), \qquad \qquad \Gamma^{\sigma_{n_{\phantom{Q}}}=\downarrow} = (0~-1), \nonumber \\ \Gamma^{\sigma_{n_Q}=\uparrow}& = (0~1/\sqrt{2})^{\mathrm{T}}, \qquad \Gamma^{\sigma_{n_Q}=\downarrow} = (1/\sqrt{2}~0)^T \end{eqnarray} \tag{ 24 }$$

as well as Λⁱ = 1. After factorizing the evolution operators $\exp (-\lambda H)$ using a second or fourth order Trotter decomposition, they can be successively applied to equation (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 [40]. One has the analytic freedom to apply any time-dependent unitary transformation

$$\begin{equation} U_Q(t) = \sum_{\sigma_{1_Q}\ldots \sigma_{{L}_Q}\atop\sigma_{1_Q}'\ldots \sigma_{{L}_Q}'} C_{\sigma_{1_Q}'\ldots \sigma_{{L}_Q}'}^{\sigma_{1_Q}\ldots \sigma_{{L}_Q}}(t) |\sigma_{1_Q}\ldots \sigma_{{\rm L}_Q}\rangle\langle \sigma_{1_Q}'\ldots \sigma_{{L}_Q}'| \end{equation} \tag{ 25 }$$

to the in principle inert auxiliary sites σ_nQ—physical quantities are determined by the trace over Q and are thus not affected by U_Q(t):

$$\begin{equation} [U_Q(t), \sigma_n] = 0. \end{equation} \tag{ 26 }$$

For example, equation (21) can be rewritten as

$$\begin{equation} \langle \Psi_T|{\rm e}^{{\rm i}Ht} A\, {\rm e}^{-{\rm i}Ht} \,{\rm e}^{{\rm i}Ht'} B\, {\rm e}^{{\rm i}Ht'} |\Psi_T\rangle = \langle \Psi_T| U^\dagger (t) A U(t) U^\dagger(t') B U(t')|\Psi_T\rangle, \end{equation} \tag{ 27 }$$

where

$$\begin{equation} U(t) = {\rm e}^{-{\rm i}Ht} U_Q(t). \end{equation} \tag{ 28 }$$

Put differently, purification is not unique. It turned out [31, 32, 40, 42, 43] that choosing

$$\begin{equation} U_Q(t) = \mathrm{e}^{+\mathrm{i}\tilde Ht},\quad \tilde H = H\big(\sigma_n\to \sigma_{n_Q},~{\rm magnetic\,\,fields\,\,reversed}\big), \end{equation} \tag{ 29 }$$

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

$$\begin{equation} {\rm e}^{-{\rm i}Ht}\,{\rm e}^{{\rm i}\tilde Ht} = {\rm e}^{-{\rm i}H\Delta t}\,{\rm e}^{{\rm i}\tilde H\Delta t}\,{\rm e}^{-{\rm i}H\Delta t}\,{\rm e}^{{\rm i}\tilde H\Delta t}\ldots. \end{equation} \tag{ 30 }$$

It is instructive to first consider the trivial case of A = B = 1 in equation (21), i.e.

$$\begin{equation} 1 =\frac{\langle\Psi_T|{\rm e}^{{\rm i}Ht}\, {\rm e}^{-{\rm i}Ht}| \Psi_T\rangle}{\langle\Psi_T|\Psi_T\rangle} = \frac{\langle\Psi_T|{\rm e}^{{\rm i}Ht}U_Q^\dagger(t)\, {\rm e}^{-{\rm i}Ht}U_Q(t)|\Psi_T\rangle}{\langle\Psi_T|\Psi_T\rangle}. \end{equation} \tag{ 31 }$$

Figure 2 shows the evolution of the bond dimension χ during the calculation of e^−iHtU_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. Choosing $U_Q(t)={\mathrm { e}}^{{\mathrm { i}}\tilde Ht}$ removes this artifact. This immediately indicates why longer time scales can be reached in DMRG evaluations of, e.g. current correlation functions e^−iHtU_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).

Zoom In Zoom Out Reset image size

An intuitive understanding of why the particular choice of $U_Q(t) = {\mathrm { e}}^{+{\mathrm { i}}\tilde Ht}$ removes the 'artifical' entanglement growth was recently provided in [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 equation (27). If A (and likewise for B) is local, most terms in the first Trotter step e^iHΔtA e^−iHΔt cancel out. As time evolves,

$$\begin{equation} \ldots {\rm e}^{{\rm i}H\Delta t}\,{\rm e}^{{\rm i}H\Delta t}A\,{\rm e}^{-{\rm i}H\Delta t}\,{\rm e}^{-{\rm i}H\Delta t}\ldots, \end{equation} \tag{ 32 }$$

entanglement can only build-up gradually around the region on which A acted due to the so-called light-cone effect in non-relativistic systems: Lieb–Robinson bounds [84] generically stipulate that correlation functions 〈A(t)B〉 of local operators A and B fall 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.

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 equation (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 non-integrable dimerized chain supports dissipationless transport or that its current decays on a hidden large scale; see [43] for further details). The steady-state current of the XXZ chain is of the simple 'black-body' form [43, 76, 77, 85]

$$\begin{equation} \lim_{t\to\infty} \langle j_{{E},n}(t)\rangle = f(T_{L}) - f(T_{R}), \end{equation} \tag{ 33 }$$

implying that it is determined solely by the linear conductance G(T) ∼ ∂_T f(T) for arbitrary large T_L − T_R:

$$\begin{equation} \lim_{t\to\infty} \langle j_{{\rm E},n}(t)\rangle \sim \int_{T_{L}}^{T_{R}} G(T)\, \mathrm{d}T. \end{equation} \tag{ 34 }$$

This agrees with a field theory result [77] at low temperatures; asymptotic low-T behavior sets in around T ≲ 0.2.

Zoom In Zoom Out Reset image size

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 equation (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 χ³.

Zoom In Zoom Out Reset image size

The same reduction of χ holds for the non-equilibrium impurity problem depicted in figure 1(b) where an interacting resonant level model defined in equation (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 [86]. 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 [86].

The discarded weight  = 10⁻⁶ 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⁻⁴ to 10⁻⁷. Moreover, the non-interacting case Δ = b = 0, λ = 1 allows for an exact evaluation of the steady-state current [76], which can be used to benchmark the accuracy of the DMRG data (a comparison can be found in [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 ${\mathrm { e}}^{-{\mathrm { i}}H\Delta t}{\mathrm { e}}^{{\mathrm { i}}\tilde H\Delta 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).

We now turn to current correlation functions in thermal equilibrium. Following equations (12) and (A.8), we need to evaluate

$$\begin{equation} {\rm e}^{-{\rm i}Ht}U_Q(t) j_{n} |\Psi_T\rangle. \end{equation} \tag{ 35 }$$

The ac conductivity σ(ω) is determined by the Fourier transform of 〈J(t)J〉 via equation (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 equation (13) is finite. D can be obtained from the long-time limit of the current–current correlation function via equation (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(a). At intermediate to high temperatures 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 [40, 44]. A fairly large  ∼ 10⁻⁵ 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For all temperatures from T = ∞ down to 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 [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).

We briefly present data for the spin–spin correlation function

$$\begin{equation} S^{zz}(n,t) = \langle S^z_{n+L/2}(t)S^z_{{L}/2}\rangle,\quad S^{zz}(k,t) = \sum_n {\rm e}^{{\rm i}kn} S^{zz}(n,t), \end{equation} \tag{ 36 }$$

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 0.1 (which for this situation corresponds to low temperature). A large discarded weight  ∼ 10⁻⁵ 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 [40] and also [33, 34] for transfer-matrix DMRG data). For finite Δ, the choice of $U_Q(t)={\mathrm { e}}^{{\mathrm { i}}\tilde Ht}$ still outperforms U_Q = 1 (except 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 [38] (see also [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.

Finally, we shortly investigate to what extend our choice of ${\mathrm { e}}^{{\mathrm { i}}\tilde Ht}$ for the time evolution of the auxiliaries is optimal (more details on the optimization of U_Q can be found in [31, 32]). To this end, we consider

$$\begin{equation} U_Q^\eta(t) = {\rm e}^{+{\rm i}\tilde Ht +\mathrm{ i}t\eta N}, \end{equation} \tag{ 37 }$$

with N begin an arbitrary Hermitian matrix. We compute the entanglement entropy

$$\begin{equation} S_{\rm ent}(t,\eta) =-2\sum_{a_n}(\Lambda^{L/2}_{a_n})^2\log_2\Lambda^{L/2}_{a_n} -2\sum_{a_{n_Q}}(\Lambda^{L/2_Q}_{a_{n_Q}})^2\log_2\Lambda^{L/2_Q}_{a_{n_Q}} \end{equation} \tag{ 38 }$$

between the left and right halves of the time-evolved state e^−iHtU^η_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

$$\begin{equation} N = 1 + 2 \sum_{n_Q=1}^{L-1} \big(S^x_{n_Q}S^x_{n_Q+1} + S^y_{n_Q}S^y_{n_Q+1}\big) \end{equation} \tag{ 39 }$$

as well as

$$\begin{equation} N = \sum_{n_Q=1}^{L-1} |\sigma_{n_Q}=\downarrow,\sigma_{n_Q+1}=\uparrow\rangle \langle\sigma_{n_Q}=\downarrow,\sigma_{n_Q+1}=\uparrow|. \end{equation} \tag{ 40 }$$

As mentioned above, Barthel et al and Barthel [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:

$$\begin{equation} \langle A(t)B\rangle = \langle A(t/2)B(-t/2)\rangle. \end{equation} \tag{ 41 }$$

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

$$\begin{equation} {\rm e}^{{\rm i}Ht/2}U_Q^\dagger(t/2) A \,{\rm e}^{-{\rm i}Ht/2}U_Q(t/2) |\Psi_T\rangle. \end{equation} \tag{ 42 }$$

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One way to evaluate equations (15) and (18) within the DMRG is to purify [41] the thermal density matrix by introducing an auxiliary Hilbert space Q:

$$\begin{equation} \rho_T = \frac{\mathrm{e}^{-H/T}}{Z_T} = {\rm Tr\,}_{Q} |\Psi_T\rangle\langle\Psi_T|. \end{equation} \tag{ A.1 }$$

At infinite temperature, where ρ_T is proportional to the unit operator, one can analytically write down the purification:

$$\begin{equation} \rho_\infty =\frac{1}{Z_\infty}=\frac{1}{2^{L}} = {\rm Tr}_{Q}|\Psi_\infty\rangle\langle\Psi_\infty| = {\rm Tr}_{Q}\prod_{n=1}^{L} |\Psi_\infty^{n}\rangle\langle\Psi_\infty^{n}|. \end{equation} \tag{ A.2 }$$

In order to exploit conservation laws within the DMRG algorithm, it is convenient to choose

$$\begin{equation} |\Psi_\infty^{n}\rangle = \frac{1}{\sqrt{2}}\left( \left|\sigma_n=\uparrow,\sigma_{n_Q}=\downarrow\right\rangle - \left|\sigma_n=\downarrow,\sigma_{n_Q}=\uparrow\right\rangle \right), \end{equation} \tag{ A.3 }$$

where σ_n = ↑,↓ denotes the eigenbasis of S^z_n, and Q is spanned by

$$\begin{equation} Q = {\rm span}\, \big\{ |\sigma_{1_Q}\ldots\sigma_{{L}_Q}\rangle \big\} . \end{equation} \tag{ A.4 }$$

One readily verifies that indeed

$$\begin{equation} \frac{1}{2} = \frac{1}{2}\sum_{\sigma_{n}} |\sigma_{n}\rangle\langle\sigma_{n}| = \sum_{\sigma_{n}}\sum_{\sigma_{n_Q}}\langle\sigma_{n_Q}|\Psi_\infty^{n}\rangle\langle\Psi_\infty^{n}|\sigma_{n_Q}\rangle\,. \end{equation} \tag{ A.5 }$$

The finite-temperature purified state |Ψ_T〉 is obtained from |Ψ_∞〉 by applying an imaginary time evolution [2],

$$\begin{eqnarray} \rho_T & =& \frac{\mathrm{e}^{-H/T}}{Z_T} = \frac{Z_\infty}{Z_T}{\rm Tr}_Q\, \mathrm{e}^{-H/2T} |\Psi_\infty\rangle\langle\Psi_\infty|\mathrm{e}^{-H/2T},\nonumber\\ & =& \frac{1}{\langle \Psi_T|\Psi_T\rangle}{\rm Tr}_Q |\Psi_T\rangle\langle\Psi_T|,~|\Psi_T\rangle=\mathrm{e}^{-H/2T}|\Psi_\infty\rangle , \end{eqnarray} \tag{ A.6 }$$

where we formally replaced

$$\begin{equation} H \to H\otimes 1_Q. \end{equation} \tag{ A.7 }$$

The correlation function of equation (15) is exactly recast as

$$\begin{equation} \langle A(t)B\rangle = \frac{\langle \Psi_\infty|\mathrm{ e}^{-\frac{H}{2T}}\,{\rm e}^{{\rm i}Ht} A {\rm e}^{-{\rm i}Ht} B\, \mathrm{e}^{-\frac{H}{2T}}|\Psi_\infty\rangle}{\langle\Psi_\infty|\mathrm{ e}^{-\frac{H}{T}}|\Psi_\infty\rangle}, \end{equation} \tag{ A.8 }$$

and for two time arguments one similarly finds

$$\begin{equation} \langle A(t)B(t')\rangle = \frac{\langle \Psi_T| \,{\rm e}^{{\rm i}Ht} A \,{\rm e}^{-{\rm i}Ht}\,{\rm e}^{{\rm i}Ht'} B\,{\rm e}^{-{\rm i}Ht'}|\Psi_T\rangle}{\langle\Psi_T|\Psi_T\rangle}. \end{equation} \tag{ A.9 }$$

In order to evaluate equation (7) for A = A^† = B, it is sufficient to compute

$$\begin{equation} {\rm e}^{{\rm i}Ht/2} A \,{\rm e}^{-{\rm i}Ht/2} |\Psi_T\rangle. \end{equation} \tag{ A.10 }$$

The expectation value of any observable A with respect to the time-dependent (non-equilibrium) density matrix of equation (18) is given by

$$\begin{equation} \langle A(t) \rangle = {\rm Tr}\,\rho(t) A = \frac{\langle \Psi_\infty|\mathrm{ e}^{-\frac{H_{\rm L}}{2T_{\rm L}}-\frac{H_{\rm R}}{2T_{\rm R}}}{\rm e}^{{\rm i}Ht} A\, {\rm e}^{-{\rm i}Ht}\,\mathrm{e}^{-\frac{H_{\rm L}}{2T_{\rm L}}-\frac{H_{\rm R}}{2T_{\rm R}}}|\Psi_\infty\rangle}{\langle\Psi_\infty|\mathrm{ e}^{-\frac{H_{\rm L}}{2T_{\rm L}}-\frac{H_{\rm R}}{2T_{\rm R}}}|\Psi_\infty\rangle}. \end{equation} \tag{ A.11 }$$

A.2.1. Initial state

It is convenient to first express the initial state |Ψ_∞〉 in terms of a matrix product state [5–8],

$$\begin{equation} |\Psi_\infty\rangle = \sum_{\sigma_n,\sigma_{n_Q}} A^{\sigma_1}A^{\sigma_{1_Q}}\cdots A^{\sigma_{{L}}}A^{\sigma_{{L}_Q}}|\sigma_1\sigma_{1_Q}\ldots\sigma_{{L}}\sigma_{{L}_Q}\rangle, \end{equation} \tag{ A.12 }$$

where

$$\begin{equation} A^{\sigma_i}_{a_ia_{i+1}} = \Lambda_{a_i}^i\Gamma_{a_ia_{i+1}}^{\sigma_i}. \end{equation} \tag{ A.13 }$$

Here and in the following σ_i is a short hand for either a physical or auxiliary degrees of freedom: σ_i∈{σ_n,σ_nQ}. The initial matrices corresponding to equation (A.3) read

$$\begin{eqnarray} \Gamma^{\sigma_{n_{\phantom{Q}}}=\uparrow}& = (1~0), \qquad\qquad \Gamma^{\sigma_{n_{\phantom{Q}}}=\downarrow} = (0~-1), \nonumber \\ \Gamma^{\sigma_{n_Q}=\uparrow}& = (0~1/\sqrt{2})^{\mathrm{T}}, \qquad \Gamma^{\sigma_{n_Q}=\downarrow} = (1/\sqrt{2}~0)^{\mathrm{T}} , \end{eqnarray} \tag{ A.14 }$$

as well as Λⁱ = 1. Normalization generally stipulates

$$\begin{equation} \sum_{\sigma_ia_i} (\Lambda^i_{a_i}\Gamma^{\sigma_i}_{a_ia_{i+1}})^\dagger (\Lambda_{a_i}^i\Gamma^{\sigma_i}_{a_ia_{i+1}'}) = \delta_{a_{i+1}a_{i+1}'} \end{equation} \tag{ A.15 }$$

as well as

$$\begin{equation} \sum_{\sigma_ia_{i+1}} (\Gamma^{\sigma_i}_{a_ia_{i+1}}\Lambda^{i+1}_{a_{i+1}})(\Gamma^{\sigma_i}_{a_i'a_{i+1}}\Lambda_{a_{i+1}}^{i+1})^\dagger = \delta_{a_{i}a_{i}'}. \end{equation} \tag{ A.16 }$$

A.2.2. Trotter decomposition

In order to successively apply the imaginary and real time evolutions operators

$$\begin{equation} \mathrm{e}^{-\lambda H} = \mathrm{e}^{-\Delta\lambda H}\,\mathrm{e}^{-\Delta\lambda H}\ldots ,~~\lambda = \Delta\lambda+ \Delta\lambda + \cdots \end{equation} \tag{ A.17 }$$

to the initial state |Ψ_∞〉, we factorize them either via a second or a fourth order Trotter decomposition [2]. The second order scheme reads

$$\begin{equation} \mathrm{e}^{-\Delta\lambda H} \approx \mathrm{e}^{-\Delta\lambda/2 H_1} \,\mathrm{e}^{-\Delta\lambda H_2} \,\mathrm{e}^{-\Delta\lambda/2 H_1}, \end{equation} \tag{ A.18 }$$

where the H = H₁ + H₂, and H_1,2 contain only terms which mutually commute:

$$\begin{equation} H_1 = \sum_{n{\rm odd}} h_n,~H_2 = \sum_{n{\rm even}} h_n. \end{equation} \tag{ A.19 }$$

In order to reduce the error, one can employ a fourth order decomposition [83],

$$\begin{equation} \mathrm{e}^{-\Delta\lambda H} \approx U(\Delta\lambda_1)^2\, U(\Delta\lambda_2)\, U(\Delta\lambda_1)^2, \end{equation} \tag{ A.20 }$$

where

$$\begin{equation} U(\Delta\lambda_i) = \mathrm{e}^{-\Delta\lambda_i H_1/2}\,\mathrm{e}^{-\Delta\lambda_i H_2}\,\mathrm{e}^{-\Delta\lambda_i H_1/2} \end{equation} \tag{ A.21 }$$

and

$$\begin{equation} \Delta\lambda_1=\frac{1}{4-4^{1/3}}\Delta\lambda,\quad \Delta\lambda_2=\Delta\lambda-4\Delta\lambda_1. \end{equation} \tag{ A.22 }$$

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].

A.2.3. DMRG step

The physical Hamiltonians of equations (8) and (10) couple matrices with indices σ_n and σ_n+1 which in the matrix product state of equation (A.12) are next-nearest neighbors—they are separated by a matrix associated with an auxiliary degree of freedom σ_nQ. Through equation (A.7), the local Hamiltonians of equations (8) and (10) are formally replaced by

$$\begin{equation} h_n \to h_n\otimes \sum_{\sigma_{n_Q}}|\sigma_{n_Q}\rangle\langle\sigma_{n_Q}|. \end{equation} \tag{ A.23 }$$

Thus, after each application of $\exp (-\Delta \lambda 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:

$$\begin{eqnarray} \Gamma^{\sigma_{0} } & \longrightarrow \tilde \Gamma^{\sigma_{0}}, \nonumber \\ \Gamma^{\sigma_{1}} & \longrightarrow \tilde \Gamma^{\sigma_{1}}, \quad \Lambda^{1}\longrightarrow \tilde \Lambda^{1},\nonumber \\ \Gamma^{\sigma_{2}} & \longrightarrow \tilde \Gamma^{\sigma_{2}}, \quad \Lambda^{2}\longrightarrow \tilde \Lambda^{2}, \end{eqnarray} \tag{ A.24 }$$

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 [2]. Firstly, one forms the three-site tensor

$$\begin{equation} \Psi_{a_0a_3}^{\sigma_0\sigma_1\sigma_2} = \sum_{a_1a_2} \Lambda_{a_0}^0\Gamma^{\sigma_0}_{a_0a_1} \Lambda_{a_1}^1\Gamma^{\sigma_1}_{a_1a_2} \Lambda_{a_2}^2\Gamma^{\sigma_2}_{a_2a_3} \Lambda_{a_3}^3, \end{equation} \tag{ A.25 }$$

which is then acted on by $\exp (-\Delta \lambda h)$:

$$\begin{equation} \Phi_{a_0a_3}^{\sigma_0\sigma_1\sigma_2} = \sum_{\sigma_1'\sigma_2'\sigma_3'} \Psi_{a_0a_3}^{\sigma_0'\sigma_1'\sigma_2'} \,\mathrm{e}^{-\Delta\lambda h}\big|_{\sigma_0\sigma_1\sigma_2}^{\sigma_0'\sigma_1'\sigma_2'}. \end{equation} \tag{ A.26 }$$

Now, two singular value decompositions (SVDs) are applied to appropriately reshaped tensors in order to obtain the new matrices $\tilde \lambda ^i$ and $\tilde \Gamma ^{\sigma _i}$:

$$\begin{eqnarray}&&\begin{array}{@{}l@{}} \Phi_{a_0a_3}^{\sigma_0\sigma_1\sigma_2} = \Phi_{(a_0\sigma_0),(a_3\sigma_1\sigma_2)},\\ = \sum_{a_1} U_{(a_0\sigma_0),a_1} \tilde\Lambda^1_{a_1} (V^\dagger)_{a_1,(\sigma_1\sigma_2a_3)},\\ = \sum_{a_1}\Lambda^0_{a_0}\underbrace{\Lambda^{0,-1}_{a_0}U^{\sigma_0}_{a_0a_1}}_{\tilde\Gamma^{\sigma_0}_{a_0a_1}} \underbrace{\tilde\Lambda^1_{a_1} (V^\dagger)_{a_1,(\sigma_1\sigma_2a_3)}}_{\tilde\Phi_{a_1,(\sigma_1\sigma_2a_3)}},\end{array} \end{eqnarray} \tag{ A.27 }$$

and likewise

$$\begin{eqnarray}&&\begin{array}{@{}l@{}} \tilde\Phi_{a_1,(\sigma_1\sigma_2a_3)} = \tilde\Phi_{(a_1\sigma_1),(\sigma_2a_3)}, \\ = \sum_{a_2} U_{(a_1\sigma_1),a_2}\tilde\Lambda_{a_2}^2(\tilde V^\dagger)_{a_2,(\sigma_2a_3)}, \\ = \sum_{a_2} \tilde\Lambda_{a_1}^1 \underbrace{\tilde\Lambda^{1,-1}_{a_1}U^{\sigma_1}_{a_1a_2}}_{\tilde\Gamma^{\sigma_1}_{a_1a_2}} \tilde\Lambda_{a_2}^2 \underbrace{(\tilde V^\dagger)_{a_2a_3}^{\sigma_2} \Lambda_{a_3}^{3,-1}}_{\tilde\Gamma^{\sigma_2}_{a_2a_3}}\Lambda_{a_3}^3.\end{array} \end{eqnarray} \tag{ A.28 }$$

In case of a real time evolution, the normalization conditions of equations (A.15) and (A.16) are preserved automatically (up to errors associated with a potential truncation of the matrices; see below). For a non-unitary $\exp (-\Delta \beta 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 [82].

A.2.4. The discarded weight