Precise evaluation of thermal response functions by optimized density matrix renormalization group schemes

This paper addresses the efficient evaluation of finite-temperature response functions

$$\begin{equation} \chi_{\skew3\hat{A}\skew3\hat{B}}(\beta,t):=\frac{1}{Z_\beta}{\rm Tr}(\mathrm{e}^{-\beta \skew5\hat{H}}\skew3\hat{B}(t)\skew3\hat{A})\quad\textrm{with}\quad \skew3\hat{B}(t)\equiv {\rm e}^{{\rm i}\skew5\hat{H} t}\skew3\hat{B} \,{\rm e}^{-{\rm i}\skew5\hat{H} t} \end{equation} \tag{ 1 }$$

for quantum many-particle systems on one-dimensional (1D) lattices in the framework of the density matrix renormalization group (DMRG) [1–3]. From the theoretical perspective, such quantities occur for example in the context of linear response theory [4, 5] and characterize the effect of a perturbation $\hat {A}$ of the system at time zero on the expectation value of an observable $\hat {B}$ at time t. Initially, the system is in thermal equilibrium, where β = 1/T is the inverse temperature and $Z_\beta ={\mathrm { Tr}} \,\mathrm {e}^{-\beta \skew5\hat{H}}$. The units are chosen such that Boltzmann's and Planck's constants are k_B = 1 and ℏ = 1. Such response functions contain important information on the many-body physics and are addressed in many experimental setups. For example, recent advances in neutron-scattering techniques make very precise measurements possible. See for example [6–10]. It is hence very important to have numerical tools at hand that allow for an efficient and precise evaluation of thermal response functions for (strongly correlated) condensed matter models in order to match theoretical models to actual materials and to gain an understanding of the underlying physical processes.

As discussed and demonstrated in several studies [11–13], finite-temperature response functions for strongly correlated 1D systems can be evaluated up to some maximum reachable time by using the time-dependent DMRG (tDMRG) [14–16] applied to a purification [17–19] of the density matrix. In the DMRG approach, those purifications are approximated with a controllable precision by matrix product states (MPS) [20–22]. To address finite-temperature states with tDMRG was suggested in [23, 24] with first applications in [25, 26]. A difficulty in the simulations of time-evolved states is the growth of entanglement with time [27–30]. In tDMRG calculations, this leads to a corresponding, typically exponential, increase of the computation cost with time and a strong limitation of the maximum reachable times, depending on the available computational resources and the desired accuracy of the simulation. The effect is much more drastic for mixed states. The focus of this paper is hence to study the computation cost of the different tDMRG schemes for the evaluation of thermal response functions and to suggest novel more efficient approaches that allow for a substantial increase in the maximum reachable times.

The results of such simulations can be Fourier transformed to study the spectral properties of the response. In order to avoid ringing artifacts from the Fourier transformation of the data on a restricted time interval, one can either use filters, which result in an artificial broadening, or use linear prediction [31, 32]. This was first employed in [33] for T = 0 and in [11] for T > 0. See also [34, 35] for further applications of linear prediction in DMRG calculations at T = 0.

Earlier alternative DMRG approaches for finite-temperature response functions have built on the transfer matrix renormalization group (TMRG) [36–39]. In a first step, analytic continuation techniques as in quantum Monte Carlo (QMC) were employed [40]. In a second development, transfer matrices comprising the imaginary- and the real-time evolution were used to access autocorrelation functions [41, 42]. In finite-temperature QMC calculations (e.g. positive-definite path integrals [43, 44] or stochastic series expansion with directed loops [45]), real-time or real-frequency response functions have to be extracted by analytic continuation from imaginary-time results [46], which is ill-conditioned and numerically challenging. Finite temperatures have also been addressed by a hybrid algorithm based on Monte Carlo and tDMRG [47, 48].

In this paper, I discuss how the simulation based on the evolution of a certain MPS purification of the thermal density matrix as described in [11, 13, 23] can actually be understood as evolving a matrix product operator (MPO) representation of the square root of the thermal density matrix. The description in terms of purifications is in this sense superfluous. A decisive observation is now that there are considerable degrees of freedom for the choice of a specific evaluation scheme. For the two specific schemes, corresponding to [11, 12] (scheme A) and [13] (scheme B), respectively, I examine the scaling of the computation cost with time. It turns out that scheme A has some advantage at low temperatures and scheme B is advantageous at higher temperatures, for which I give an explanation on the basis of quasi-locality [49, 50]. The maximum reachable times t^A_max(β) and t^B_max(β) differ, for the same computation cost and accuracy, by a factor of order 1. It is then shown that an optimized evaluation scheme, for which one optimizes over the aforementioned degrees of freedom, yields maximum reachable times that are at least twice as large as for scheme B, specifically t^opt_max(2β) ≳ 2t^B_max(β). In all cases one finds that t_max(2β) ≈ t_max(β) as t_max varies slowly as a function of log β. We can devise a new scheme C that does not require any optimization, but outperforms scheme B by a factor of 2 in the maximum reachable time. It is only outdone by scheme A at very low temperatures. To go beyond scheme C, one can study the computation cost for calculating MPO approximations of operators ${\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} s}\,\mathrm {e}^{-\alpha \skew5\hat{H}}\hat {B} \,{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} s'}$ to determine optimized evolution schemes. Finally, the connection between scheme A and an earlier TMRG approach [41] is explained. For each of the tDMRG schemes, one can devise a corresponding TMRG variant.

Due to the substantially increased reachable times and the simplicity (no need for optimization), the novel scheme C, which was introduced recently in [51] and is explained here in detail, can be expected to be the method of choice for future applications. In [51], it allowed us to show that the thermal spectral functions of 1D bosons in the quantum critical regime with dynamic critical exponent z = 2 follow a universal scaling form and to compute the corresponding scaling function precisely.

As described for example in [11, 13, 23], the (unnormalized) thermal density matrix

$$\begin{equation*} \hat\rho_\beta:=\mathrm{e}^{-\beta \skew5\hat{H}},\quad Z_\beta:={\rm Tr} \,\mathrm{e}^{-\beta \skew5\hat{H}} \end{equation*}$$

for a chain of L sites can be obtained in the form of an MPS purification

$$\begin{equation} |\rho_\beta\rangle \in {\mathcal{H}}\otimes{\mathcal{H}}_{{\rm aux}},\quad {\mathcal{H}}_{{\rm aux}}\simeq {\mathcal{H}}, \end{equation} \tag{ 2a }$$

$$\begin{equation} |\rho_\beta\rangle = \sum_{{\boldsymbol{\sigma}},{\boldsymbol{\sigma}}'} A^{\sigma_1, \sigma'_1}_1 A^{\sigma_2, \sigma'_2}_2\ldots A^{\sigma_L, \sigma'_L}_L \underbrace{|\sigma_1\sigma_2\ldots \sigma_L\rangle}_{=:|{\boldsymbol{\sigma}}\rangle}\otimes \underbrace{|\sigma'_1 \sigma'_2\ldots \sigma'_L\rangle_{{\rm aux}}}_{=:|{\boldsymbol{\sigma}}'\rangle_{{\rm aux}}}, \end{equation} \tag{ 2b }$$

${\mathcal {H}}$

${\mathcal {H}}_{{\mathrm { aux}}}$

${\mathcal {H}}$

${\mathcal {H}}_{{\mathrm { aux}}}$

In order for |ρ_β〉 to be a purification of $\hat \rho _\beta$, it needs to be constructed such that

$$\begin{equation*} {\rm Tr}_{{\rm aux}}|\rho_\beta\rangle\langle\rho_\beta|\equiv \sum_{{\boldsymbol{\sigma}}{\boldsymbol{\sigma}}'{\boldsymbol{\sigma}}''}|{\boldsymbol{\sigma}}\rangle\langle {\boldsymbol{\sigma}}{\boldsymbol{\sigma}}'|\rho_\beta\rangle\langle\rho_\beta|{\boldsymbol{\sigma}}''{\boldsymbol{\sigma}}'\rangle\langle {\boldsymbol{\sigma}}''| =\hat\rho_\beta. \end{equation*}$$

This can be achieved by choosing the infinite-temperature purification to be

$$\begin{equation*} |\rho_0\rangle=|\mathbbm{1}\rangle\quad\textrm{with}\quad |\mathbbm{1}\rangle := \sum_{{\boldsymbol{\sigma}}}|{\boldsymbol{\sigma}}\rangle\otimes|{\boldsymbol{\sigma}}\rangle_{{\rm aux}} = \bigotimes_{i=1}^L\left(\sum_{\sigma_i} |\sigma_i\rangle\otimes|\sigma_i\rangle_{{\rm aux}}\right) \end{equation*}$$

which can be written as an MPS (2) with bond dimensions M_i = 1. In this so-called ancilla approach, finite-temperature purifications can be calculated by imaginary-time evolution,

$$\begin{equation*} |\rho_\beta\rangle = {\rm e}^{-\beta \skew5\hat{H}/2}\otimes \mathbbm{1}_{{\rm aux}}\cdot |\mathbbm{1}\rangle , \end{equation*}$$

and response functions are evaluated after a subsequent real-time evolution,

$$\begin{equation} \chi_{\skew3\hat{A}\skew3\hat{B}}(\beta,t) = \frac{1}{Z_\beta}\big[\langle\rho_\beta| {\rm e}^{{\rm i}\skew3\hat{H} t}\big]\skew3\hat{B} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} |\rho_\beta\rangle\big].\end{equation} \tag{ 3 }$$

The evolution of the MPS can for example be implemented by decomposing the propagators into circuits of local gates by a Trotter–Suzuki decomposition [14–16] as described in section 6. The square brackets in equation (3) indicate which parts of the expression are represented as MPS after the evolution. Note that Z_β, which is required in equation (3), can be determined as $\langle \rho _\beta |\rho _\beta \rangle ={\mathrm { Tr}}\;\mathrm {e}^{-\beta \skew5\hat{H}} =Z_\beta$. Equivalently, one can work with purifications |ρ_β〉 that are normalized to one.

In each step of the tDMRG, the evolved states |ψ〉 = |ψ(β,t)〉 are approximated by an MPS with bond dimensions M_i = M_i(β,t) that are as small as possible for a given constraint on the desired precision of the approximation [3]. This is achieved through truncations. For every splitting of the system into a left and a right part, one does a Schmidt decomposition $|\psi \rangle =\sum _{k=1}^{{\skew3\tilde{M}}}\lambda _k|k\rangle _{\mathrm {L}}\otimes |k\rangle _{\mathrm {R}}$ of the state [19], which boils down to doing singular value decompositions of the tensors A_i. The corresponding reduced density matrices are $\sum _k\lambda _k^2|k\rangle _{\mathrm {L}}\langle k|_{\mathrm {L}}$ and $\sum _k\lambda _k^2|k\rangle _{\mathrm {R}}\langle k|_{\mathrm {R}}$. The bond dimension is then reduced from ${\skew3\tilde{M}}$ to some value $M{<}{\skew3\tilde{M}}$ by retaining only the M largest Schmidt coefficients and truncating all smaller ones:

$$\begin{equation} |\psi\rangle=\sum_{k=1}^{{\skew3\tilde{M}}}\lambda_k|k\rangle_{\mathrm{L}}\otimes|k\rangle_{\mathrm{R}}\mapsto |\psi_{{\rm trunc}}\rangle=\sum_{k=1}^{M}\lambda_k|k\rangle_{\mathrm{L}}\otimes|k\rangle_{\mathrm{R}}. \end{equation} \tag{ 4 }$$

The precision is in each step of the algorithm controlled by bounding the truncation weight

$$\begin{equation} \epsilon=\left(\frac{\Vert \psi_{{\rm trunc}} - \psi\Vert}{\Vert \psi\Vert}\right)^2 = \frac{\sum_{k>M}\lambda_k^2}{\sum_k\lambda_k^2}. \end{equation} \tag{ 5 }$$

As a matter of fact, the notation of this procedure as an algorithm on purifications is in a sense superfluous. Due to the fact that $\mathcal {B}({\mathcal {H}})$, the space of linear maps on ${\mathcal {H}}$, and the tensor product space ${\mathcal {H}}\otimes {\mathcal {H}}$ are isomorphic, everything can be rewritten in terms of MPOs in $\mathcal {B}({\mathcal {H}})$ instead of MPS purifications in ${\mathcal {H}}\otimes {\mathcal {H}}$:

$$\begin{equation} \langle{\boldsymbol{\sigma}}{\boldsymbol{\sigma}}'|X\rangle \equiv \langle{\boldsymbol{\sigma}}|\hat X|{\boldsymbol{\sigma}}'\rangle\quad\forall{{\boldsymbol{\sigma}}{\boldsymbol{\sigma}}'}. \end{equation} \tag{ 6 }$$

Examples for this one-to-one relation are

$$\begin{eqnarray*} &&\hspace{44pt}|\mathbbm{1}\rangle \longleftrightarrow \mathbbm{1},\\ &&|\rho_\beta\rangle \longleftrightarrow \mathrm{e}^{-\beta\skew5\hat{H}/2},\\ &&\hspace{44pt}\hat Y\otimes\hat Z_{{\rm aux}} |X\rangle \longleftrightarrow \hat Y\hat X \hat Z^{\mathrm{T}}, \end{eqnarray*}$$

where the transposition is to be executed in the $|{\boldsymbol{\sigma}}\rangle$ basis. In particular, the counterpart of the MPS purification (2) is the MPO

$$\begin{equation} \sum_{{\boldsymbol{\sigma}}{\boldsymbol{\sigma}}'} A^{\sigma_1, \sigma'_1}_1 A^{\sigma_2, \sigma'_2}_2\ldots A^{\sigma_{\mathrm{L}}, \sigma'_{\mathrm{L}}}_{\mathrm{L}} |{\boldsymbol{\sigma}}\rangle\langle{\boldsymbol{\sigma}}'|. \end{equation} \tag{ 7 }$$

See for example [24, 52, 53] for discussions of MPOs.

Applying evolution operators $\mathrm {e}^{-\Delta \beta \skew5\hat{H}}$ or $\mathrm {e}^{\pm \mathrm {i}\Delta t \skew5\hat{H}}$ to the left or right of an MPO $\hat X$ can be done by tDMRG in exactly the same way as for MPS, for example, by a Trotter–Suzuki decomposition as described in section 6. The truncation weight, which is kept below a certain bound in order to control the precision, is then given by

$$\begin{equation} \epsilon=\left(\frac{\Vert \hat X_{{\rm trunc}} - \hat X\Vert_2}{\Vert \hat X\Vert_2}\right)^2, \end{equation} \tag{ 8 }$$

where $\Vert \hat X\Vert _{2}=\sqrt {{\mathrm { Tr}}\,\hat X^{\dagger} \hat X}$ is the Schatten 2-norm. This is the natural norm within the DMRG framework.

Evaluation scheme A as described and employed in [11, 12] corresponds to equation (3). In the MPO notation it reads

$$\begin{equation} \chi^A_{\skew3\hat{A}\skew3\hat{B}}(\beta,t) = \frac{1}{Z_\beta}{\rm Tr}\left(\big[\mathrm{e}^{-\beta\skew5\hat{H}/2} \,{\rm e}^{{\rm i}\skew5\hat{H} t}\big]\skew3\hat{B} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\big]\right). \end{equation} \tag{ 9 }$$

The square brackets indicate which parts of the expression are represented as MPOs after the evolution. Figure 2 shows the algorithm diagrammatically. The partition function Zβ can be obtained from the Schatten 2-norm $Z_\beta =(\Vert [\mathrm {e}^{-\beta \skew5\hat{H}/2}]\Vert _{2})^2$. In the DMRG method, multiplications of the matrices Aσi,σ'ii and singular value decompositions, required for the truncations of the MPOs, dominate the computation costs. Hence, the computation costs for a single evolution step for the first MPO (bond dimensions Mi) and the second MPO (bond dimensions Mi'), and for the final evaluation of equation (9) scale as

$$\begin{equation} \sum_{i=1}^L M_i^3,\quad \sum_{i=1}^L (M'_i)^3\quad\textrm{and}\quad \sum_{i=1}^L {{\rm max}}\left(M_i(M'_i)^2,M_i^2 M'_i\right), \end{equation} \tag{ 10 }$$

respectively²; see figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For a non-singular operator $\hat T\in \mathcal {B}({\mathcal {H}})$, the value of (9) is invariant under transformations

$$\begin{equation} \big[\mathrm{e}^{-\beta\skew5\hat{H}/2} \,{\rm e}^{{\rm i}\skew5\hat{H} t}\big] \quad\longrightarrow\quad \big[\hat T \,\mathrm{e}^{-\beta\skew5\hat{H}/2} \,{\rm e}^{{\rm i}\skew5\hat{H} t}\big], \end{equation} \tag{ 11a }$$

$$\begin{equation} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\big] \quad\longrightarrow\quad \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\hat T^{-1}\big] \end{equation} \tag{ 11b }$$

$\hat T$

$\hat T$

$\hat T=\mathrm {e}^{-\beta '\skew5\hat{H}}{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t'}$

The modified evaluation scheme B (see figure 2) employed in [13, 56] corresponds to the choice $\hat T={\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}$ and reads in the operator notation

$$\begin{eqnarray} \chi^B_{\skew3\hat{A}\skew3\hat{B}}(\beta,t) &=& \frac{1}{Z_\beta}{\rm Tr}\left(\big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\,\mathrm{e}^{-\beta\skew5\hat{H}/2} \,{\rm e}^{{\rm i}\skew5\hat{H} t}\big]\skew3\hat{B} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\,{\rm e}^{{\rm i}\skew5\hat{H} t}\big]\right)\nonumber\\\ms\ms &=& \frac{1}{Z_\beta}{\rm Tr}\left(\big[\mathrm{e}^{-\beta\skew5\hat{H}/2}\big]\skew3\hat{B} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\,{\rm e}^{{\rm i}\skew5\hat{H} t}\big]\right). \end{eqnarray} \tag{ 12 }$$

In order to study the computational costs of evaluation schemes A and B (equations (9) and (12)), as functions of time and temperature, let us choose as an example the spin-1/2 XXZ Heisenberg model

$$\begin{equation} \hat H = \sum_{i=1}^{L-1}(\hat{S}^x_i\hat{S}^x_{i+1}+\hat{S}^y_i\hat{S}^y_{i+1}+J_z \hat{S}^z_i\hat{S}^z_{i+1}) -h\sum_{i=1}^L\hat{S}^z_i. \end{equation} \tag{ 13 }$$

The phase diagram [57–59], as derived by the Bethe ansatz, is shown in figure 3. The simulations are carried out for vanishing magnetic field h = 0, system size L = 128 and coupling constants J_z = 0 and 1, where the model is gapless, and at J_z = 3, where the model is in its gapped antiferromagnetic phase. For the time evolution, a fourth-order Trotter–Suzuki decomposition with step sizes Δβ = Δt = 1/8 is used. The truncation weights in the imaginary-time and real-time evolutions were fixed to values _β and _t, which are specified in the captions of the corresponding figures.

Zoom In Zoom Out Reset image size

Figure 4(a) shows the computation cost per time step for the operators $\skew3\hat {A}=\skew3\hat {B}^{\dagger} =\hat {S}^+_{L/2}$ in equations (9) and (12), i.e. spin flips in the middle of the system⁴. The density plots show, as a function of time and temperature, the number of operations needed for a single time step as measured by $\sum _i \left (M_i(\beta ,t)\right )^3$. For scheme A, the cost is dominated by the cost for the computation of the MPO

$$\begin{equation} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\big]. \end{equation} \tag{ 14 }$$

For scheme B, the computation of the MPO

$$\begin{equation} \big[{\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\,{\rm e}^{{\rm i}\skew5\hat{H} t}\big] \end{equation} \tag{ 15 }$$

dominates the numerical costs. Hence, the bond dimensions M_i(β,t) of those operators were used in the diagrams. Note the logarithmic scale for β = 1/T and the cost, i.e. the color coding and the distances of the contour lines (see footnote ²). The contour lines correspond to maximum reachable times for certain computational resources and a predefined precision of the simulation, determined by _β and _t. Figure 4(b) shows, for the same simulations, the maximum bond dimensions max_i∈[1,L]M_i(β,t).

Zoom In Zoom Out Reset image size

First of all, the results show that the computation costs increase, for fixed temperature, exponentially with time t. This corresponds to a linear increase of the entanglement entropy in the evolution of pure states after a quench [27–29]. For the non-critical case, J_z = 3, the costs per time step become β-independent for temperatures that are sufficiently below the excitation gap. At higher temperatures, the cost for scheme B is systematically smaller than that for scheme A. The effect is strongest for the non-interacting system at J_z = 0. For J_z = 1 and 3, the increase in the maximum reachable times is more moderate with a factor of ≲1.4. At lower temperatures, scheme A is more efficient than scheme B. This trend strengthens when the accuracy of the simulation is reduced (higher truncation weights _β and _t), as documented in the left part of figure 5. Figure 4(b) shows that the maximum bond dimensions max_i M_i(β,t) evolve almost in the same way as the computation costs. For the non-critical case J_z = 3, one sees however that the maximum bond dimensions occurring in scheme A are for all temperatures smaller than those occurring in scheme B, whereas the computation cost for scheme B is lower at higher temperatures. Also, for J_z = 1, the maximum bond dimensions occurring in the two schemes are quite similar.

Zoom In Zoom Out Reset image size

These findings can be explained as follows. In [13] it was pointed out that the operator (15) of scheme B is time independent for the simple case $\hat {A}=\mathbbm {1}$, whereas the computation cost for the MPO (14), occurring in scheme A, can increase with time, even in this trivial case. More generally, the following argument applies for all operators $\skew3\hat {A}$ with finite spatial support ${{\mathrm { supp}}}(\hat {A})$. The typical condensed matter systems are quasi-local [49, 50], i.e. the spatial support of operators such as ${\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A} \,{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}$, occurring in scheme B, grows only linearly with time. More precisely, outside a certain space–time cone originating from ${{\mathrm { supp}}}(\hat {A})$, the evolved operator acts almost like the identity and hence does not change the entanglement in that region; see figure 6. In mathematical terms, quasi-locality means

$$\begin{equation} \left\Vert {\rm e}^{-{\rm i}\skew5\hat{H} t}\skew3\hat{A} \,{\rm e}^{{\rm i}\skew5\hat{H} t}-{\rm e}^{-{\rm i}\skew5\hat{H}_V t}\skew3\hat{A} \,{\rm e}^{{\rm i}\skew5\hat{H}_V t}\right\Vert \leqslant C\Vert \skew3\hat{A}\Vert\mathrm{ e}^{vt-{\rm dist}\left(\partial V,{{\rm supp}}(\skew3\hat{A})\right)}, \end{equation} \tag{ 16 }$$

where all terms $\hat h_i$ in the Hamiltonian $\skew5\hat{H}=\sum _i\hat h_i$ are required to be short range and norm bounded, $\skew5\hat{H}_V=\sum _{i\in V}\hat h_i$ is the Hamiltonian truncated to a vicinity V of the spatial support of the operator $\hat {A}$ (${{\mathrm {supp}}}(\hat {A})\subset V$), ∂V is the boundary of V and $v\propto {\mathrm { sup}}_i \Vert \hat h_i\Vert$ is the Lieb–Robinson velocity [49, 50]. So the norm distance between the exactly evolved operator and the operator evolved on the subsystem V only, decays exponentially with increasing distance of ∂V to a space–time cone defined by the Lieb–Robinson velocity v, and originating from ${{\mathrm { supp}}}(\hat {A})$ at time zero. As the evolved operator ${\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A} \,{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}$, hence, behaves up to exponentially small corrections like the identity outside the specified space–time cone, it does not alter the operator $[\mathrm {e}^{-\beta \skew5\hat{H}/2}]$ in that region—in particular, not the MPO bond dimensions M_i—for the product (15) occurring in scheme B. This effect is very visible in the plots for J_z = 1 and β = 3 shown in figure 6. For both tDMRG schemes, the maximum bond dimensions occur in the middle of the system, where $\hat {A}=\hat {S}^+_{L/2}$ acts, and grow (exponentially) with time t. The maxima have comparable values. Whereas in scheme B the M_i remain unchanged outside the Lieb–Robinson space–time cone due to the quasi-locality, they grow for all bonds in scheme A. This implies, according to equation (10), an increased cost for scheme A at those temperatures.

Zoom In Zoom Out Reset image size

Nevertheless, as figures 4–6 exemplify, scheme A is often advantageous at low temperatures, especially for non-critical systems. For low temperatures, the quasi-locality is not essential as $\mathrm {e}^{-\beta \skew5\hat{H}/2}$ limits the effect of the real-time propagators $\mathrm {e}^{\pm \mathrm {i}\skew5\hat{H} t}$ to a low-energy subspace. Sufficiently far away from the support of $\hat {A}$, the resulting dynamics is then very restricted and cannot lead to a large increase of the bond dimensions. This is true for both schemes. Furthermore, in the middle of the system, where $\hat {A}=\hat {S}^+_{L/2}$ acts, the growth of the bond dimensions is slower for scheme A, which can be attributed to the fact that one acts with one propagator in equation (14) instead of acting with two propagators in equation (15). As the plots for J_z = 1,3 and β = 40,48 in figure 6 exemplify, the effect is much more pronounced for gapped systems.

The described properties have also been found for other local operators $\skew3\hat {A}$ and $\skew3\hat {B}$ such as $\hat {S}^z_{L/2}$, etc. The behavior for non-local operators such as $\hat {S}^+_{k=0}:=\sum _i\hat {S}^+_i$ is also quite similar, as exemplified in figure 5 for J_z = 0,3. Note that $[\hat {S}^+_{k=0},\skew5\hat{H}]=0$ in the isotropic case J_z = 1, for which the evolution of $\hat {S}^+_{k=0}(t)$ is hence trivial.

Schemes A and B, equations (9) and (12), are typically far from optimal. One has a lot of freedom in designing a scheme that is as efficient as possible. The choice for a specific measure for the efficiency of a scheme can depend on the available computational resources. In this paper, the computation cost per time step as quantified by $\sum _i \left (M_i(\beta ,t)\right )^3$ (equation (10)) and, as an alternative, the maximum bond dimension max_i M_i(β,t) are considered as useful measures. With a non-singular operator $\hat T$, the most generic scheme involving two MPOs is $Z_\beta ^{-1}{\mathrm { Tr}}([{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}\hat {B}\hat T] [\hat T^{-1}\,{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A}\, \mathrm {e}^{-\beta \skew5\hat{H}}])$. As pointed out in section 3, to optimize efficiently with respect to generic $\hat T$ is, in general, not possible. Thus, it is reasonable to constrain ourselves to a certain subclass of schemes and to optimize the efficiency over the remaining degrees of freedom of that class. Let us consider the class of schemes

$$\begin{equation} \chi^{\beta't't''}_{\skew3\hat{A}\skew3\hat{B}}(\beta,t)= \frac{1}{Z_\beta}{\rm Tr}\left(\big[{\rm e}^{{\rm i}\skew5\hat{H} t'}\,\mathrm{e}^{-\beta'\skew5\hat{H}}\skew3\hat{B} \,{\rm e}^{-{\rm i}\skew5\hat{H} t''}\big] \big[{\rm e}^{-{\rm i}\skew5\hat{H} (t-t'')}\skew3\hat{A} \,\mathrm{e}^{-(\beta-\beta')\skew5\hat{H}}\,{\rm e}^{{\rm i}\skew5\hat{H}(t-t')}\big]\right).\end{equation} \tag{ 17 }$$

As before, the two factors in square brackets are to be approximated by MPOs (see footnote ³). For given inverse temperature β, time t and accuracy as quantified by the truncation weight (equation (8)), one can now optimize the chosen efficiency measure with respect to β', t' and t''. Scheme A corresponds to the choice (β',t',t'') = (β/2,t,0) and scheme B to (β',t',t'') = (β/2,0,0). In actual applications, one can first study the computation cost, for obtaining operators $[{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} s}\,\mathrm {e}^{-\alpha \skew5\hat{H}}\hat {B} \,{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} s'}]$ and $[{\mathrm { e}}^{-{\mathrm { i}}\skew3\hat{H} s}\hat {A} \,\mathrm {e}^{-\alpha \skew3\hat{H}}\,{\mathrm { e}}^{{\mathrm { i}}\skew3\hat{H} s'}]$. From the results, one can then conclude on optimal values for β', t' and t'' in the evaluation of the response function according to equation (17). For the demonstrational purposes of this paper, it is sufficient to not explore the full three-dimensional parameter space, but to constrain ourselves in the following to the subspace with t' = t'', i.e. to the schemes

$$\begin{eqnarray} \chi^{\beta't'}_{\skew3\hat{A}\skew3\hat{B}}(\beta,t) &= \frac{1}{Z_\beta}{\rm Tr}\left(\big[\mathrm{e}^{-\beta'\skew5\hat{H}}\skew3\hat{B}(t')\big] \big[\skew3\hat{A}(t'-t)\, \mathrm{e}^{-(\beta-\beta')\skew5\hat{H}}\big]\right)\nonumber\\ &= \frac{1}{Z_\beta}{\rm Tr}\left(\big[{\rm e}^{{\rm i}\skew5\hat{H} t'}\,\mathrm{e}^{-\beta'\skew5\hat{H}}\skew3\hat{B} \,{\rm e}^{-{\rm i}\skew5\hat{H} t'}\big] \big[{\rm e}^{-{\rm i}\skew5\hat{H} (t-t')}\skew3\hat{A} \,\mathrm{e}^{-(\beta-\beta')\skew5\hat{H}}\,{\rm e}^{{\rm i}\skew5\hat{H}(t-t')}\big]\right). \end{eqnarray} \tag{ 18 }$$

Figure 7 shows, with $\hat {B}^{\dagger} =\hat {A}$ (see footnote ⁴), the computation cost per time step and the maximum bond dimensions for the optimized scheme. The maximum reachable times t^opt_max of the optimized scheme are at least twice as large as the reachable times t^B_max(β) for scheme B, specifically t^opt_max(β) ≳ 2t^B_max(β/2). See also figure 8. The factor 1/2 in the temperature argument is inessential as one finds that t_max(β) ≈ t_max(β/2). This is due to the fact that t_max varies slowly as a function of log β for all response functions considered in this paper. Let us discuss the possible scenarios for the optimum values of β' and t' for the case $\hat {B}^{\dagger} =\hat {A}$. The computation cost for scheme B is dominated by the tDMRG computation of the operator $[{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A}\, \mathrm {e}^{-\beta '\skew5\hat{H}/2}\,{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}]$ in equation (12). The corresponding maximum time, reachable for some given computational resources, truncation weight , and inverse temperature β', is denoted by t^B_max(β'). The maximum reachable time for the scheme (18) as a function of β is then given by

$$\begin{equation} t^{{\rm opt}}_{{\rm max}}(\beta) = {{\rm max}}_{0\leqslant \beta'\leqslant\beta}\left( t^B_{{\rm max}}(\beta') + t^B_{{\rm max}}(\beta-\beta') \right). \end{equation}{} \tag{ 19 }$$

There are two possible scenarios, as depicted in figure 9. If t^B_max(β') is concave on the interval [0,β], the optimal scheme corresponds to β' = β/2. Otherwise, there will be an optimum β' ≠ β/2. In the examples studied here, t^B_max(β') was found to be almost concave in the relevant temperature ranges.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One can hence devise a corresponding scheme C that is near-optimal among the schemes (18) (at least for $\hat {B}^{\dagger} =\hat {A}$), does not require any optimization, but outperforms scheme B by a factor of two in the maximum reachable times. It is only outdone by scheme A at very low temperatures. The scheme corresponds to the choice β' = β/2 in equation (18):

$$\begin{eqnarray} \chi^C_{\skew3\hat{A}\skew3\hat{B}}(\beta,t_A+t_B)&=& \frac{1}{Z_\beta}{\rm Tr}\left(\big[\mathrm{e}^{-\beta\skew5\hat{H}/2}\skew3\hat{B}(t_B)\big] \big[\skew3\hat{A}(-t_A)\, \mathrm{e}^{-\beta)\skew5\hat{H}/2}\big]\right)\nonumber\\ &=& \frac{1}{Z_\beta}{\rm Tr}\left(\big[{\rm e}^{{\rm i}\skew5\hat{H} t_B}\,\mathrm{e}^{-\beta\skew5\hat{H}/2}\skew3\hat{B} \,{\rm e}^{-{\rm i}\skew5\hat{H} t_B}\big] \big[{\rm e}^{-{\rm i}\skew5\hat{H} t_A}\skew3\hat{A} \,\mathrm{e}^{-\beta\skew5\hat{H}/2}\,{\rm e}^{{\rm i}\skew5\hat{H} t_A}\big]\right). \end{eqnarray}{} \tag{ 20 }$$

After the imaginary-time evolution that yields $[\mathrm {e}^{-\frac {\beta }{2}\skew5\hat{H}}]$, one runs two real-time tDMRG simulations to obtain MPOs $[{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t_B}\,\mathrm {e}^{-\frac {\beta }{2}\skew5\hat{H}}\hat {B} \,{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t_B}]$ and $[{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t_A}\hat {A} \,\mathrm {e}^{-\frac {\beta }{2}\skew5\hat{H}}\,{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t_A}]$. With equation (20) one then obtains $\chi ^C_{\hat {A}\hat {B}}(\beta ,t_A+t_B)$; see figure 9. The accuracy of the MPOs should be kept under control during the whole simulation, for example, as described in section 2 by bounding the truncation error (see footnote ³). If this is done properly, it is of minor importance which specific t_A and t_B = t − t_A are chosen to evaluate $\chi _{\hat {A}\hat {B}}(\beta ,t)$ for a given time t. The maximum reachable time t is of course given by the sum of the maximum reachable t_A and t_B. For the case $\hat {A}=\hat {B}^{\dagger}$ (or cases where $\hat {B}^{\dagger}$ is for example simply a translate of $\hat {A}$, see footnote ⁴), the maximum reachable times for t_A and t_B are equal, and the total maximum reachable time with this scheme is then twice as large as the maximum time of scheme B. This makes many more physical applications accessible.

As pointed out above, scheme C is optimal among the classes of schemes corresponding to equation (18) if t^B_max(β') is a convex function. For the cases studied here it was indeed found to be convex or almost convex in the considered temperature ranges. One can expect different behavior, for example, for systems with dynamics on different energy scales. In such cases one can achieve considerable performance improvements by optimizing among the class of schemes (17) or (18). In exceptional cases, evolved operators $\hat {A}(t)$ have an MPO representation with a time-independent bond dimension. This implies t^B_max(0) = ∞ and one can compute the response function for arbitrary times by using β' = β and t' = 0 in equation (18). One such case is the operator $\hat {S}_j^z(t)$ for the XY model (J^z = 0), which can for all times be written as an MPO of bond dimension M_j = 4 [60].

The time evolution of matrix product states or density operators, as employed in the presented evaluation schemes for thermal response functions, can be implemented within the tDMRG framework in several different ways. The results for the computation costs of the different evaluation schemes are essentially independent of the chosen time evolution algorithm. The currently most common choice [14–16] employs a Trotter–Suzuki decomposition [61–65]. Alternatively, one can use, for example, the Arnoldi method, the Runge–Kutta method or other Krylov subspace approaches [66–68].

Let us consider in the following a Hamiltonian $\skew5\hat{H}=\sum _{i=1}^L\hat {h}_i$ with nearest-neighbor interactions $\hat {h}_i$, operators $\hat {A}$ and $\hat {B}$ with finite spatial support (for simplicity, single sites), and implementing the time evolution with a Trotter–Suzuki decomposition. This case allows for an alternative evaluation of the thermal response functions on the basis of the transfer matrix renormalization group (TMRG) [36–39]. The Trotter–Suzuki decompositions of real- or imaginary-time propagators read

$$\begin{eqnarray} \mathrm{e}^{-\tau\skew5\hat{H}}&=&\left(\mathrm{e}^{-\Delta\tau(\skew5\hat{H}_{{\rm odd}}+\skew5\hat{H}_{{\rm even}})}\right)^N\nonumber\\ &=&\left(\mathrm{e}^{-a_{m+1}\Delta\tau\skew5\hat{H}_{{\rm odd}}}\prod_{n=1}^m\left(\mathrm{e}^{-b_{n}\Delta\tau\skew5\hat{H}_{{\rm even}}}\,\mathrm{e}^{-a_{n}\Delta\tau\skew5\hat{H}_{{\rm odd}}}\right) + {\mathcal{O}}\left((\Delta\tau)^{p+1}\right) \right)^N, \end{eqnarray} \tag{ 21 }$$

where τ = β or τ = ± it, respectively, Δτ = τ/N is the time step and $\skew5\hat{H}_{{{\mathrm { odd}}},{{\mathrm { even}}}}$ contain the Hamiltonian terms on even and odd bonds, respectively. The coefficients a_n and b_n sum up to one ($\sum _n a_n,\sum _n b_n=1$) and define different decompositions of order p and number of stages m per time step. A common choice is the leapfrog algorithm which has m = 1 stages and order p = 2 with a₁ = a₂ = 1/2 and b₁ = 1. In this work, a decomposition with m = 5 stages and order p = 4 was used. When the operator exponentials occurring in the formula $Z^{-1}_\beta {\mathrm { Tr}}(\mathrm {e}^{-\beta \skew5\hat{H}/2}\, {\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}\hat {B}\, {\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A} \,\mathrm {e}^{-\beta \skew5\hat{H}/2})$ for the response function in scheme A, equation (9), are decomposed by such a Trotter–Suzuki decomposition, one ends up with the task of contracting a two-dimensional (2D) tensor network (quantum cellular automaton) of local gates as displayed in figure 10 with the space direction being horizontal and the time direction being vertical. Similar tensor networks result for the other evaluation schemes.

Zoom In Zoom Out Reset image size

The tDMRG approach scheme A corresponds to starting with a trivial MPO $[\mathbbm {1}]$ that represents the identity with bond dimensions M_i = 1 ∀ i∈[1,L]. Beginning at the bottom, one then contracts one row of local gates after another to that MPO until reaching the row that contains operator $\hat {B}$, yielding $[{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} t}\hat {A}\, \mathrm {e}^{-\beta \skew5\hat{H}/2}]$. Multiplying, in every step of this iteration, one row of gates with the current MPO in an exact manner gives a resulting MPO with increased bond dimensions. A subsequent truncation step reduces the bond dimensions again. The degree to which the bond dimensions are reduced in the truncation steps is controlled by the predefined truncation weight (8), which determines the accuracy of the computation. Similarly, starting from the top, one contracts rows of gates to obtain $[\mathrm {e}^{-\beta \skew5\hat{H}/2}{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} t}]$. Finally, one computes the response function by evaluating the Hilbert–Schmidt scalar product (9) as visualized in figure 1.

In the TMRG approach, the columns of local gates in figure 10 are interpreted as transfer matrices $\hat T=\hat T(\beta ,t)$ (also $\hat T_{\hat {A}}$ and $\hat T_{\hat {B}}$) and operate on a Hilbert space of dimension d^{β/Δβ+2t/Δt}, where d denotes the dimension of the single-site Hilbert space. The last column corresponds to a pure state |ψ_L〉 and the first column to a dual state 〈ψ₀|. Those transfer matrices and states have matrix product representations with bond dimensions M_i ⩽ d² ∀ i. Starting from the left and the right with states ψ₀ or ψ_L, respectively, one transfer matrix after another is applied with intermediate truncation steps, for example, until one reaches the column containing operator $\hat {B}$. The response function is then given by the overlap of the two resulting MPS, $[\langle \psi _0|\hat T\ldots \hat T\hat T_{\hat {A}}\hat T\ldots \hat T][\hat T_{\hat {B}}\hat T\ldots \hat T|\psi _{\mathrm {L}}\rangle ]$.

The advantages of the tDMRG approaches, pursued in this work, are that (i) the number of matrices in the MPOs is fixed by the lattice size L instead of growing linearly with time t and inverse temperature β in the TMRG approach, (ii) the operators $\hat {A}$ and $\hat {B}$ can have non-local support and (iii) one can evaluate the spectral function for all times t, up to the maximum reachable time, by two tDMRG runs, whereas in the TMRG approach one has to carry out a new calculation for every time t of interest.

Concerning (ii), it should be pointed out that non-local operators $\hat {A}$ and $\hat {B}$ would also be accessible in the TMRG approach as long as they are MPOs with sufficiently small bond dimensions. In [41, 42] the focus was on autocorrelation functions of local operators in the thermodynamic limit L → ∞. In this case, one does not need to contract one column after another, but can evaluate the response function directly from the left and right transfer matrix eigenvectors with maximum eigenvalue. When one uses the so-called infinite-system DMRG [1] to obtain those eigenvectors (a single build-up sweep for the finite imaginary-time lattice; not to be mistaken with the algorithm presented in [69] which generates translationally invariant MPS for infinite systems), problem (iii) does not occur and autocorrelation functions can be evaluated for all accessible times t in a single run. One should however be careful in applying infinite-system DMRG only, as this algorithm has several pitfalls [3]. At least for larger times, finite-system DMRG (multiple sweeps) should be necessary.

In this paper, I have studied and explained the computation costs of different tDMRG schemes for the efficient and precise evaluation of finite-temperature response functions of strongly correlated quantum systems. Simplifying the notation to some extent by formulating everything in terms of MPOs elucidated the effects of quasi-locality on the costs. The new class of optimizable evaluation schemes, equations (17) and (18), typically outperforms the earlier schemes from the literature in terms of the maximum reachable times by a factor of ≳2. This gives access to many more physical applications. The novel scheme C, which requires no additional optimization and is near-optimal in many typical examples, can be expected to be the method of choice for future applications. For more complex models, such as systems with dynamics on different energy scales, scheme C still outperforms the older tDMRG methods, but one can achieve further substantial performance gains by first studying the costs for the computation of operators $\big [{\mathrm { e}}^{{\mathrm { i}}\skew5\hat{H} s}\,\mathrm {e}^{-\alpha \skew5\hat{H}}\hat {B} \,{\mathrm { e}}^{-{\mathrm { i}}\skew5\hat{H} s'}\big ]$, and then using this information to determine the optimal parameters for the scheme (17). Finally, it has been argued that the tDMRG schemes are in some respects favorable to corresponding TMRG variants. As a first application, scheme C has been used in [51] to calculate, for 1D bosons in the quantum critical regime with dynamic critical exponent z = 2, the universal scaling function for the thermal spectral function.

Discussions with J Barthel, C Karrasch, A Kolezhuk, I P McCulloch, T Nishino, S Sachdev, J Sirker and U Schollwöck as well as financial support through DFG FOR 801 are gratefully acknowledged.