{'fraction_non_alphanumeric': 0.03557445173068147, 'fraction_numerical': 0.009199875093079676, 'mean_word_length': 4.553154595171402, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Minimally entangled typical thermal states versus matrix product purifications for the simulation of equilibrium states and time evolution

For the simulation of equilibrium states and finite-temperature response functions of strongly-correlated quantum many-body systems, we compare the efficiencies of two different approaches in the framework of the density matrix renormalization group (DMRG). The first is based on matrix product purifications. The second, more recent one, is based on so-called minimally entangled typical thermal states (METTS). For the latter, we highlight the interplay of statistical and DMRG truncation errors, discuss the use of self-averaging effects, and describe schemes for the computation of response functions. For critical as well as gapped phases of the spin-1/2 XXZ chain and the one-dimensional Bose-Hubbard model, we assess the computation costs and accuracies of the two methods at different temperatures. For almost all considered cases, we find that, for the same computation cost, purifications yield more accurate results than METTS -- often by orders of magnitude. The METTS algorithm becomes more efficient only for temperatures well below the system's energy gap. The exponential growth of the computation cost in the evaluation of response functions limits the attainable timescales in both methods and we find that in this regard, METTS do not outperform purifications.

Finite-temperature correlation and response functions of quantum many-particle systems are of great interest. They provide insights into the many-body physics and allow to compare theoretical models to experimental results. However, their precise computation remains challenging. For many relevant models, one has to rely on the development of efficient numerical techniques. The most successful method for the study of strongly correlated one-dimensional (1D) systems is the density-matrix renormalization group (DMRG) which is based on matrix product states (MPS) [1][2][3]. While DMRG was originally designed to study ground states of 1D systems, its extension to the time evolution of quantum states within tDMRG [4][5][6] allows for the simulation of quenches and response functions. Based on this extension, quite different methods for the simulations at finite temperatures have been developed. One of them rests on a purification [44] of the density matrix [7][8][9] which can be encoded in matrix product form [10,11]. First, this was successfully applied to study static finite-temperature properties of quantum spin chains [12,13]. The combination with real-time tDMRG allows for the precise evaluation of finite-temperature response functions and can be applied to compute spectral functions and to study a variety of experimentally relevant systems [14][15][16][17][18][19][20][21][22]. Despite their success, simulations based on such purifications are often limited with respect to the reachable inverse temperatures and/or times, due to a growth of entanglement which is accompanied by a corresponding growth of computation costs. The search for comple-mentary approaches led to an algorithm that avoids the direct encoding of the mixed states: Instead of purifying the density matrix, one can sample cleverly chosen pure states, so-called minimally entangled typical thermal states (METTS) [23,24]. While they can be efficiently encoded in matrix product form as their entanglement is relatively low, they represent well the thermal properties of the system at hand. The METTS algorithm has been successfully applied to study static properties and quantum quenches at finite temperature [25][26][27]. However, a thorough analysis of its accuracy and efficiency compared to computations using matrix product purifications of the density matrices was lacking.
In this work, we discuss how self-averaging can be used to moderately reduce statistical errors in the METTS algorithm and introduce schemes for the evaluation of response functions using METTS. We compare the accuracies and computation costs of the METTS and purification approaches for the evaluation of finite-temperature correlation and response functions. We focus on two paradigmatic models of interacting quantum systems, namely the spin-1/2 XXZ chain [28][29][30] and the 1D Bose-Hubbard model [31,32] at critical as well as non-critical points of their phase diagrams. In contrast to indications and expectations expressed in the earlier literature [23,27,33], for almost all cases considered here, we find that, for the same computation cost, the purification approach yields more accurate results than METTS -often by orders of magnitude. METTS become more efficient only for temperatures well below the energy gap of the system. It would be interesting to investigate further whether other approaches as, for example, computations based on the ground state and a few excited states, which can be determined variationally, could outperform the METTS approach for such very low temperatures. For the comparisons, we always use equal total computation costs for both methods, ignoring that METTS simulations can be parallelized more easily than purification simulations by generating independent Markov chains on different computing nodes. This can be taken into account by keeping in mind that, to reduce the presented METTS errors by an order of magnitude, one needs to increase the number of employed computing nodes by at least a factor of 100.
The article is structured as follows. Sections II and III shortly review the algorithms for computing static observables with purifications and METTS, respectively, and discuss the interplay of statistical and truncation errors for METTS (Sec. III B) as well as self-averaging (Sec. III C). In section V, we introduce for METTS a simple scheme and two more elaborate schemes for the computation of response functions which are analogous to corresponding schemes based on purifications [14,16,17]. The main objective of the paper, the efficiency comparison of METTS and purifications, is presented in section IV for static correlations functions in the spin-1/2 XXZ chain and the Bose-Hubbard model, and in section VI for response functions in the XXZ model. Some technical issues are described in appendices. We summarize and conclude in section VII.
A state |Pρ ∈ H ⊗ H aux is called a purification of the density matrixρ on H if Tr aux |Pρ Pρ| =ρ. (1) Choosing the auxiliary Hilbert space H aux isomorphic to the physical Hilbert space H, i.e., H H aux , it is simple to give a purification of the infinite-temperature statê where |σ i are orthonormal basis states for lattice site i, and |σ i aux for the corresponding lattice site of the auxiliary system. For the orthonormal basis {|σ = i |σ i } of H, let |X ∈ H ⊗ H aux denote the vectorization of an operatorX on H such that In this notation, we have that |ρ β/2 ∈ H ⊗ H aux is according to equations (1) and (3) a purification of the density matrixρ β . Because |ρ 0 ≡ |Pρ 0 , as given in Eq. (2), is a product state, it can be encoded as an MPS with matrices of bond dimension one (cf. appendix A). With |ρ 0 as the initial state, one can employ imaginary-time evolution, to obtain purifications |ρ β/2 for finite-temperature stateŝ ρ β , To this purpose one can employ the time-dependent DMRG algorithm (tDMRG) [5,6] or the almost identical time-evolved block decimation (TEBD) [4]. Specifics of our simulations are summarized in appendix A. Exploiting that Z β = ρ β/2 |ρ β/2 andρ β = Tr aux |ρ β/2 ρ β/2 |, thermal expectation values can be computed by where both physical and auxiliary degrees of freedom are summed over.

A. Algorithm for static observables
The strategy employed in the minimally entangled typical thermal states (METTS) algorithm is to approximate thermal expectation values Ô β by sampling pure quantum states that have two favorable properties. They represent well the physical properties of the system for the given temperature, and the entanglement of the states is relatively low which makes DMRG calculations efficient. Expressing the trace for the thermal expectation value using some orthonormal basis {|n } of product states where |n i are (arbitrary) orthonormal basis states for lattice site i, we have Defining the METTS |φ n and their probabilities P n , the thermal average reads Thus, by sampling the states |φ n according to the probability distribution P n /Z β , we can approximate Ô β by averaging over φ n |Ô|φ n . The computation cost of DMRG is directly related to the entanglement of the quantum state (see appendix C). Hence, product states (6) are a natural choice because their entanglement entropy is zero and it remains reasonably low during imaginary-time evolution. The sampling is accomplished efficiently by generating a Markov chain of METTS as illustrated in Fig. 1. An arbitrary initial product state |n is evolved in imaginary time to obtain the METTS |φ n . (See appendix A for details on the tDMRG evolution.) Then, the METTS is collapsed through a projective measurement of the operator n |n n |, yielding a new product state |n with probability p n n := | n |φ n | 2 from which one subsequently computes |φ n and so on. The transition probabilities obey detailed balance p n n P n = | n |e −βĤ/2 |n | 2 = p nn P n (10) such that the desired distribution P n is indeed the fixed point of this Markov process. Note that the projective measurement, |φ n → |n , can be carried out sequentially, site by site. Starting at some site i, we go from |φ n to |n i n i | · tDMRG tDMRG |φ n / π(n i ) with probability π(n i ) := n i |φ n 2 . Measuring subsequently on site j, we go to |n i n j n i n j | · |φ n / π(n i )π(n j |n i ) with probability π(n j |n i ) := n i n j |φ n 2 /π(n i ) such that, in the end, we arrive at state |n indeed with probability p n n = π(n i )π(n j |n i )π(n k |n i n j ) . . . = | n |φ n | 2 .
Due to this, the projective measurement of MPS |φ n can be done efficiently in a single sweep through the lattice [23,24].
In order to ensure ergodicity and reduce autocorrelation times, it is useful to switch between different measurement bases {|n (k) } during the sampling. Details on this and our corresponding choice are described in appendix B.

B. Errors: Statistics, truncations, and Trotter
There are two error sources for the evolution of MPS in the framework of tDMRG as described in appendix A. The first is due to truncations of low-weight terms in the Schmidt decomposition of the wavefunction. This error is well-controlled by the truncation threshold (we call it P for purifications and M for METTS) which bounds, in every time step, the two-norm deviation ψ trunc − ψ of the truncated MPS from the exactly evolved state. We implement tDMRG using fourth-order Trotter-Suzuki decompositions of the evolution operators with time steps ∆t of size 0.125 for purifications and 0.05 for METTS ( = 1). These are the second error source. The resulting errors of order ∆t 5 can only become relevant for very large times and we made sure that they are never dominant for the presented data. Additionally, the accuracy of the METTS sampling algorithm is influenced by a third error source -the statistical error that depends on the number of samples N used for averaging. Fig. 2 illustrates the interplay of statistical errors and truncation errors in a METTS computation of the correlator S − 3 S + 0 β for the exactly solvable 1D XX model [∆ = 0 in Eq. (11)] at inverse temperature β = 4. The exact solution is shortly described in appendix D. The top panel shows the convergence of the METTS algorithm for different truncation thresholds M as a function of the number of samples. For low sample numbers, the statistical error dominates and the truncation error is negligible. Autocorrelation times between subsequent samples are short and the statistical error is to a good approximation proportional to 1/ √ N and independent of the truncation threshold. Once the sample number reaches a certain M -dependent threshold, the statistical error has reduced to a magnitude that is comparable to the error induced by the truncations. The curves begin to level out as the relative contribution of the truncation error grows. At a certain point, further samples will not enhance the accuracy of the simulation as the truncations prevent further convergence to the correct value. Of course, the less we truncate (the lower M is), the longer the 1/ √ N -convergence persists. The truncation affects the accuracy in two ways. As every METTS is approximated by a truncated MPS, the values one obtains for each sample are not exact. Additionally, the produced samples will not correspond exactly to the correct probability distribution P n of Eq. (9). This is because the transition probabilities depend on the samples and are thus also affected by the truncations.
While lowering the truncation threshold M yields more accurate results, it also increases the computation cost per sample. In the lower panel of Fig. 2, we present the same errors of the METTS algorithm as in the top panel, but here as a function of the total computation cost which we quantify in an implementation-independent way as described in appendix C. This shifts the simulations with lower truncation thresholds and thus more costly samples to the right. Whereas, for a fixed number of samples (top panel), the accuracy is a monotonic function of M , this is not necessarily so for fixed total computation cost (lower panel). When plotted against the total computation costs, accuracy curves of METTS simulations with different M have crossings. For practical simulations, it is therefore important to choose the truncation threshold such that the two error sources are balanced. In the lower panel of Fig. 2, one can easily read of the truncation threshold that is optimal for a given computation cost. While the optimal truncation threshold depends on the specific system studied and the observable that is evaluated, it generally shifts towards lower values of M with increasing total computation cost.

C. Exploiting self-averaging
If the considered model is translation invariant, we are free to choose an arbitrary position in the lattice for the evaluation of an observable. While the average of these expectation values will converge to the correct result independent of the position, a single METTS sample itself is not translation invariant. We can thus exploit selfaveraging to reduce statistical errors in the METTS algorithm. For a correlation length ξ, averaging a local observable Ô x0 β over a block of L sites x 0 corresponds for high and intermediate temperatures to O(L /ξ) statistically independent samples. Therefore, we can expect that the statistical METTS errors reduce by a factor of order ξ/L . For systems with open boundary conditions, one has to restrict the averaging to sites x 0 with a sufficient distance from the boundaries.
We illustrate the effect for the spin-1/2 XXZ chain in  The observed error reduction is of the expected order of magnitude. For ∆ = 3 and β = 4, the impact of selfaveraging is rather small for short distances i. This is due to the fact that, in this case, the temperature is already well below the gap (cf. Table I) and all excitations occurring in the METTS are of long wavelength. Hence, short-range correlations in the METTS are almost translation invariant.
As the additional computation cost for the spatial averaging of time-local observables (in equilibrium or quench dynamics) is negligible, it is advisable to enhance the METTS accuracy through the self-averaging whenever finite-size effects are well-controlled. Spatial averaging in the evaluation of response function would however require additional real-time evolutions and seems hence not useful.

IV. COMPARISON FOR STATIC CORRELATORS
A. Procedures to compare efficiencies In this section we compare the efficiencies of METTS and matrix product purifications, by studying accuracies of static thermal correlation functions for fixed computation costs. To this purpose, accuracies are quantified by the deviations of the obtained expectation values from exact or quasi-exact data as described in appendix D. For the error of N METTS, we generate several sets of N subsequent samples, and take the root mean square of the average deviations from the reference data in each set. The computation cost is quantified in a largely implementation-and platform-independent way as a function of the MPS bond dimensions D i = D i (β, t) as detailed in appendix C. For METTS, we average the computation costs of the sample sets.
To take the interplay between statistical and truncation errors (section III B) into account in the assessment of the performances of METTS and purifications, we proceed as follows. We choose a truncation threshold P for the purification and determine the theoretical computation cost of the simulation. Then we produce METTS samples, using different truncation thresholds M , and fix the sample set size N for each M by the quotient of the purification cost and the average cost per METTS sample such that, for each truncation threshold, the total computation costs of both methods are equal. In the figures, we only present the results for the truncation thresholds M that yield the best results, i.e., those M that approximately minimize the error for fixed computation cost, as well as the results for two nearby values of M .
Setting the total computation costs of both methods equal, ignores that METTS simulations can be parallelized more easily than purification simulations, by generating independent Markov chains on different computing nodes. However, this can be easily taken into account. To reduce the METTS errors, as presented in the following, by an order of magnitude, one has to increase the number of employed computing nodes by at least a factor of 100. This assumes that the extent of thermalization phases at the beginnings of the Markov chains is negligible, and the factor 100 is a lower bound because one also needs to decrease M (hence, increasing the computation cost per METTS) when the statistical error is being reduced.

B. Spin-1/2 XXZ chain
For the spin-1/2 XXZ chain [28][29][30] with Hamiltonian let us consider three values of the anisotropy parameter: the exactly solvable non-interacting case ∆ = 0 and the isotropic Heisenberg antiferromagnet at ∆ = 1, which are both critical (gapless), as well as the point ∆ = 3 in the gapped Néel phase. The lattice has size L and site i = 0 is at the center of the chain. We apply the METTS and purification algorithms, as described in sections II and III A, to compute the static correlation function Ŝ −   Table II. the errors of the purification simulations with P = 10 −10 and 10 −12 and the errors of several METTS simulations with different truncation thresholds M . All curves that appear within a panel and refer to the same temperature are based on simulations of equal total computation costs. For each panel, the METTS truncation thresholds M are chosen such that, for the largest M , the truncation error dominates; for the lowest, the statistical error dominates; and for the intermediate M , statistical and truncation errors are balanced in an optimal way such that the error is (approximately) minimized for the given computation cost.
For almost all parameters, the matrix product purification simulations yield more accurate results than the best METTS computations. For β = 4, the errors of the methods differ by up to a few orders of magnitude. When lowering the temperature, entanglement and correlation lengths increase. The increased absolute value of the correlation function is reflected in a correspondingly larger absolute error for the simulations based on purifications. The efficiency of the METTS sampling can increase when lowering the temperatures, especially when it gets suffi-ciently below the energy gap ∆E. In this case, the dimension of the relevant state space to be sampled by METTS is strongly reduced and hence is the statistical error. This is confirmed for the lower temperature β = 16 in Fig. 4. For the critical (gapless) systems, the purification approach is still more precise. For the gapped system (∆ = 3), the temperature is with β = 16 already substantially below the energy gap, β∆E ∆=3 ∼ 10 (cf .  Table I), and METTS can in this case indeed outperform the matrix product purification.
That METTS can become more efficient than the purification is most obvious for the limit β → ∞. In this case, every METTS |φ n for which the initial state |n has nonzero overlap with a ground state will simply be this ground state, |φ n = |gs (up to truncation errors). The purification, on the other hand, evolves to the purification |gs ⊗ |gs aux of the ground state density matrix |gs gs|. As the tensor product of two MPS of bond dimension D is an MPS with bond dimension D 2 , for the same accuracy, β → ∞ computation costs for the METTS are reduced by roughly a factor of D 3 [45].
The error of the purification changes significantly as a  Table III. function of the distance i. Here, the deviation from the (quasi)-exact reference data is mainly due to the truncations. Generally, when we evaluate a correlator based on a truncated matrix product state, the error can grow at short distances, reaching an P -dependent maximum, before starting to decay as the absolute value of the correlator itself becomes very small. On the other hand, Spin-1/2 XXZ chain  METTS simulation errors often remain constant at large distances. Even for the exponentially decaying correlation functions, the METTS error stays well above zero. This is clearly a signature of the statistical error. Cases where the METTS error decays with distance, are usually situations where the truncation error dominates over the statistical error. For the computation costs chosen in our study, this is seen for temperatures well below the energy gap, but one can also observe this behavior when choosing very large truncation thresholds.
Finally, let us shortly discuss an ergodicity issue in the gapped Néel phase. For β = 16, the system is essentially in its ground state. Hence, the weight of the two degenerate Néel states in the thermal state becomes significant. The purified state still obeys the Z 2 symmetry. Every METTS sample is some linear combination of the ground states. Depending on the choice of the collapse basis {|n }, severe ergodicity problems can occur in the METTS algorithm. As can be shown, even for the rotation-symmetric random bases that we use (described in appendix B) and which seem to have very good ergodicity properties at first sight, the transition probabilities from one of the degenerate ground states to the other decay exponentially in the system size L. Several remedies are available. One can for example avoid this issue by symmetrizing every METTS before measurements. However, in the case of response functions (see section V), this would require an additional real-time evolution and, hence, roughly a doubling of the costs. Alternatively, one could, e.g., use theŜ x and theŜ y eigenbases for the METTS collapse. Here, we decided to keep the random collapse bases and implicitly average over the two Néel states by using symmetrized observables; for Ŝ − i (t)Ŝ + 0 (0) β , according to such that the measurement is not sensitive to the broken symmetry.

C. 1D Bose-Hubbard model
The Hamiltonian of the one-dimensional (1D) Bose-Hubbard model [31,32] is given bŷ = δ ij and the number operatorsn i =b † ib i , hopping J, onsite interaction U , and chemical potential µ.
We set µ = 0.5 and U = 1, and consider three values of the hopping parameter. J = 0.09375 corresponds to the Mott-insulating phase with one particle per site, J = 0.125 places the system close to the phase boundary between the Mott-insulating and the superfluid phase, and J = 0.25 lies in the superfluid phase. In Fig. 5, we show the comparison between METTS and purifications for the Bose-Hubbard model at β = 4 and β = 16. Again, all curves in a given panel that refer to the same temperature share the same computation cost. In general, the results are similar to those for the XXZ model. However, the differences between the METTS and purifications errors are larger in the comparison for the Bose-Hubbard model. This can be explained with the smaller energy gaps, compared to those in the spin-1/2 XXZ model. See Table I. For the Bose-Hubbard model at the specified points in the phase diagram, errors of the METTS simulations exceed those of purification simulations by two to three orders of magnitude.

V. ALGORITHMS FOR RESPONSE FUNCTIONS
The matrix product purification and METTS methods, as described in sections II and III A, can be extended to withX(t) ≡ e iĤtX e −iĤt . For purifications, several schemes have been suggested and analyzed [14,[16][17][18]. They have different properties concerning the dependence of computation costs on β and t. According to the naming introduced in Refs. [17,18], we address below computation schemes A, B, and C, introduce corresponding schemes for the METTS framework, and discuss their properties. In Refs. [17,18], further optimizable classes of schemes have been studied. In comparison to the near-optimal scheme C, they allow to substantially reduce computation costs, for example, for systems with separated energy scales. As there is probably no useful adaption of them for METTS, they will not be addressed here. In related work, Pižorn et al. [34] recently discussed the option of working, for matrix product purifications (equivalently, one can think in terms of matrix product operators [18]), in the Heisenberg picture, i.e., to compute the vectorization |X † (t) of the evolved operatorX † (t) and the vectorization |ρ β of the thermal stateρ β to then obtain X (t)Ŷ β = X † (t)|Ŷ |ρ β /Z β . One may notice that this is just a special case of the class of optimizable schemes suggested in Eq. (9) of Ref. [17] or Eq. (18) of Ref. [18]. We will not discuss it further because this  Heisenberg picture scheme is in general computationally suboptimal and has no direct METTS equivalent. Specifically, the cost analysis in Refs. [17,18] suggests that the Heisenberg picture scheme will typically reach half the maximum times that can be reached by the optimized schemes or the near-optimal scheme C to be described below.
A. The simple scheme A Starting from the matrix product purification |ρ β/2 [Eq. (4)] of the density matrixρ β , according to evaluation scheme A [14] for the response function (13), we first compute matrix product representations of e −iĤt |ρ β and e −iĤtŶ |ρ β using tDMRG. The response function is then given by the matrix element The vectorization (3) As done in Eqs. (14) and (15), the evolved MPS or MPOs that are used for the evaluation of the response function are always indicated by square brackets in the following.
The METTS equivalent of scheme A is illustrated in Fig. 6. For every sample |φ n , we compute φ n |e iĤt X e −iĤtŶ |φ n (16) using real-time evolution and then average over the values obtained for each sample. In practice this means that we have to carry out two independent tDMRG simulations, using |φ n andŶ |φ n as initial states, up to some maximum time. At any intermediate time-point that we are interested in, we can evaluate the response function X (t)Ŷ β -even for a set of operators {X i } if we wish.

B. Scheme B is not useful for METTS
In the context of matrix product purifications, we have the alternative scheme B [16] which reads in the MPO representation . (17) In comparison to scheme A (15), it corresponds to shifting exp(iĤt) from the first to the second MPO. In Ref. [18], it was explained why this scheme has some advantages at higher temperatures and disadvantages at lower temperatures. The METTS equivalent would be to compute φ n | X e −iĤtŶ e iĤt |φ n .
As exemplified in Fig. 7 and explained in the following, its computation costs are unfortunately strictly higher than those of scheme A (16). Both schemes share the cost for computing [e iĤt |φ n ]. In scheme A, the evolution of [e −iĤt Ŷ |φ n ] for t = 0 → t has usually about the same cost, and the entanglement in both states increases typically linearly with t everywhere in the system, reaching some value E t . In scheme B (18), the required evolution of [e −iĤt Ŷ e iĤt |φ n ] for t = 0 → t starts with the high entanglement value E t . Due to quasi-locality [35,36], the entanglement will then reduce in regions of the lattice that are at sufficient distance from the spatial support of operatorŶ . It will however remain high or even increase in the vicinity ofŶ .

C. Reaching longer times with scheme C
Among the classes of optimizable evaluation schemes studied in Refs. [17,18], the near-optimal scheme C was found to be very useful and reaches about twice the maximum times reachable with schemes A and B. In the MPO representation, it reads In a similar way, we can compute the response function with METTS, by averaging over which makes it possible to reach, for the same computational resources, times approximately twice as big as in scheme A.

D. Comparison
Due to the cyclic property of the trace and the (trivial) fact that the evolution operators commute with density matricesρ β , the results of all three schemes converge up to differences in truncation and Trotter errors (see appendix A) to the same thermal average (13). As the truncations in the real-time evolutions are well controlled, the errors of the simulations based on different schemes do not differ significantly if all other parameters are left unchanged. Statistical errors can however differ somewhat [46]. Fig. 7 compares the computation costs of the different evaluation schemes in the METTS algorithm for an autocorrelation function in the isotropic Heisenberg antiferromagnet [∆ = 1 in Eq. (11)]. While scheme B proves to be the least efficient for the parameters considered here, scheme C indeed reduces the computation costs significantly and thus allows for the evaluation of longer maximum times with METTS.
Nevertheless, we choose to use scheme A in this work because it allows to evaluate response functions X i (t)Ŷ β for a whole set of operators {X i } and all times t (up to the maximum reachable time) with only two tDMRG runs per METTS sample. This is very useful for the analysis of response functions like Ŝ − i (t)Ŝ + 0 (0) β , studied in section VI. In contrast, scheme C requires, for every METTS sample, separate tDMRG runs for every required time t and operatorX i . This is also a drawback in comparison to scheme C for purifications [Eq. (19)], where one obtains results for all times t with only two tDMRG runs. Fig. 8 compares the accuracies of the response function Ŝ − i (t)Ŝ + 0 (0) β in the XXZ chain as computed via METTS and purification simulations for times t = 2.5 and 10. The structure of the plots is as in Fig. 4, i.e. the columns again refer to the three values of the anisotropy parameter ∆ = 0, 1, 3 already considered in section IV B and all curves that appear within a panel and refer to the same temperature are based on simulations with the same total computation costs. As the thermal response functions are complex-valued, we show both the real and the imaginary part of the (quasi-)exact results.

VI. RESPONSE IN THE SPIN-1/2 XXZ CHAIN
Generally, the error curves closely resemble the results we obtained for the computation of static correlators in Fig. 4. The purifications are more efficient than METTS except for the gapped system (∆ = 3) at the low temperature β = 16. Fig. 9 displays the errors and computation costs of both methods as a function of time t in the evaluation of response functions for ∆ = 0 at inverse temperature β = 16. Each column refers to a truncation threshold M for the METTS algorithm and shows the errors (top) and computation costs (bottom) for different numbers of samples. The purification errors and costs of three different truncation thresholds P are shown in each column.
The convergence of the METTS errors is clearly reminding of the behavior shown in Fig. 2. As a function of the number of samples, the error is approximately proportional to 1/ √ N , until it saturates to an M -dependent threshold that can be lowered by decreasing the truncation threshold.
The exponentially growing computation costs in both the METTS approach as well as the matrix product purification approach limit the reachable maximum times, i.e., at least qualitatively, METTS seem to have no favorable properties in this respect.

VII. CONCLUSIONS AND DISCUSSION
We have studied properties of the METTS algorithm. There is an interplay of statistical and DMRG truncation errors that one should take into account for efficient simulations. While the optimal truncation threshold M for METTS depends on the specific system studied and the observable that is evaluated, it generally shifts towards lower values of M with increasing total computation cost. As demonstrated, one can also exploit self-averaging of (approximately) translation invariant systems to reduce statistical errors in METTS simulations for static observables.
We have presented a simple scheme for the evaluation of response functions using METTS and two more elaborate schemes, one of which gives access to longer maximum times but needs a separate simulation for every required point in time.
For spin-1/2 XXZ chains and the 1D Bose-Hubbard model, we have compared the accuracies and computation costs of the METTS and purification approaches for the evaluation of finite-temperature correlation and response functions. For almost all cases considered here, we found in contrast to indications and expectations expressed in the earlier literature that, for the same total computation cost, the purification approach yields more accurate results than METTS -often by orders of mag- nitude. METTS become more efficient only for temperatures well below the energy gap of the system (β∆E 10 in our case). For both methods, we have discussed the temperature dependence of the accuracies, and for correlators also their distance dependence.
It would be interesting to investigate further whether other DMRG approaches for the low-temperature regime, such as computations based on the ground state and a few excited states or a sampling that is restricted to the complement of the ground state space, could outperform the METTS in these cases. ends of the chain are one such that the matrix product ( i A σi i ) evaluates to a scalar. When a time-step evolution operator such as exp(−a k ∆τĤ even ) is applied to the MPS, the bond dimensions D i increase toD i ≥ D i -partly for technical reasons, because the resulting state needs to be brought to MPS form, and partly because the entanglement in the state can grow. States with higher entanglement generally require larger bond dimensions. Hence, it is necessary to compress the matrix product representation, i.e., to make a controlled approximation such that the bond dimensions are again reduced to some extent. This can be achieved through a Schmidt decomposition |ψ = D j=1 λ j |j L ⊗ |j R of the state [9] which boils down to doing singular value decompositions of the tensors A i . The corresponding reduced density matrices for the left and right parts of the system are j λ 2 j |j L j| L and j λ 2 j |j R j| R , respectively. The bond dimension for the given bipartition of the lattice is then reduced from D to some value D <D by retaining only the D largest Schmidt coefficients λ j and truncating all smaller ones.
The resulting two-norm error is In our simulations, we control errors due to truncations by discarding in every step of the algorithm only components with density matrix eigenvalues λ 2 j below a predefined truncation threshold. In this way, the bond dimensions become functions of (inverse) temperature and time, D i = D i (β, t). For purifications, the truncation threshold is denoted by P and by M for the METTS. In section III A, we reviewed the METTS algorithm for a fixed orthonormal basis {|n }, used in the projective measurements to generate a new product state |n from the current METTS state |φ n ≡ e −βĤ/2 |n / √ P n . In order to ensure ergodicity and reduce autocorrelation times, it is useful to switch between different measurement bases {|n (k) } during the sampling [24]. This is possible, because it simply corresponds to using an overcomplete basis in the representation of the trace in Eq. (7), i.e., to starting instead from the expression where π k ≥ 0 are probabilities ( k π k = 1) according to which we choose basis {|n (k) } to collapse a METTS. The rest of the derivation remains unchanged except for extending the detailed balance to the choice of the measurement bases.
In our simulations, we found the following procedure to be robust and efficient. Before every collapse of the wavefunction, we randomly determine measurement bases {|n (a n,σ + ib n,σ )|σ with a n,σ and b n,σ drawn from the standard normal distribution. Using the Gram-Schmidt process, these d states are subsequently orthonormalized to obtain {|n (k) i }.

Appendix C: Quantification of computation costs
The main contribution of this work is to compare the efficiency of the METTS and matrix product purification algorithms. In order to make this comparison as independent as possible from specific details of the implementation and computer architecture, we approximate the computation costs associated with the imaginary-and realtime tDMRG evolutions by integrals over a function of the bond dimensions D i = D i (β, t) which were described in detail in appendix A. As the elementary steps in the employed DMRG algorithms are matrix multiplications and singular value decompositions for the MPS matrices A σi i [Eq. (A2)], the computation cost scales in leading order with the third power of the bond dimensions. How- ever, small bond dimensions can lead to substantial subleading contributions corresponding to copying matrices etc. In a realistic assessment of the actual computation cost, these need to be taken into account, especially as optimal METTS bond dimensions can be quite small in some cases. Therefore, we model the computation cost c(D), associated with a single local tDMRG step as a function of the bond dimension D, by c(D) = D 3 + aD 2 + bD. (C1) As shown in Fig. 10, a fit to actual measured computation times yields coefficients a = 7.6 and b = 17.6. While the exact values of a and b depend to some extent on the particular implementation at hand and the employed computer architecture, deviations of the coefficients do not decisively change the results of our study. In essence, the coefficients only depend on the cost scaling of the employed linear algebra routines that are implemented in the highly optimized LAPACK library [39] for the given computer architecture.
Based on this, we can model the computation costs for an imaginary-time evolution from τ = 0 to β/2 and a subsequent real-time evolution to time t, as they occur in the different algorithms (see sections II, III A, and V), by where i runs over all bonds of the lattice. For schemes that require a certain sequence of imaginary-and realtime evolutions, we sum the corresponding contributions to quantify the total cost. For METTS, we use the aver-age computation cost per METTS and multiply for comparisons by the considered number of METTS. Fig. 11 shows that the theoretical computation cost (C2) provides a realistic description of the DMRG computation cost for all simulations, independent of the specific system parameters. Over several orders of magnitude, the agreement is very good; only for extremely short computation times, one obtains slight deviations. This additional overhead is mainly due to the computation time needed for the collapse procedure and the evaluation of observables.
Appendix D: Exact versus quasi-exact reference data In our assessments of efficiencies, we need to quantify the errors of observables and for this purpose, reliable exact or quasi-exact data is indispensable. For the noninteracting point ∆ = 0 of the XXZ model (11), the XX model, we can evaluate all observables exactly. As the system is for this case in a Gaussian state, all correlation functions are determined by the single-particle Green's function due to Wick's theorem [40]. Exploiting this, we can obtain static correlations as well as response functions by evaluating Pfaffian determinants of matrices that contain elements of the single-particle Green's function [41][42][43]. The colored curves in Fig. 12, marked with triangles and circles, show the errors of matrix-product purification simulations for different truncation thresholds P as compared to such exact results.        Although the XXZ model (11) is Bethe-Ansatz integrable [28] for all values of ∆, no exact expressions for longer-ranged static correlators or response functions are available. The same is true for the Bose-Hubbard model (12) which is not (known to be) integrable. For these cases, we use very accurate purification simulations to obtain quasi-exact reference data. In particular, we used the truncation threshold P = 10 −14 during imaginary-time evolutions and P = 10 −12 during (subsequent) realtime evolutions. ( P has been defined in appendix A.) To justify this procedure, we show in Fig. 12 with black solid lines (no symbols) also the errors of the other simulations compared to this quasi-exact data. These curves deviate indeed only very slightly from the colored curves that correspond to the comparison to the exact results.