Exponential Thermal Tensor Network Approach for Quantum Lattice Models

We speed up thermal simulations of quantum many-body systems in both one- (1D) and two-dimensional (2D) models in an exponential way by iteratively projecting the thermal density matrix $\hat\rho=e^{-\beta \hat{H}}$ onto itself. We refer to this scheme of doubling $\beta$ in each step of the imaginary time evolution as the exponential tensor renormalization group (XTRG). This approach is in stark contrast to conventional Trotter-Suzuki-type methods which evolve $\hat\rho$ on a linear quasi-continuous grid in inverse temperature $\beta \equiv 1/T$. In general, XTRG can reach low temperatures exponentially fast, and thus not only saves computational time but also merits better accuracy due to significantly fewer truncation steps. We work in an (effective) 1D setting exploiting matrix product operators (MPOs) which allows us to fully and uniquely implement non-Abelian and Abelian symmetries to greatly enhance numerical performance. We use our XTRG machinery to explore the thermal properties of Heisenberg models on 1D chains and 2D square and triangular lattices down to low temperatures approaching ground state properties. The entanglement properties, as well as the renormalization group flow of entanglement spectra in MPOs, are discussed, where logarithmic entropies (approximately $\ln\beta$) are shown in both spin chains and square lattice models with gapless towers of states. We also reveal that XTRG can be employed to accurately simulate the Heisenberg XXZ model on the square lattice which undergoes a thermal phase transition. We determine its critical temperature based on thermal physical observables, as well as entanglement measures. Overall, we demonstrate that XTRG provides an elegant, versatile, and highly competitive approach to explore thermal properties in both 1D and 2D quantum lattice models.


I. INTRODUCTION
Efficient simulations of interacting quantum many-body systems are crucial for a better understanding of correlated materials. In particular, accurate computation of thermodynamic quantities including magnetization, heat capacity, magnetic susceptibility, etc., enables a direct comparison to experiments and helps to identify relevant microscopic models. The exotic quantum matter includes Luttinger liquids in one (1D) [1,2] and spin liquid materials in two dimensions (2D) [3][4][5][6][7][8]. Besides, the exploration and understanding of the rich and diverse behavior of quantum many-body physics at different energy or, equivalently, temperature scales are interesting from a theoretical perspective. One example are thermal states near 1D quantum critical point which show universal entropy in the partition function due to emergent conformal symmetry [9][10][11] in the low-energy regime. Another prominent example is the thermal fractionalization in the honeycomb Kitaev spin liquid (KSL) at finite temperature [4], which has been experimentally observed [12] in proximate KSL material α-RuCl 3 [13,14].
At a first glance, the simulation of thermal many-body states seems a task more than challenging. There exist exponentially many excited states in the energy spectrum, many of which possess volume-law entanglement and deny any efficient representation in classical computers. However, it turns out that the ensemble density operator, say e −βH with β ≡ 1/T being the inverse temperature, can be efficiently expressed and manipulated in terms of thermal tensor network (TTN) states. The matrix product operator (MPO) is a very natural TTN for describing 1D quantum systems at finite temperature [i.e., (1+1)D], due to the "entanglement" area law in thermal states of both gapped and gapless systems with local interactions. Intuitively, thermal fluctuation effectively "opens" an excitation gap and introduces a finite correlation length in mixed states, rendering an area law in terms of total correlation [15] (as well as operator space entanglement [16]). However, it was estimated that the required MPO bond dimension has an upper bound scaling as D ∼ e β [17] which Finite-size spectra Es with low-energy level-spacing δE ∼ 1/L. By requiring T > δE for thermal averaging, this suggests L β, and therefore a thermal correlation length ξ ∼ β. (b) A large or infinite system has an effective thermal cutoff ξ ∼ β in system length when measuring local observables. Therefore, provided ξ L finite systems can be used, as a very good approximation, to simulate thermal properties in the thermodynamic limit.
still seems to pose a severe barrier towards obtaining low-T properties.
To estimate the computational cost in thermal simulations, one can introduce a formal entanglement entropy in the TTN, e.g., in the MPO representation of the mixed state density matrix, as introduced in Refs. [16,31]. It has been revealed that this MPO entanglement saturates for gapped systems and scales logarithmically (as c 3 ln β) for quantum critical spin chains [16,31]. Very recently, two independent works [32,33], deployed conformal field theory (CFT) arguments to show on general grounds that the Renyi entropy of thermal states of effective 1D systems scales as S (n) E ∼ c 6 (1 + 1 n ) ln β. In the limit n → 1 this implies that also the von Neumann entropy scales like S E ∼ c 3 ln β for thermal states in 1D on general grounds.
Intuitively, this scaling can be understood simply by considering finite-size spectra with many-body low-energy level spacing δE ∼ 1/L, as schematically depicted in Fig. 1. In order to sample a thermal average, it must hold that at the very least T δE, or equivalently L β [in other words, finite β introduces an effective cutoff ξ ∼ β of system size; see Fig. 1(b)]. Now given that low-energy states violate a strict area law via logarithmic corrections, one has the block entropy S E c 6 ln L + const [34] for individual low-energy states using open boundary conditions (OBC). By choosing L = aβ with fixed constant a 1, and by going from individual low-energy pure states |s to a thermal state with weights ρ s , i.e., |s → s ρ s |s s|, the block entropy for the outer product |s s| acquires another factor of 2. Thermal averaging does not change this scaling due to the subadditivity of the von Neumann entropy [see Appendix A], given the constraint L = aβ with fixed a. By thermal averaging over a similar set of low-energy states, the MPO block entropy at the center of the system saturates by further increasing L aβ at a finite value, i.e., is cut off by and, importantly, becomes independent of L. This block entanglement entropy of the thermal state scales similar versus β to the block entropy of a ground state calculation versus L using periodic boundary conditions (PBC). The above intuitive argument fits the holographic picture in terms of thermal multiscale entanglement renormalization ansatz (MERA) [35], where the minimal surface (of half the system) in thermal MERA, as well as the corresponding (bipartite) entanglement entropy, is argued to be proportional to ln β [36].
Furthermore, the argument of translating the scaling of the entropy in ln L to that of ln β is completely consistent with the notion within CFT that β and L are equivalent directions connected via a modular S transformation. This has direct consequences for conformal TTN framework in (1+1)D, i.e., with one spatial (horizontal) and one imaginary time (vertical) axis. The horizontal transfer matrix e −τ H across different temperatures has the ground state of H as its dominant eigenvector which thus contains logarithmic entanglement. For thermodynamics of an infinite-size quantum chain (L β 1), we therefore expect that the vertical transfer matrix (across different real-space sites) also has a dominant eigenvector with entanglement entropy S E ∼ ln(β). In addition, however, by definition of the partition function, it acquires intrinsic PBC in the direction of temperature, which therefore doubles the prefactor in entanglement entropy scaling, in agreement with the earlier arguments.
This logarithmic growth of entropy versus β provides a tight upper bound in efficient thermal simulations [32]. This, together with the constant entanglement for gapped systems, suggests that the bond dimension D only needs to scale algebraically (constantly) as β increases for critical (non-critical) quantum chains, respectively.
Conversely, it directly follows from the above logarithmic scaling of the entanglement entropy that β needs to change significantly on a relative and not an absolute scale, in order to see sizeable effect on the entanglement entropy. This suggests for simulations that the numerical grid in β should be logarithmically discretized. In particular, as depicted in Fig. 2 by doubling β → 2β, one can therefore design an exponentially faster cooling procedure, in contrast to current standard simulation techniques which linearly evolve the full density matrix A single step in XTRG evolution by projecting MPO ρn (at β = 2 n τ ) to itself. Following common notation, tensor networks are graphically depicted by blocks (i.e., tensors) connected by lines which are to be contracted.
Here vertical lines indicate physical state spaces, whereas horizontal lines indicate virtual or bond state spaces. The exploitation of symmetries, quite generally, mandates directed lines, hence each line carries an arrow. [18][19][20][21][22][23][24] or the typical sampling states [28][29][30] in imaginary time. We note that essentially a similar, even though much more involved strategy was pursued in Refs. 26,27, and 37. Their approach was based on a dimensional reduction via a nested contraction of linear Trotter gates, followed by a variational optimization of coarse-graining transformations [38]. In contrast, our approach does not rely on Trotter gates, and hence is straightforwardly applicable to arbitrary Hamiltonians (1D and 2D). Overall, it represents an extremely simple yet also very efficient approach.
In this work, inspired by the logarithmic MPO entanglement entropy, we propose a one-way exponential tensor renormalization group (XTRG) scheme along imaginary time. Interestingly, this allows to draw parallels to the concept of energy scales in the Numerical Renormalization Group (NRG) [39][40][41]. There also, with every new iteration the energy scale is reduced by a factor, typically √ 2. Consequently, this also zooms into the low-energy regime in an exponential fashion, while dealing only with a very manageable linear number of iterations.
We benchmark our results with conventional linear evolution schemes. The results show that, by following the entanglement structure and exploiting the logarithmic temperature scale, one can obtain more accurate results with less cost. By implementing non-Abelian symmetries in the MPO, we can even simulate 2D clusters down to low temperatures with high precision, and investigate thermodynamics and related entanglement properties.
The model systems considered here are (anisotropic) spin-half Heisenberg XXZ models both, in 1D Heisenberg chains (∆ = 1) of length L, as well as in the 2D square lattice (∆ = 1 and 5) for systems of width W and length L, and thus with a total of N = W L sites, using open (OBC) as well as periodic (PBC) boundary conditions. We only include nearest neighbor couplings as indicated by the sum ·, · . For the purpose of benchmarking, we also consider the spin-half XY-chain with J z = 0, i.e., as this can be mapped to a fermionic tight-binding chain. The XXZ model in Eq. (2) possesses an U(1) symmetry, which restores a larger SU(2) symmetry when ∆ = 1. Symmetries, whether non-abelian or abelian, are fully exploited, throughout. We also set J := 1 as the unit of energy, unless specified otherwise. Furthermore we use units k B = = 1. The rest part of the paper is organized as follows. In Sec. II, we introduce the XTRG scheme with symmetries implemented, as well as an improved series-expansion thermal tensor network (SETTN) method [42] based on a pointwise Taylor expansion algorithm that also exploits the logarithmic temperature scale. The performances of these methods in the simulations of both 1D and 2D quantum many-body system, are presented and compared in Sec. III. In Sec. IV, the entanglement properties of MPO are investigated, where logarithmic entanglement entropies versus β in the Heisenberg chain and also square lattice models are discussed. XTRG is also employed to study the finite-temperature phase transition of 2D Heisenberg XXZ model, where we demonstrate that XTRG can accurately pinpoint the critical temperature.

II. SYMMETRIC THERMAL TENSOR NETWORKS IN LOGARITHMIC TEMPERATURE SCALE
By construction, a thermal density matrix ρ = e −βH is a scalar operator, and thus shares exactly the same symmetries as the Hamiltonian. For example, symmetries are preserved for Trotter-Suzuki type TTNs [22,23], where every local tensor (storing Boltzmann weights) is symmetric. Similarly, in the series-expansion TTNs, it is also clear that arbitrary H n 's have exactly the same symmetry as H (any unitary symmetry transformation that leaves H intact also leaves H n intact), and so does the resulting tensor network representation of ρ(β).
Concepts such as spontaneous symmetry breaking apply to individual low-energy (eigen)states, but not to a thermal state. Hence the full exploitation of all symmetries of the Hamiltonian, abelian and non-abelian alike, are very natural in XTRG. For finite systems with open boundary, in particular, XTRG shares the same benefits in efficiency as DMRG in quasi-1D. There is a notable difference, however: as long as the thermal correlation length is clearly smaller than the system size under consideration, local thermal properties in the center of the system can be regarded as in the thermodynamic limit.

A. Symmetric Matrix Product Operator
The explicit implementation of non-Abelian symmetries has been regarded as a standard technique in ground state DMRG simulations (T = 0) [43], which has many important applications including exploring quantum spin liquids in frustrated quantum magnets [44,45], and is also shown to be useful in METTS-type thermal simulations [46,47]. However, the implementations of non-Abelian symmetries in MPO for finite-temperature simulations are still absent. Here by virtue of the flexible and versatile QSpace framework [48], we implement non-Abelian SU(2), as well as Abelian U(1), symmetry in the MPO algorithm and thus realize a very efficient thermal renormalization group (RG) algorithm that can also be applied to 2D problems.
In our MPO-based thermal algorithm, we start by constructing an SU(2) invariant MPO representation of H. As this involves reduced matrix elements in the spirit of Wigner-Eckart theorem [48], we refer to this as the reduced MPO in contrast to the full MPO when not exploiting non-Abelian symmetries. By switching from a state-based to a multiplet-based description, we can reduce the overall bond dimension from D states to D * < D multiplets. The reduced MPO representation of H can be constructed by automata method [49][50][51] for 1D Hamiltonians and MPO sum-and-compress scheme [52] for more complicated 2D lattice models. For a Heisenberg chain with nearest-neighbor interactions, a full MPO requires D H = 5 bond states. As these correspond to two singlets and one triplet, i.e., 1 2 ⊕ 3 1 where d n specifies n multiplets of dimension d each, the reduced SU(2) invariant MPO only involves D * H = 3 multiplets. For the Heisenberg model on a 2D square lattice, we map a system of width W to a 1D snake-like chain with "long-range" interactions (up to distance 2W − 1). Then e.g., using OBC, the full MPO requires D H = 3W + 2 bond states, while the reduced SU(2) invariant MPO has a significantly more compact representation with only D * H = W + 2 multiplets (1 2 ⊕ 3 W ). More details on the symmetric MPO representation of total Hamiltonian can be found in Appendix B.
The computational cost in a tensor network algorithm typically scales with some power O(D m ) aside other factors concerning number of sites etc., where for the MPO structure of this paper we encounter m = [3, . . . , 6]. By exploiting non-abelian symmetries, the computational cost can be effectively reduced to O((D * ) m ) which leads to a gain in numerical efficiency by O D * D m . For a single SU(2) symmetry, it roughly holds, on average, D/D * ∼ 3 . . . 4 for spin-1/2 systems. Note also that multiplet dimensions are typically somewhat larger in thermal MPO as compared to matrix product ground states which renders us even greater numerical gain of symmetry implementation. The underlying reason for this is that an MPO has two physical indices associated with the same site. Therefore their direct product already also leads to an enlarged effective local spin. With this in mind, for the sake of readability, we will generally quote estimates in numerical efficiency in terms of D since after all, O(D m ) = O (D * ) m with the overall scale factor (D * /D) m 1 absorbed into the definition of O( * ).

B. Exponential Tensor Renormalization Group
For one-dimensional critical systems, the entanglement entropy in the MPO of a thermal state diverges only logarithmically in β. Therefore to see a sizeable effect in the properties of a thermal state, β must change significantly on a relative and not an absolute scale. E.g., a change β → aβ → a 2 β → . . . with some constant a > 1 will change the entanglement by linear increments. This strongly suggests to scale β on a logarithmic and not on a linear scale.
We can take fully advantage of this by a novel approach, which we refer to as exponential tensor renormalization group (XTRG), to simulate quantum many-body systems at finite temperatures, with high efficiency and accuracy. We start by preparing an MPO of the (unnormalized) thermal state ρ(τ ) = e −τ H at exponentially small τ , i.e., at very high temperature. Then we can proceed to cool down the system exponentially by multiplying the thermal state with itself, Feeding the last MPO iteratively into the next step, with τ n ≡ 2 n τ and therefore τ 0 ≡ τ , we obtain, This directly implies an exponential acceleration to reach low temperatures. Importantly, in the present XTRG scheme we can easily start from exponentially small τ . For example, for τ J = 10 −3 with J := 1 the largest local energy scale here given by the Heisenberg coupling strength, we can use an efficient series expansion scheme [cf. Sec. II C]. For τ as small as 10 −6 even simplest lowest-order linear expansion of e −τ H can suffice, which extremely simplifies initialization even for longerranged Hamiltonians which become cumbersome for Trotterlike decompositions, or for 2D Hamiltonians in the effective 1D-MPO setup. In the latter setup, with minor modifications, the MPO of the Hamiltonian already encodes the essential structure of the thermal state using the same bond dimension. A detailed comparison of different initialization strategies, including the Trotter-Suzuki decomposition, SETTN, and simple linear initialization for small τ 0 is provided in Appendix C.
Given an MPO representation for ρ(τ n = β/2), we can compute the (unnormalized) thermal state at temperature T = 1/β via ρ(β) = ρ(τ n ) † ρ(τ n ), i.e., by contracting ρ(τ n ) with its conjugate. This guarantees positivity of the thermal state ρ(β) even in the presence of truncation of the MPO for ρ(τ n ). Furthermore, we can also compute the partition function at β = 2τ n via Z(β) = Tr ρ(τ n ) † ρ(τ n ) , and thus gain another factor of two to reach lower temperatures. The latter can be simply obtained by computing the Frobenius norm squared of ρ(τ n ). Not incidentally, many of the features above are directly related with common procedures within the setup of a purified thermal state. [18-22, 24, 26, 27, 53].
In case the grid of inverse temperature values is too sparse, intermediate values can be easily obtained by shifting the initial value of a procedure that is entirely analogous to z-shifts within the NRG. In order to obtain a uniform logarithmic grid over n z shifts, one may simply choose Different "z-shifts" can be computed completely independently from each other and can therefore be efficiently parallelized. Truncation errors are still kept minimal by moving to large β as quickly as possible in an accurate manner. Alternatively, one can obtain intermediate values of β also by computing ρ n,n ≡ ρ n−1 · ρ n −1 for various n ≤ n.

C. Series Expansion Thermal Tensor Networks
Series-expansion thermal tensor network (SETTN) is a "continuous-time" RG approach for the accurate simulation of quantum lattice models at finite temperature [42]. By exploiting the series expansion of density matrix in Eq. (7a), SETTN is essentially free of discretization errors, making it distinct from previous Trotter-Suzuki type RG methods including TMRG [19,20], finite-temperature DMRG [21], LTRG [22][23][24], and METTS [28], etc. The efficient MPO representations of H n is the key for the algorithm to work, and both OBC and PBC chain systems can be equally well dealt with in SETTN (here, specifically, we simply use one longrange bond for the simulation of PBC). Being free of Trotter errors, SETTN has better controllable and uniformly higher accuracy, compared to conventional thermal RG methods.
To initialize ρ 0 for small τ 0 , a series expansion yields [cf. Fig. 3(a)] The required cutoff order of the expansion is N c ∼ N τ 0 , i.e., proportional to the total number of sites N . In practice, N c is determined automatically by only allowing a negligibly small expansion error (< 10 −15 ). Therefore for sufficiently small τ 0 , the initialization of ρ(τ 0 ) above is well-controlled and accurate, typically resulting in N c 10.
The high-temperature Maclaurin expansion in Eq. (7a) can be employed not only in the intialization stage, but also for simulating low-temperature thermal states, as shown in Ref. 42. Despite its competitive performance, this method still leaves room for further improvement. Since Eq. (7a) expands ρ around the infinitely high temperature, i.e., β = 0, the power series in H n involves large N c ∝ N β for large system size N and low temperature 1/β. The precision of SETTN is limited by the truncation in H n [see Sec. II D], which generally increases as n increases [42]. In this sense, a pointwise Taylor expansion can help reduce the expansion order N c and improve the accuracy, i.e., Equation (7b) expands the density operator around an arbitrary but fixed τ n . For generality, the initialization in Eq. (7a) may be viewed as iteration n = −1 having ρ −1 = I and τ −1 = 0. Now given the density operator ρ(τ 0 ) obtained by initialization via Eq. (7a), ρ(aτ 0 ) with a > 1 can be obtained via Taylor expansion around ρ(τ 0 ). In particular, this also hold for a = 2 which thus may serve as a complimentary scheme to the XTRG above. For example, alternative to the plain doubling scheme above, SETTN may be employed to cool down the system and obtain the MPO form of density operators at the inverse temperature grid τ n . Given ρ n , the MPO representation of ρ n+1 can be expanded as in Eq. (7b). Each term in the summation there can be obtained iteratively by projecting H onto (H n−1 ρ n ) and compressing the product. For the overall sum then [cf. Fig. 3(b)] we also employ variational optimization (see Appendix D) to finally arrive at the MPO representation of ρ n+1 ≡ ρ(τ n+1 ).
By repeating this procedure, we also can follow the XTRG protocol to cool down the system along the inverse temperature grid τ n = 2 n τ 0 . In contrast to the plain doubling scheme in XTRG, however, in case of SETTN the step size δτ ≡ τ n+1 − τ n can be chosen continuously. In this sense, SETTN is more flexible as it permits the flexible exploration of thermal properties in the immediate vicinity of temperature τ n with only modest cost.
Note that for this improved SETTN, as we will refer to it, using an exponentially increasing τ n series does not acquire exponential acceleration as XTRG does, since in the case of SETTN one still needs to perform projection and compression operations N c ∝ βN times. Nevertheless, from the point of view of SETTN, the exponential τ n series is computationally preferable to, say, a linear τ n series as expansion points, since the former can reduce expansion overhead and thus save computational time, in practice, without losing any accuracy (see Appendix E for a detailed comparison).
a XTRG cost is set as the time unit.

D. MPO compression and numerical cost
In SETTN we start with a reduced SU(2) invariant MPO for H. Then we iteratively apply the projections H onto (H k−1 ρ n ) to obtain H k ρ n , with Eq. (7a) represented by ρ n=−1 = I. These projections need to be combined with a compression algorithm to reduce numerical cost in a controlled manner. In the present context, however, truncation by discarded weight is dangerous since small weights for small τ 0 can affect the accuracy for large τ n . Hence we truncate by number of multiplets, throughout. For this, we introduce the control parameters D * n,k which stand for the maximum number of multiplets D * to be kept in the k-th iterative term when computing ρ n . For simplicity, we set this parameter constant, i.e., D * n ≡ D * n,k , which also stands for the bond dimension of the target state ρ n . Furthermore, we choose constant D * ≡ D * n>0 but, for the sake of the analysis, may use a different value for D * 0 for the initialization in Eq. (7a) if specified.
For an extremely small τ 0 (say, as small as 10 −4 to 10 −8 ), the initialization of ρ(τ 0 ) = e −τ0H can be simplified to lowest-order, i.e., linear expansion ρ(τ 0 ) 1 − τ 0 H. Having N c = 1 in Eq. (7a), the result shares the same bond-dimension D * 0 = D * H as H itself. In constrast, when expanding around finite τ n as in Eq. (7b), the bond dimension D * in ρ(τ n ) typically needs to grow significantly, and therefore is fixed to some specified D * D * 0 . The compression of the SETTN projections above can be achieved either by a singular value decomposition (SVD) technique quite similar to that in Ref. [42], apart from the fact that the MPOs here have SU(2) symmetry, or by a variational optimization which can greatly improve numerical efficiency (see Appendix D for more details on related MPO compression techniques). Within SETTN, the cost of either compression scheme scales like O D 3 . We tested both and found comparable numerical accuracy. Finally, we variationally add up the MPOs for H k ρ n with coefficients as in Eqs. (7) to obtain ρ n+1 (cf. Appendix D 2).
In contrast, XTRG projects ρ n onto itself in Eq. More explicitly, we summarize in Tab. I the time costs of the three algorithms involved in the current discussions. The numerical cost of SETTN for an entire run up to inverse tem-perature β scales as O βN 2 D 3 D H assuming N c ∝ βN for large β, N with β = τ n in Eqs. (7), whereas XTRG scales as O ln(β)N D 4 . We can thus estimate the relative run time of SETTN over XTRG as q S ≡ βN D ln β D H . For practical simulations as in Figs. 4 and 5, we find that XTRG calculations are faster than SETTN by more than one order of magnitude. In 1D critical systems, since the required bond dimension scales as D ∼ e S ∼ β λ (λ 1 for c = 1 CFTs, say, spin-1/2 Heisenberg chain, see Ref. 32). Thus, q S N ln β D H , with q S 1 for N large and β > 1 (in units of 1/J). Similarly, also Trotter-Suzuki type linear thermal RG methods, like the finite-temperature DMRG [21] and LTRG [22,24], with scaling O β τ N D 3 (last column of Tab. I with W = 1), are typically much slower by a factor q L 1 τ ln β 1 as compared to XTRG.
It is also revealing to compare the efficiency of XTRG with currently most efficient scheme in 2D systems, i.e., Trotter-Suzuki decomposition plus swap gates [30]. The numerical (time) cost of the latter scheme scales like O β τ N D 3 W , where the additional factor W stems from the number of required swap gates which is proportional to the width W . For 2D, however, typically D β. Therefore the relative cost of XTRG scales like q 2D ≡ τ D W β ln β, resulting in q 2D ∼ O(1) [e.g., with W = 8, τ = 0.05, β = 50, and D * ∼ 2000 (correspondingly D 8 × 10 3 ∼ 10 4 ) in SU(2) simulations, or D ∼ 2000 in U(1) calculations, one obtains q 2D ∼ 0.98]. Nevertheless, XTRG is still clearly advantageous over Trotter gates due to the far fewer truncation steps involved. Besides, XTRG can be simply and efficiently parallelized based on zshifts [cf. Eq. (6)]. For thermal simulations that are dominated by Trotter error and swap gates in LTRG schemes to reach low temperatures, and not necessarily by the truncation error due to entanglement growth, XTRG may be crucial to reach the lowest temperature scales e.g. in systems with more than one well separated physical energy scales. As a very interesting example, in the Kitaev honeycomb model the gauge flux excitation peak in the specific heat curve locates at very low temperature, to be referred to as T l , compared to the hightemperature peak at T h . To see the T l peak in the specific heat, one needs to cool down the system till extremely low temperatures (T l /T h 10 −2 , depending on the coupling constants and boundary conditions, [4]) where the exponential acceleration in XTRG can play a very important role. A similar scenario is observed in the 2D triangular Heisenberg lattice which we will discuss in more detail below. huge challenge to numerical simulations.

A. Thermodynamic quantities
Here we briefly summarize the thermodynamic quantities that will be computed and analyzed in detail below. An equilibrium thermal state is described by the partition function Z(β) ≡ Tr e −βH ≡ Tr [ρ(β)]. Typical interesting thermodynamic quantities, which constitute important tasks for the TTN algorithms to compute, including etc. The computation of the free energy f and energy density u, are straightforward, where H and ρ(β) are expressed as MPO, and the calculations amount to efficient contractions of tensor networks consisted of these MPOs. The linear derivatives in β for the specific heat as well as the energy density, however, are not very natural for XTRG, which obtains the thermal data on a uniform logarithmic β grid. Therefore it is more suitable to use ∂ ∂β = ∂ β∂(ln β) , i.e., instead. This is also more stable numerically for small temperatures since, the quotient of numerical differences is divided by T for the specific heat in Eq. (9b) and not by T 2 as in Eq. (8c), and is multiplied by T for the internal energy in Eq. (9a). This formula is used to compute the specific heat in Figs. 6, 7. In order to reduce numerical differential errors, independent calculations with slightly different initial τ values are run in parallel, e.g., using n z = 16 in Eq. (6b) which produces interleaved data points with δz = 1/16, i.e., δ ln β = δz ln 2 0.0433. Finally, we note that the LTRG approach adopted in this work, e.g., for the data in Figs. 4 and 5 below, has been streamlined with the remainder of the TTN procedures used in this work. It differs from the original LTRG algorithm in Refs. [22,24] in that it successively projects the MPO for ρ(τ 0 ) to the density operator ρ(β) to increase β linearly. Since the Trotter-Suzuki decomposition is not involved in the procedure, it is thus free of Trotter error.

B. Heisenberg chain
Firstly, we benchmark XTRG results with conventional linear evolution in Fig. 4, where an 18-site spin-1/2 Heisenberg (PBC) chain [cf. Eq. (2)] is calculated up to β 82. From Fig. 4(a) it is clear that the accuracy in the low-T regime gets continuously improved as D * increases in XTRG. Starting from a fixed τ = 0.01, XTRG reaches a precision as good as 10 −8 also at the lowest temperatures for D * = 500 as compared to exact diagonalization (ED) data. By keeping the same bond dimensions, LTRG and SETTN have almost the same accuracy in all temperature regimes. When compared to XTRG, they are of similar accuracies only at high temperatures (β 1), but are clearly less accurate in the low-T region where truncation errors dominate. These remarkable results suggest that as XTRG targets the low-T properties much faster than LTRG as well as SETTN, due to the (much) fewer truncation steps and its algorithmically much simpler setup, it also gains better results.
The MPO entanglement, as defined in Sec. IV and measured in the center of the chain, is plotted in Fig. 4(b), which offers a quantitative estimate of computational complexity. The entanglement data suggests that truncation errors start to develop for β 1 and the simulation errors stop increasing due to the convergence of entanglement for β 10. This is also clearly reflected in the overall error in physical quantities such as the free energy in Fig. 4

(a).
Besides the Heisenberg chain, we also benchmark XTRG, LTRG and SETTN for an XY chain [cf. Eq. (3)] with size L = 50, where analytical solutions are available (Appendix F). As shown in Fig. 5, again XTRG gets better results than LTRG and SETTN, and the accuracy in the low-temperature regime also improves continuously as we increase bond dimensions D. This simulation on longer XY chain again confirms that increasing β exponentially fast not only improves the efficiency but also gains in accuracy.
Besides the free energy, we also calculate the specific heat of a spin-1/2 Heisenberg chain of length L = 300, utilizing the XTRG algorithm with bond dimension up to D * = 250. As shown in Fig. 6, in the low-T region, the specific heat of the system shows a universal linear relation versus temperature, as indicated by the polynomial fitting (purple dashed line) with η 1. In addition, the fitted slope is also in perfect agreement with the well-known value 2/3 from CFT prediction [54], from which we extract the central charge c 1. On the other hand, in the high-T regime, the specific heat is also universal, and decays as 1/T 2 (a polynomial fit in the log-log scale, depicted by the yellow dashed line yields an exponent µ = −2.01). This exponent can be confirmed by a high-temperature expansion up to the second order, which approximates the energy as where Z 0 = Tr(I), with the high-T limit c V ∼ 1/T 2 .

C. Square lattice Heisenberg model
Symmetric TTN methods, including the XTRG and the SETTN, can be conveniently employed to calculate 2D systems, with minor adaptations. We map the 2D clusters into a 1D snake shape, and prepare the MPO representation of this Hamiltonian (with "long-range" interactions) as elaborated in 1, a universal linear behavior versus T is observed, i.e., cV = αT η with fitted exponent η = 0.996 and slope α 2/3, in the regime T ≤ 0.025. Exploiting the fact α = πc 3v (v = π/2 for spin-1/2 Heisenberg chain), we extract the central charge c 1. For large T , the specific heat shows a universal 1/T 2 temperature dependence (the fit shown was performed for T > 15). The inset shows the same data on a linear vertical scale.
Appendix B. Other than that, one follows exactly the same line as in 1D simulations and can represent the density matrix of the 2D systems accurately in terms of MPO.
Here we perform calculations on 2D clusters and benchmark the calculations with ED for small (OBC) systems (4×4) in Figs. 7(a-c) and quantum Monte Carlo (QMC) for larger systems (16 × 5) in Figs. 7(d-f). With non-Abelian symmetries implemented in the highly efficient XTRG algorithms, we obtain high quality data till quite low temperatures, which was not accessible before by other thermal RG algorithms.
In Figs. 7(a-c), we show the free energy, energy and specific heat results of a 4×4 square lattice Heisenberg (SLH) model. Very nice agreement between XTRG and ED data is observed in all three plots. As seen in the inset of Fig. 7(a), the relative accuracy is quite high, i.e., 10 −4 (10 −5 ) for D * = 200 (400) at low temperatures (T ≤ 0.05). The energy density shown in Fig. 7(b) is obtained by taking derivatives of interleaved XTRG free energy data [cf. Eq. (9a)]. The error in the energy density is small even down to T ≤ 0.05, as seen in the inset of Fig. 7(b) by zooming in the low-T region. XTRG data differs from the ED results in the fourth digit for D * = 200, and is bounded by numerical differentiation error for D * = 400. In Fig. 7(c), we show our results of the specific heat, which was calculated by taking derivatives of energy data as in Eq. (9b). Inset plots c V on a log-log scale, from which we again observe an algebraic behavior (1/T 2 ) at high temperatures. In the low-T region, it shows a very rapid (exponential) decay versus the temperature. This can even be observed in the dashed region of Fig. 7(c), although there the XTRG data departs from ED results due to lack of accuracy.
For a 16×5 SLH which is far beyond the scope of ED cal- culations, we compare our XTRG results to those of quantum Wang-Landau (QWL) simulations [55] in Figs. 7(d,e,f). We run the calculation down to T = 0.025. For the smallest temperature T = 0.1 for which with have well-converged QWL reference data at comparable numerical cost, the error in the D * = 600 data for the free energy is ∼ 2 × 10 −4 . Since the truncation errors is generally larger on 16 × 5 lattice, we extrapolate the free energy in Fig. 7(d) to 1/D * → 0 and observe a perfect agreementwith QWL data, with the error further reduced by about an order of magnitude. Besides the free energy, in Figs. 7(e-f) we also show that the energy density u as well as specific heat c V all have very good accuracy. In the high-T region, c V again shows a 1/T 2 relation, as shown in the inset of Fig. 7(f) and as already discussed with Fig. 6. In the very low temperatures [dashed region of Fig. 7(f)], though, our accuracy becomes somewhat limited as the derivative to obtain c V develops minor wiggles at the lowest temperatures, seen on the log-scale in the inset. However, the upturn does not appear to be due to numerical inaccuracies, but appears to be physical, in the sense that, similar to Fig. 7(c), it is a precursor before the finite-size spectral gap sets in.
In Figs. 7(b, e), for comparison we also included METTS data exploiting U (1) symmetry only [30]. The METTS results also show good agreement with our other methods in both cases, apart from the fact that the METTS energy data is not strictly variational, i.e., could be even lower than the (quasi) exact value. As for the 16 × 5 plot in Fig. 7(e), note that T = 0.1 is currently the typical lowest temperatures that 2D METTS simulations can reach [30] at computational resources that are comparable to XTRG. However, the current 2D METTS involves many swap gates and needs at least a few hundreds of samples. In contrast, our XTRG method, with SU(2) symmetry implemented, is much more efficient and can reach much lower temperatures with great accuracy in 2D.
of TLH were hindered by lack of sufficiently powerful numerical approaches.
Besides non-frustated systems such as the SLH, XTRG can also be applied to frustrated magnets like TLH. In Fig. 8, we present exemplary XTRG results of TLH on width W = 4 systems [see inset of Fig. 8(d)], including open strips (OS) and Y cylinders (YC). To be specific, we consider YC4 ≡ YC(W=4) geometry with various aspect ratio L/W = 1 ∼ 2.5, as well as a 4 × 4 OS for comparison.The YC geometry we adopt in practice is shown in the inset of Fig. 8(d), and also in Figs. 9(f-i) after a proper transformation to restore the corresponding triangular lattice geometry. Our realization is essentially equivalent to conventional YC in previous DMRG studies [73] at T = 0. We cool down the system from high temperatures to as low as T /J = 0.03, and D * = 400 multiplets are kept in the following calculations, to insure convergence vs. bond dimensions.
Although YC4 shows clear finite-size effects in the ground state at T = 0, nevertheless, we can observe how correlations gradually develop as we lower temperatures, with relevant finite-T physics in the thermodynamic limit at larger T down to temperatures where finite-size effects set in. In Fig. 8(a) we plot the specific heat of a 4 × 4 OS and various YC4 lattices. By comparing Fig. 8(a) to corresponding c V in Fig. 7(c) of the SLH, we observe quite distinct features in the TLH case. Already on the 4 × 4 lattice, either OS or YC, the specific heat c V exhibits a peak, whose location will be denoted by T l , and a shoulder structure at a higher temperature T h ∼ 0.5 ∼ J. For YC4 lattices with increasing length L = 6, 8, 10, these features in c V develop into a pronounced two-peak structure where the plateau-like feature for L = 4 develops into a broad peak located around the same stable value of T h . Note that T h is also about the same temperature scale where the data from high-temperature expansion (HTSE) shows a round peak [56]. In stark contrast, the peak at T l moves towards lower temperatures. Its value scales like T l ∼ 1/L as seen in inset of Fig. 8(a) and therefore is clearly linked to the finite system size. Overall, the stark qualitative difference between the data for TLH in Fig. 8(a) as compared to SLH in Fig. 7(c) may be ascribed to finite-temperature effects of magnetic frustration.
To gain a better understanding, we analyze the static structure factor where r 0j ≡ r j − r 0 with r j the lattice location of site j.
In practical calculations, site 0 is fixed in the system center, while j runs over the whole lattice. Therefore, by inversion symmetry of the TLH, S(q) is a real number. Note that also the structure factor above can be conveniently and efficiently obtained via the expectation value of single q-dependent MPO of bond dimension D * S = 2. The static structure factor for the THL is analyzed at specific values of q in Fig. 8(b). However, before discussing these in detail, let us look at the static structure factor over the entire Brillouin zone, as presented in Fig. 9 for a YC4 system at L = 8 along with the bond energy texture, at various temperatures (including T = 0, obtained by DMRG). In agreement with the previous discussion, Figure 9(a) indicates the existence of three regions, separated by the two temperature scales T l and T h as introduced with Fig. 8(a). For T T h , the system is in a "spin gas" paramagnetic phase, with featureless structure factor and bond texture in Figs. 9(b) and (f), respectively. In Fig. 9(c) we can observe that the intensities of representing the development of strong 120 • -ordering correlations.
In addition to that, there also emerges an extended region with considerable intensity around q = M ≡ ±2π(0, 1 , near which roton-like excitations were reported both theorectically [57,65,66,68] and experimentally [74]. The (anomalous) enhancement of S(q) around these M points for temperatures T ∼ T h , which can be more clearly observed still in Fig. 8  with roton excitations. By association with T h , these are also linked to the round peak/shoulder in the specific heat c V . The corresponding bond texture in Fig. 9(g) reveals much stronger spin-spin correlation than those in paramagnetic phase, while it still maintains a uniform pattern that respects the TLH lattice symmetry, like a liquid state. Thus we dub the finite T region between T l and T h anomalous quantum liquid phase.
As T is lowered further down to less than T l ∼ 0.14, the system in Fig. 9 undergoes a rapid crossover. Due to the finite YC circumference, tightly bound RVB rings of length W form as can be clearly observed in the bond textures in Fig. 9(h,i). There emerges a vertical strip pattern around the cylindrical circumference at T T l [ Fig. 9(h)], which becomes strongly dominating for lower temperatures [ Fig. 9(i)]. Note that the low-temperature (T /J = 0.03) bond texture [left half of Fig. 9(h)] is in perfect agreement with T = 0 DMRG results [right half in Fig. 9(h), glued together at dashed center line]. The agreement between bond energies is better than 0.1% for most bonds, indicating the XTRG already effectively reached the T = 0 regime. This is further reflected in the average energy per site e 0 = −0.53027(1), for given YC4 system at L = 8 and T /J = 0.03, which is in perfect agreement with, i.e. just above the DMRG ground state energy e 0 = −0.53034 [with D = 4000 U(1) states retained]. Similarly, also the structure factor is in very good agreement, e.g., with S(M ) = 1.82 at T = 0 DMRG to be compared to S(M ) = 1.79 at T /J = 0.03 XTRG.
As a consequence of the tightly bound rings around the circumference of the YC4 cylinders, the correlations become significantly enhanced along a line connecting two of the three initially equivalent M points in the Brillouin zone [dashed line in Fig. 9(d)]. This strongly competes the triangular magnetization associated with the K points [ Fig. 9(d)], such that the structure factor S(M ) at these points eventually dominates over S(K) [ Fig. 9(e), and also Fig. 8(b)]. This gives rise to the low-T peak of the specific heat curve in Fig. 8(a). From Fig. 8(b), it can be observed that almost exactly at T l , the two correlations cross and S(M ) becomes larger than S(K) for T < T l . Here the symmetry breaking across the initially equivalent M points is due to the finite-size cylindrical structure.
The results of magnetic susceptibility χ of various YC4 lattices are shown in Fig. 8(c). It can be observed that W = 4 data can produce results in agreement with the HTSE [56], bold diagrammatic Monte Carlo (BDMC) [58], as well as Pade approximation results [57], down to T ∼ 0.4 (for YC 10 × 4), despite a very limited circumference. This observation is quite surprising, which suggests a very small finite-size effect and correlation length ξ ( 1) for T /J ≥ 0.4, in accordance to previous studies [56].
In Fig. 8(d), we plot the bipartite entanglement entropy S E in the purified state. Since S E differs from bond to bond, we show there the maximal value amongst all bonds in the effective 1D chain structure adopted in the calculation [inset in Fig. 8(d)]. It can be clearly seen in Fig. 8(d) that this low-temperature entanglement reaches a relatively small value around S E = 2.16. DMRG simulation shows that the corresponding entanglement S E 1.1 which, as expected, is roughly one-half of the entanglement in low-T mixed state.
The tightly bound stripes in the low-T regime in given YC4 system lead to a strongly reduced entanglement entropy when cutting the system exactly in between the stripes, i.e. vertically the inset in Fig. 8(d), where the entanglement is reduced e.g. down to S E ∼ 0.9 at T = 0.03 . Therefore at low temperatures T < T l the system is close to a direct product of tightly bound uniform 1D rings around the circumference of the YC4 cylinder. These four-site rings favor a well-known RVB ground state |RVB = 1 constitutes a valence bond, and the lowest excited state is separated by a significant energy gap ∆ ∼ J. These tightly bound rings are disrupted when opening the boundary. Therefore in the OS geometry, e.g., in the 4×4 data in Fig. 8(b), S(K) remains strong, and never crosses with S(M ) at M = ±2π(0, 1 √ 3 ) as T is lowered. However, for the 4 × 4 OS, the dominant weight of the structure factor turns out to be located at one out of three types of Néel antiferromagnetic order, i.e., S(M ) with M = ±(π, −π/ √ 3) [data not shown in Fig. 8]. The entanglement also increases monotonically on the OS lattice, until it saturates due to finite system size [ Fig. 8(d)].
Recently, BDMC has been employed to explore thermal properties of TLH, and it was found that the "extrapolated" ground state is disordered, via a particular quantum-toclassical mapping based on their best thermal data (with temperatures down to T /J 0.375) [58]. The anomalous thermodynamic behavior is in contradiction to the ground-state 120 • ordering, which reveals that the true low temperature regime has not been reached. To fully understand the finitetemperature anomaly and to resolve the above apparent contradiction, a more extensive XTRG survey of the TLH e.g. on wider cylinders down to low temperatures is required. This is beyond the scope of this paper, and thus will be reported elsewhere [75]. For the purpose of this paper, nevertheless, we have demonstrated that XTRG provides a highly competitive approach that allows one to tackle complex and challenging problems.

A. Thermal Entanglement Renormalization Group Flow
The entanglement measure in a thermal state is more complicated as compared to a ground state due to the interplay of classical correlation and quantum entanglement. Among various definitions, we take a very natural and most relevant measure of the entanglement in practise, i.e., entanglement in the normalized "superstate" [76] which vectorizes the MPO for e −βH/2 . In other words, the MPO is simply transformed into a matrix product state (MPS) with doubled local state spaces. Then the partition function is equivalent to the overlap of the unnormalized superstate Z(β) = Ψ (β)|Ψ(β) , whereas Ψ(β)|Ψ(β) = 1.
Note that this definition, is a specific (and most natural) choice of purification [21,77], which in some other context is also called the thermofield double (TFD) state |Ψ(β) = En |n,n , where E n is the eigen-energy of eigenstate |n and its duplicate |n in the auxiliary state space [53,[78][79][80]. Here, for simplicity of notation, by the entanglement entropy or the entanglement spectrum (ES) of the MPO or the thermal state, we refer to precisely these quantities obtained from the underlying TFD, or equivalently, the purified and normalized state in Eq. (11). Specifically, the entanglement spectrum is given by the eigenspectrum of H ES ≡ − ln R where R is the 'super'-density-matrix of the purified thermal state |Ψ(β) as in Eq. (11). Note that the MPO entanglement analyzed here is not directly related to the entanglement of purification, which is defined as the minimal value amongst various purification schemes [81,82]. Nevertheless, through the optimal truncation via orthogonal state spaces in the XTRG (and also LTRG), one is simultaneously optimizing the super-state overlap (i.e., partition function), as well as this MPO (TFD) entanglement. Therefore, this MPO entanglement, as well as some other measures, like the mutual information, quantifies the resources required to perform efficient thermal simulations and, thus, have attracted recent interest [32,83,84].
We start by analyzing the entanglement spectra of the thermal state for a spin-1/2 Heisenberg chain, from which we can also compute its entanglement block entropy S E . By lowering the temperatures, one generates a RG flow that directly reflects different physical regimes of the system at various temperatures, i.e., energy scales. In Fig. 10, we show the entanglement RG flow over a very wide range of temperatures. We can vary β over 7 orders of magnitude, which thus reaches far beyond Trotter-Suzuki type calculations. The RG flow reveals three distinct regimes, demarcated by vertical dashed lines in Fig. 10: (i) a low entanglement region β 1, (ii) an intermediate region 1 β 100 where entanglement rises quickly, and (iii) the saturation region for β 100 where the ES flows to a fixed point, either converging to the ground state of a physically gapped state in the thermodynamic limit, or resolving the gap of finite size level spacing.
When approaching the low-energy fixed point ES, lines systematically merge into groups with larger degeneracy as seen in the inset of Fig. 10. Given that the entropy already clearly converges to a finite value at the lowest temperatures, this suggests that the low-energy fixed-point spectrum must be related to the tensor product space of two copies of the ground state (bra and ket) which naturally results in systematically enlarged degeneracies. For the remainder of this section, we focus on the entanglement entropy both in 1D chain and 2D lattice models over a wide range of temperature scales. Interesting logarithmic behaviors are observed, which intimately relate to (gapless) lowenergy excitations, yet also suggest efficient computational complexity of thermal simulations for the specific model systems considered.

B. Universal entanglement behavior in (1+1)D conformal thermal states
In Fig. 11 we plot the entanglement entropy in a spin-1/2 Heisenberg chain. By simulating twice the length (L = 200) as in Fig. 10, we can observe a logarithmic divergence in the low temperature region (L β 1) The logarithmic entropy was already observed in the past and related to the computational complexity of finite temperature simulations [16,31]. More recently, this was further analyzed by numerical simulations on Renyi entropy and also conformal field theory (CFT) analysis [32]. Notably, the finite-temperature entanglement calculation in Fig. 11 provides a convenient and accurate way to extract the central charge c of CFT: without going into ground states calculation at T = 0 [85][86][87], one can fit the MPO entanglement at finite temperatures. In Fig. 11, we observe that, by fitting L = 200 data, the estimate of central charge c is already very accurate (c = 0.999 1). For this we fit the data to the CFT prediction S E = c 3 ln β + const [32]. Importantly, by having the system length sufficiently large, the physics of the thermal state in the center of the system is effectively short-ranged by a thermal correlation length. In this sense, the simulation of the central charge in Fig. 11 does not yet see the finite open boundary condition. This is in stark contrast to the evaluation of central charge using ground state properties.
Universal features of the entanglement property appear also at large temperatures. As seen in the left inset of Fig. 11, the MPO entanglement shows a power-law behavior for β 1. The slope γ on the log-log plot at very high temperatures is analyzed in the right inset in Fig. 11, which suggests a powerlaw exponent γ ≈ 2 for β 1. The growth in the entanglement S E , however, slows down strongly for S E 1, i.e., β 1, where γ drops significantly below 1 as seen in the right inset of Fig. 11. The entanglement behaviors at high temperatures can be understood from a lowest, i.e., first-order expansion of density operator, ρ(τ ) = I − τ H. The singular value spectrum of this MPO is given by the vector s = [1, ατ ], with α another numerical vector. The resulting normalized "density matrix" of the supervector in Eq. (11) has eigenvalues r i = s 2 i / i s 2 i with lowest-order thermal contributions ∝ τ 2 . With the von Neumann entropy S E = − i r i ln(r i ) and one obtains that (12b) For simplicity, one may consider a single value for the vector α, resulting in the two normalized weights with von Neumann entropy A subsequent one-parameter fitting of S(τ ; α) with respect to α to the actual entanglement entropy S E for τ = β 1 nicely reproduces the high-T entanglement data, as shown in the right inset in Fig. 11 for α = 0.05. The slope γ decreases monotonically, starting from γ = 2 and undergoing a sharp decrease around β ∼ 1. Note, however, that the convergence towards the power-law exponent of γ = 2 for small τ is extremely slow, as also clearly supported by the simple asymptotic analysis above. For example, for τ = 10 −3 , one only has γ 1.94.
Nevertheless, it follows from the generality of the above asymptotic argument, that the exponent γ = 2 for infinitesimal τ is universal. It should hold for any Hamiltonian, and therefore, in particular, also in arbitrary dimensions. Furthermore, given that by construction, the exponent γ = 2 only holds for β 1 where S E 1, the growth of the entropy of the MPO may be considered sublinear in this regime, in the sense that the entropy grows slower than linear for infinitesimal τ , having lim β→0 + dS E dβ = 0.

C. Logarithmic entanglement in thermal states of 2D Heisenberg model
The low temperatures entanglement of the thermal state saturates for gapped quantum chains and grows only polynomially for critical ones. This directly implies excellent numerical efficiency in 1D quantum systems, since the required bond dimension grows at most polynomially with inverse temperature β [16,31], rather than exponentially as originally estimated [17]. We take this as a motivation to explore the MPO entanglement of the Heisenberg magnet on the square lattice, and take it also as an indicator of computational complexity for the latter.
In Fig. 12 we plot the entanglement property versus inverse temperature β, for the SLH of various system sizes, ranging from width W = 2 to W = 6. As shown in the inset of Fig. 12, in the high temperature regime one still recovers the universal power-law γ 2, c.f. Eq. (12). Specifically, the data is very well fitted by the function in Eq. (13) with exactly the same parameter α = 0.05 as in Fig. 11 for the 1D quantum chain. For the inset of Fig. 12, since the entanglement entropy satisfies area law, i.e., S E ∝ W , we divide the entropy S E by the width W . For high temperatures β 1, this collapses the data for different system widths on top of each other, indeed, demonstrating universal area law in this regime. For intermediate and low temperatures, β > 1, deviations from the strict scaling collapse of the area law can be observed.
In the main panel of Fig. 12, the entropy S E changes gradually into a logarithmic divergence versus ln β as W increases for both, even and odd widths. This suggests that simulations are also efficient with increasing β in the 2D setting, while bearing in mind an additive constant term to the entropy that is proportional to the width (note that in order to satisfy area law, the entropy data for 2 ≤ β ≤ 20 is separated roughly by equal vertical offsets when incrementing the width for W ≥ 3). Since the calculations of the entropy in the system center are not fully converged for the wider systems, we also extrapolate the width W = 5 and 6 systems in 1/D * → 0. This actually further reinforces the regime of logarithmic increase of S E vs. ln β for 2 ≤ β ≤ 20.
A similar additive logarithmic scaling of the entanglement entropy (∝ ln W ) in 2D Heisenberg model has been also found numerically via QMC calculations in the ground state [80,[88][89][90][91][92]. The coefficient of logarithmic correction was argued to be universal [93] (proportional to the number of Goldstone modes in the system). In the 2D Heisenberg model here, according ED and DMRG studies of the tower of states (ToS) in the energy and the entanglement spectra [94,95], respectively, the relevant low-energy ToS has characteristic level spacing that scales as 1/N with N = W L the total number of sites. This is in contrast to spin wave excitations (Goldstone modes) which have characteristic level spacing that scales with inverse linear system size, i.e., 1/L. These 1/N ToS excitations are responsible for the logarithmic entanglement at T = 0, and possibly also relate to the ln β scaling of the entropy observed in the present study. In Fig. 12, we restrict the length to be as small as L = 10, which suggests that the magnon excitations are gapped out at temperatures as low as T = 1/30 -1/10. Therefore the relevant energy scale are likely only ToS modes with energy level spacing ∝ 1/N ∼ 1/50, which is smaller than the temperatures in the regime where logarithmic scaling S E ∼ ln β is observed.
The logarithmic growth in entanglement in Fig. 12, as well as the low-T specific heat behavior, are quite remarkable. As we will show shortly, they differ qualitatively from the anisotropic case, ∆ = 1. In the isotropic case, the entanglement curves do not show any sign of singularity at any finite T , which suggests the absence of phase transition at T = 0. The low-temperature specific heat c V gets significantly enhanced when the width W is increased from 4 to 5, but it still grows monotonically with increasing β, e.g., with no sign of a singularity at finite β. This is, of course, in complete agreement with the celebrated finite-temperature Mermin-Wagner theorem [96].

D. XTRG and thermal phase transitions in 2D
Lastly, we analyze the 2D anisotropic XXZ model on the square lattice, i.e., Eq. (2) with J = −1, and ∆ = 1. There exists a finite-temperature phase transition at the critical temperature T c towards a gapped low-energy ferromagnetically ordered phase. Note that while Trotter-like methods have to necessarily work straight through a thermal phase transition point, which becomes singularly hard already with the first one encountered, our XTRG can jump across phase transition points [26] by reducing temperatures by a factor of 2!
In the following, we show that XTRG can be employed to simulate such model with nonzero T c with high accuracy. In particular, we determine the phase transition point using various quantities including the block entanglement entropy, specific heat, as well as the Binder ratio in spin fluctuations. We choose the same model parameters as in Ref. 30, i.e., ∆ = 5, where the critical temperature was estimated as T * c /∆ = 0.56 ± 0.01. This was in agreement with the exact result T c /∆ = 0.56 in Ref. [97] from QMC for much larger system sizes, which we will also take as reference for our data below.
Our results are presented in Fig. 13, where we also make comparison to QMC data explicitly generated by the ALPS looper code [55]. Figure 13(a) shows the landscape of the block entanglement entropy across each bond of the MPO density matrices with decreasing temperature. Because the low-temperatures phase is gapped, prominent peaks are present around the critical temperature T c = 0.56, with only weak dependence on the snake-like serial ordering of our 2D system, otherwise. Figure 13(b) shows a top view of the same data where the temperatures T S c of the maxima of the entropy data S E (dashed dotted line) agrees very well with the exact QMC value T c (vertical dashed line).
To see this more clearly, we also take a cut in the center of the system [while shifted by half a column to avoid local minimal, see caption of Fig. 13(c)]. In Fig. 13(c) we show S E vs. temperature T for various system widths W . From this we observe that the peak positions T S c in S E vs. T quickly approaches T c as W increases. Specifically, for W = 8 we already have T S c 0.552 with an error |T S c −T c | < 0.01. Similar calculations are also performed using cylindrical boundary conditions, as show in Appendix G [Fig. A.6], where we also see a quite accurate agreement with T c for rather small W .
We also resort to other more standard thermal quantities including specific heat c V [Fig. 13(c)] and Binder ratio [ Fig. 13(d)] to pinpoint the critical temperature. Our results for c V by XTRG are in perfect agreement with QMC simulations, as shown in Fig. 13(c). The specific heat also exhibits a peak near T c . As is well known from numerous QMC simulations, however, similar to other physical observables such as the spin susceptibility, the specific heat suffers significant finite-size corrections. So this often only provides a first rough estimate for T c in practical numerics, where even the finitesize extrapolation to the thermodynamic limit often has larger error bars still, as is also the case here [cf. inset to Fig. 13(c)].
The lowest order, finite size corrections in standard thermal quantities, however, can be eliminated by taking ratios of expectation values. In practice, a common way to pinpoint T c more precisely, is the Binder ratio [98,99] where S z tot = i S z i is the total magnetization in z-direction, with S z tot = 0. We calculate U 4 with both, XTRG and QMC the looper code, and find perfect agreement. Note that the Binder ratio U 4 can be conveniently calculated using an exact MPO representation for the total magnetization (operator) S z tot = i S z i of bond dimension D = 2. By taking product of S z tot and performing compression (without any essential truncations), one can obtain compact representations of (S z tot ) 2 and (S z tot ) 4 with D = 3 and D = 5 MPOs, respectively. With these MPOs one can evaluate the finitetemperature expectation values . β of the two operators required for the Binder ratio in Eq. (14).
From the data as in Fig. 13(d), we collect the crossing points T X c for two consecutive curves for widths W and W + 1. The data is summarized in the inset to Fig. 13(d). By considering the data for the largest system to be the most accurate given the aspect ratio L/W = 2, we obtain T X c 0.554 from both the XTRG and QMC data. By considering the minor trend still with increasing system size, we find agreement between this T X c and T c to within 1% which thus serves as a very good estimate for the thermodynamic limit.
From the above we conclude that XTRG can be used to determine T c of thermal phase transition very accurately. Fur- thermore, the maximum in the entanglement entropy S E itself can already provide a good estimate for the critical temperature. Here, in particular, it clearly outperforms conventional thermal quantities like the specific heat c V , as seen in the upper inset in Fig. 13(d). Nevertheless, since S E grows linearly with W for T > T c but stays constant for T < T c in given system (the system nearly becomes a product state for T T c ), the peak around T S c is rather round. This might lead to systematic offsets in T S c in the thermodynamic limit, but needs further studies.

V. CONCLUSIONS AND OUTLOOK
Inspired by the logarithmic growth of entanglement of purified (thermofield double) states, we propose an exponential speed up of thermal simulation. This thermal tensor network algorithm employs an MPO form of the density operator and proceeds via doubling of the density matrix ρ(β) along the imaginary time evolution. We show that this exponential tensor renormalization group (XTRG) method gains both accuracy and efficiency, in thermal simulations of the Heisenberg models. We also implement this idea of logarithmic temperature setup in a pointwise series-expansion thermal tensor network (SETTN) algorithm. Also there we get more efficient and accurate results than previous Maclaurin SETTN [42].
We apply XTRG and SETTN to efficiently simulate thermal states of 1D and 2D Heisenberg spin models, obtain accurately the thermodynamic quantities including free energy, energy, and specific heat, and study their low-and hightemperature behaviors. We have also investigated the temperature dependence of entanglement properties S E in the MPO, and observed logarithmic entropies S E a + b ln β with constants a and b at low-temperatures not only in gapless quantum chains, but also in the SLH at fixed system size due to gapless ToS modes.
We applied XTRG to a 2D Heisenberg models with a thermal phase transition , as well as the TLH with spin frustration. The results demonstrate the efficiency of the present algorithm, which is capable to show rich finite-T physics, especially for those system with spin frustration. It will be interesting to see how useful the present algorithm can be in exploring more 2D challenging systems, such as the kagome and the J 1 -J 2 Heisenberg models on the square lattice, as well as interacting fermionic systems.
With XTRG, however, we are not limited to high temperatures but can simulate down to much lower temperatures than previously anticipated [30]. The present MPO algorithms may be improved in several directions, including combining them with METTS samplings at low T , or linked-cluster expansion to reduce finite-size effects, etc. Moreover, XTRG may also be straightforwardly combined with efforts to reduce block entanglement entropy further by operating with disentanglers on the auxiliary state space in the purified scheme [81,82], all of which certainly deserves further exploration.

ACKNOWLEDGMENTS
The authors thank Benedikt Bruognolo for nicely providing the METTS data and for very constructive discussions. The block entropy in the center ∼ L/2 of an individual low-energy pure state |s scales like S E (|s ) c 6 log +const [86,87]. The MPO block entropy for the outer productρ s ≡ |s s| then acquires a factor of 2, i.e. S E (ρ s ) c 3 log + const. Now by going to a thermal stateρ = s ρ s |s s| with weights s ρ s = 1, this scaling does not change. More precisely, due to the subadditivity of entanglement entropy, one obtains an upper bound S E (ρ) c 3 log + const. More explicitly, the block entropy of the thermal density matrix changes most for the worst case that its spectrum is altered from having used s ρ s = 1, as well as i (s) i = 1 for all s. Here S({ρ s }) is the entropy of the weight distribution ρ s . Now, Eq. (A1) provides an upper bound. As argued with Fig. 1, in order to sample thermal averages, one requires that the system size must be larger than the thermal correlation length. For physical systems whose finite size spectra have a low-energy level spacing that scales like inverse system size, this is achieved by choosing L = aβ with a 1 sufficiently large but constant. Therefore one can substitute = L/2 → β in Eq. (A1), resulting in the overall block-entropy of the thermal state The block-entropy of the thermal state thus saturates, since the system length is effectively cut off by the thermal correlation length ξ, resulting in ∼ min(L, ξ) in the thermodynamic limit L → ∞. See also [32,33] for a more rigorous derivation based on CFT arguments. Furthermore note, that even Eq. (A2) can be still considered an upper estimate, since in the purification scheme, the thermal state allows one to minimize block entanglement entropy by disentangling operations on the auxilliary state space [81,82].

Appendix B: Symmetry invariant matrix product operator for Hamiltonians
In this Appendix, we discuss our approach to the implementation of both, abelian as well as non-abelian symmetries into the MPO representation of a given Hamiltonian. Conceptually, non-abelian symmetries proceed the same way as abelian symmetries, as we explain below. The actual implementation is based on the framework of the tensor libary QSpace [48] that can deal with abelian and non-abelian symmetries such as SU(N ) or the symplectic symmetry Sp(2N ) on a generic footing.
In order to emphasize the generality of the argument, we will frequently use the notation q for a label of a generic irreducible representation (irep) of a given symmetry. Here, specifically, it may either stand for the spin-projection S z or spin S label in the case of a U(1) [an SU(2)] symmetry, respectively. For this reason, we also refer to q = 0 only as the scalar representation even if for U(1) symmetry all symmetries multiplets are actually one-dimensional and in that sense, scalars. Examples for scalar operators are the full Hamiltonian, as well as all of its terms in its sum including local 1-site terms. Vacuum states transform like a scalar multiplet.
The dual representation q * of some given irep q is defined by the unique representation that allows to form a scalar, i.e., with Clebsch-Gordan coefficients (CGCs) (q, q * ; 0). These CGCs when properly normalized define a unitary matrix U [q] . This will be referred to as 1j-symbol by analogy e.g., to 3j symbols for the SU(2) spin symmetry, with the difference, that here only a single irep label is concerned. While SU(2) symmetry is self-dual, i.e., q * = q, abelian U(1) symmetries are not, since one has q * = −q such that q + q * = q properly adds up to the scalar representation q = 0. For self-dual symmetries, such as SU(2), however, one must be careful in that U [q] , when written simply reduced to a unitary matrix of rank-2, it becomes indistinguishable from (U [q] ) † = U [q * ] which, however, may differ by a sign, e.g., for half-integer spins in the case of SU (2). Importantly, 1j-symbols allow to revert arrows in lines in a tensor network by inserting I = U [q] † U [q] .

Automata approach
Firstly, we briefly recapitulate the automata approach for constructing MPOs of the Hamiltonian [49][50][51]. Consider, for example, the quantum Ising chain H = i S x i S x i+1 − hS z i , which provides a simple example of a Hamiltonian with a single nearest-neighbor interaction term together with a local term (here with magnetic field strength h). We need to compute and store the matrix elements of the spin operators {S x ,S z }. Together with the identity operator, I, these form a basis of local operators that enter a rank-4 tensor T α,α ,σ,σ , where by rank we refer to the number of indices (or legs in a graphical depiction) of a given tensor. The tensor T is the elementary local tensor of the MPO, with σ the local state space of a given site, and α the virtual bond states that tie together the MPO. The tensor T has the same form for every site due to the translational invariance of the Hamiltonian [assuming open boundary condition, the open virtual indices of the T -tensors for the first and last site are contracted ("capped") with a start and a stop state, respectively; see below].
To be concrete, each tensor T contains D 2 H d 2 matrix elements, where d is the dimension of the local state space σ, and D H is the bond dimension of the virtual state space α. Every matrix element of T in the indices (α, α ) is linked to a local operator with matrix indices (σ, σ ). It is therefore natural to group the relevant local operators into an (orthogonal) set, that we will also index below. For the Ising model above, for example, the relevant set of local operators is given by {I,S x ,S z }.
The virtual bond state space is given by a start state [α = 1, or equivalently, (1, 0, 0, . . .) T ], a stop state [α = 2, or equivalently, (0, 1, 0, . . .) T ], followed by α = 3, . . . , m int + 2 which assigns an index position to every one of the m int interaction terms in the Hamiltonian that stretches across a given bond in the MPO (strictly speaking, m int corresponds to the number of operators that need to be stored across a given bond which may be less than the number of elementary interaction terms in the Hamiltonian if interaction terms can be grouped by factorizing out specific operators). Hence the dimension of the virtual bond state space is given by D H = 2 + m int . For example, for the Ising model above, the i-th bond in the system in between sites i ≤ i and j ≥ i + 1 carries the single interaction term S x i S x i+1 , hence m int = 1 and D H = 3. The general strategy then for setting up the MPO w.r.t the specific example of the Ising model is as follows: starting from the left end, the bond state space, i.e., the automaton is initialized in the start state (α = 1). This is carried through the MPO (therefore T 1,1 = I) until an interaction term in the Hamiltonian occurs, say at site i, which brings the automaton into the state α = 3 (therefore T 1,3 = S x ). Having only nearest neighbor terms, the subsequent T tensor e.g., at site The values α1, α2, . . . on given A-tensors are reserved for this very specific interaction term, where in general the indices αi > 2 are not all the same. Note that the arrows on the red line are reversed w.r.t. site j which thus need to incorporate the Clebsch-Gordan coefficients (1, 1; 0) that combines the S = 1 multiplet of the spinor S with its dual (also S = 1) into a scalar. It is this CGC on the A-tensor of site j within the super-MPS that takes properly care of the dagger in the scalar productŜi ·Ŝ † j . Similarly, also the coupling strength J is encoded with the A-tensor in the super-MPS.
i + 1 immediately brings down the automaton to the end state α = 2 (therefore T 3,2 = S x ). By having completed the interaction term, the automaton stays in that state (hence T 2,2 = I).
Overall, what has been encoded this way was simple the interaction term With the same line of arguments, the local transverse field term, say at site i, is described by T 1,2 = −hS z , which directly brings the automaton from the start into the end state. By translational invariance, there is nothing special about site i, though. Therefore all of the matrix elements of the tensor T specified above must hold for every site. Given these local tensors in MPO, the summation (trace) over geometric indices α is equivalent to adding up all the interaction terms in the total Hamiltonian.

From super-MPS to MPO
In the presence of global continuous symmetries, all state spaces must be organized into symmetry multiplets. Naturally, this also implies a directedness of lines in a tensor network. From the point of view of a given tensor, the direction on its lines indicate bra or ket nature of these state spaces which, in pictorial language, is equivalent to legs (lines) entering or leaving a given tensor, respectively. For U(1) symmetries it implies that the sum of all charges that enter a tensor must be exactly equal to the sum of all charges leaving it. For SU (2) symmetries, the fusion of all ingoing lines must result in a symmetry sector that exactly matches a symmetry sector resulting from the fusion of all outgoing lines. If all lines are ingoing, the tensor must be scalar in that the (skipped) outgoing line transforms as a singleton index that transforms like the vacuum state (and vice versa, if all lines are outgoing).
In contrast to the Ising model above which has no simple continuous symmetry for h = 0, let us continue with the model system of interest in this work, the (anisotropic) Heisenberg model [compare Eq. (2), using J : where we temporarily introduce hats on top of operators, in order to differentiate them from symmetry labels (e.g.,Ŝ z vs. S z ). The model in Eq. (B1) is U(1) symmetric as it preserves S z tot . In the isotropic case, ∆ = 1, it becomes SU(2) spin symmetric. Then the spin operators need to be grouped into a spinor,Ŝ such that Eq. (B1) can be rewritten in SU(2) invariant form, Note that the relative weights and signs in Eq. (B2) are important for consistency with standard conventions on SU(2) spin multiplets. In particular, the operators in the spinor in Eq. (B2) exactly represent, top to bottom, the states S z = (+1, 0, −1) of an S = 1 spin multiplet (e.g., see Ref. [48]). In contrast, the Hermitian set of operators (Ŝ x ,Ŝ y ,Ŝ z ) does not. However, using Eq. (B2), the dagger on one of the spinors in Eq. (B3) is important.
In general, a local operator acting on some physical site is a spinor, i.e., a collection of operators that transforms like some multiplet q [cf. Fig. A.1(a)]. This can be written as the irreducible operator (IROP)X [nq;qz] where the composite index (nq; q z ) naturally specifies entire state spaces [48], or here an operator space: the index n differentiates between local IROPs that transform according to the same irreducible representation q. By definition of an index, n = 1, 2, . . ., we therefore also introduces an arbitrary but fixed order to the local operators. The label q z , finally, fully differentiates the operators within a given spinor [48]. For example, within SU(2), q z simply stands for S z .
The matrix elements of IROPs are determined via the Wigner-Eckart theorem. For a generic spinor, this IROP acquires a third dimension, which indexes the operators in the irreducible set. Scalar operators then are special. With one in-and one outgoing index, the third index having q = 0 is a trivial singleton dimension that may safely be skipped. In this sense, scalar operators can be reduced to rank-2, and are block-diagonal. For U(1) spin symmetry, for example, scalar operators are the identity operatorÎ or the spin projection op-eratorŜ z . In contrast, the operatorsŜ ± carry q = ± 1 2 , hence switch between symmetry sectors, and therefore are not considered scalar operators.
All local operators eventually can be combined into a single rank-3 tensorX [cf. Figs. A.1(a,b)]. The third index then represents the state space |nq; q z [48] of the "supervectors" X [nq;qz] . For efficiency, the set of local operators should be orthogonal in the sense with arbitrary normalization, otherwise. This is also in the spirit of an orthogonal local (super-) state space of a (super-) MPS. Conversely, assuming that the tensor T of the MPO is given, note that the intermediate supervector index that connects the super-MPS with the local operatorsX [cf. Fig. A.1(c)] may also be generated by the reverse operation of splitting off the local state space (σ, σ ) from the tensor T via SVD. Then by construction, the operators inX would be orthonormal. For both, conceptual and implementational transparency, we can construct an MPO as a super-MPS of operators. By this we mean, that the local state space of the super-MPS are "superstates" that actually refer to a set of orthogonal local operators [e.g., see j † acting on site j that transforms according to exactly the same irreducible representation (typicallyŶ =X; here we also ignore 3-or more-site interactions). This observation holds both, for abelian and non-abelian symmetries.
The construction of the super-MPS that encodes the MPO is greatly simplified by the simple bilinear structure of 2-site interactions as in Eq. (B1) or Eq. (B3). In partiuclar, the super-MPS can be built completely analogous to the automata approach above, while paying simple attention to symmetry sectors. The start (α = 1) and the stop (α = 2) state on the virtual bonds transform like scalars (i.e., have q = 0), whereas the bond states α > 2 directly inherit the symmetry labels from the underlying IROPs in the 2-site interactions [cf. Fig. A.1(e)].
For the Heisenberg model in Eq. (B1), the set of local operators is given byX = {Î,Ŝ z ,Ŝ + ,Ŝ − } for the U(1) symmetric setup, and byX = {Î,Ŝ} for the SU(2) symmetric setup. In either case, the set of local operators is orthogonal as in Eq. (B4). Note also that while in the U(1) symmetric case, also the daggered operator (Ŝ + ) † =Ŝ − appears in the set, this is not the case for the SU(2) symmetric case, since SU(2) is self-dual, and thereforeŜ † i ·Ŝ j =Ŝ i ·Ŝ † j . Specifically, with U [1] ∝ (1, 1; 0) a unitary transformation that corresponds to the Clebsch-Gordan coefficients which combine a spin S = 1 with its dual (again S = 1) into a singlet [cf.
1j-symbol earlier], the spin-spin interaction can be written aŝ Therefore the action of the dagger on one of the spin operators can be transferred via the unitary 1j-symbol into the definition of the super-MPS itself, proper sign-convention on U [1] implied.
For more complicated cases, like the snake MPO representation of 2D Heisenberg Hamiltonian, longer range interactions need to be included. This is straightforward in the automata construction above, yet requires that the bond dimension D H increases [see also Fig. A.1(e)]. The period of the translational invariance of the MPO also increases from 1 to the width W of the system and hence requires at least W different A-tensors in the super-MPS [cf. Fig. A.1(c-d)].
Once the super-MPS is obtained, one can use standard MPS techniques to check whether it can be compressed. An important ingredient here is that the local supervector space is orthogonal, indeed [cf. Eq. (B4)]. If the bond-dimension D H can be reduced at no cost, i.e., by discarding singular values that are strictly zero, the super-MPS and subsequently the MPO contains inefficiencies that may be simply removed with an improved setup of the super-MPS itself. On the other hand, for long-ranged systems the bond dimension D H may simply become too large, in practice, for an exact representation of the Hamiltonian. In this case, standard MPS truncation techniques may be employed on the level of the super-MPS itself. Here a uniform normalization of the supervector space (i.e., the local operators) is advised such that standard MPS techniques are directly applicable without any further ado. Alternatively, one may truncate on the level of the MPO, either by SVD or variational techniques. The latter is unavoidable for MPO products or sums in any case, as will be discussed next.

Appendix C: Initialization in the XTRG algorithm
In this section, we compare three different initializations of ρ(τ 0 ) in the XTRG algorithm. The quality of our initial ρ(τ 0 ) for small τ 0 is measured by estimating the relative error of the free energy |δF (β)/F (β)|, starting from exponentially small β = 2τ 0 (i.e., the first data point after initialization at β = τ 0 ) down to intermediate temperatures β 10. Interestingly, first-order Trotter-Suzuki initialization manages to arrive at an MPO with bond dimension that is lower than what is required for the actual representation of the Hamiltonian itself. But as a consequence, the overall errors are also larger. For comparison, nevertheless, we also show data initialized via SETTN at the same D * 0 = 2 (D 0 = 4; green data). The resulting errors are much lower 10 −10 and 10 −15 for τ 0 = 0.1 and 0.01, respectively.
By increasing the initial bond dimension to D * H and slightly above, i.e., D * 0 = 3, 4 (D 0 = 5, 8), as seen in Fig. A.2, SETTN initialization as in Eq. (7a) can offer generally better accuracy. The initial inaccuarcy represents a systematic error that also propagates along the XTRG procedure towards lower temperatures. e.g., for β ∼ 10, the SETTN initialization with D * 0 = 2 is still an order of magnitude more accurate as compared to the Trotter-initialized data, and orders of magnitude more accurate if D 0 is only marginally increased to D * 0 = 3, 4 relative to the D * = 150 of the subsequent XTRG.

Linear initialization ρ(τ0) I − τ0H
By sweeping over several orders of energy scales, the XTRG algorithm permits a simple linear initialization ρ(τ 0 ) I − τ 0 H at basically infinitesimal, i.e., exponentially small τ 0 . In the following, we analyze the effect of the initial τ 0 in more detail.
In Fig. A.3, we compare the initialization of ρ(τ 0 ) with SETTN [panel (a)] and linear expansion [panel (b)], respectively. We benchmark the accuracy by analyzing the relative errors of the free energy δF/F , with various initial τ 0 = 10 −8 , . . . , 10 −2 . In Fig. A.3(a) the initial state ρ(β = τ 0 ) is constructed via SETTN [cf. Eq. (7a)] and constrained to the bond dimension D * 0 for β = τ 0 only. For all other points at lower temperatures β n ≡ 2 n τ 0 with n = 1, 2, . . . we use D * = 200. With this setup, the data in Fig. A.3(a) at the lowest temperatures (largest β) exhibits a similar level of accuracy, irrespective of the choice of the initial τ 0 10 −4 for D * 0 ≥ D * H = 6. Note that D * 0 = D * H is the minimal bond dimension to represent the lowest order linear expansion ρ(τ 0 ) I − τ 0 H. The long extended straight slopes for τ 0 = 10 −6 or τ 0 = 10 −8 are simply ∝ β, as indicated by the guide to the eye (purple dashed line), which just indicates that the accuracy is limited by accumulated double precision error of the calculation. The strong upturn for β 0.1 then is where truncation error sets in. Using larger initial τ 0 , e.g., τ 0 = 10 −4 (yellow curve) error accumulates more strongly, which implies that D * 0 = 6 already starts to affect the accuracy at larger β. Clearly, for τ 10 −4 a larger D * 0 , i.e., higherorder terms in the series expansion are required to maintain accuracy. For example, the data for τ 0 = 10 −2 and D * 0 = 30 again shows good accuracy in comparison.
Once τ 0 is small enough, plain simple lowest-order linear expansion suffices. This is analyzed in Fig. A.3(b) where we replace series expansion for the initialization of ρ(τ 0 ) by the MPO for I − τ 0 H. This is an extremely convenient initialization that can be simply derived from the MPO for the Hamiltonian. In particular, the MPO for this initial ρ(τ 0 ) can be exactly represented with bond dimension D * 0 = D * H . Starting with β = τ 0 ≪ 10 −2 , XTRG can be used to exponentially decrease temperature down to values β ≈ 10 −2 which still may be considered part of the initialization of ρ at large temperatures before actual physical energy scales set in. Therefore while we always have D * 0 = D * H for the very first step, by definition, we can alo constrain D * 0 to the values specified for the entire range β ≤ 10 −2 [vertical dashed line in Fig. A.3(b)]. For β > 10 −2 we allow the MPO to grow up to dimension D * = 200 to capture the quickly growing MPO entanglement.
Using D * 0 = D * H up to β = 10 −2 [squares in Fig. A.3(b)], the data already significantly deteriorates at the largest β 1 by about two orders of magnitude. This shows that D * 0 = D * H up to β = 10 −2 is simply too small, as it introduces systematic errors. By increasing D * 0 modestly, i.e., D * 0 = 30 which is still about an order of magnitude smaller than what is required for large β, the systematic errors reduce dramatically. Starting with τ 0 10 −4 , the accuracy at large β competes with a careful higher-order expansion for finite τ 0 [for reference, we replotted the data from the D * 0 = 30 series-expansion initialization starting with τ 0 = 10 −2 from panel (a)]. Similar to panel (a), for our data sets with smallest initial τ 0 , we see wide ranges where the numerical error is simply ∝ β, and hence given by accumulated double precision error and, in particular, not truncation error.
The initialization procedure in Fig. A.3(b) above started from infinitesimally small β = τ and worked its way up exponentially using XTRG to β = 10 −2 . Using D * H < D * 0 D * (final) , the overall numerical lost of the entire calculation is strongly dominated by the simulation of β 1. For example, for D * 0 = 6 (30) and τ 0 = 10 −7 , the cost for β ≤ 10 −2 is about 7% and 11% of the total calculation, respectively. In this sense the procedure above is an extremely simple, efficient, and accurate initialization for finite temperature calculations. The algorithm works for arbitrary topologies of Hamiltonians, including long-range interactions or higher-dimensional systems.
Nevertheless, since the series expansion scheme serves as a systematic way to provide accurate initialization at finite τ 0 < 0.1, in this work, we sticked to initialize of the density operator ρ(τ 0 ) with our already existing codes on series expansion. Appendix D: Compression of matrix product operators MPO compression is of key importance and frequently used in the XTRG and SETTN algorithm, to compress the product or sum of two MPOs. The overall procedure follows standard MPS strategies, where the implications of abelian or nonabelian symmetries can be largely put aside as an extremely convenient benefit of using the QSpace tensor library.
Overall, the variational method is preferable due to its higher efficiency, while the direct SVD compression is also useful as long as the bond-dimension D is manageable. In the following, we focus on the variational compression, but also provide details of the SVD compression along the way.

Compression of MPO product
Consider the product of two MPSs A * B = C, with bond dimensions D a , D b , and D c , respectively. The representation of the product is exact if D c = D a D b which, how-ever, is typically numerically costly. Therefore C needs to be compressed in an efficient and numerically controlled manner. A typical example in this paper is the representation of C := H n X = H * (H n−1 X), with X = I or ρ(β) in SETTN, where we need to project A := H to a previously e.g., iteratively obtained MPO for B := H n−1 X and then compress the "fat" MPO to bring down its bond dimensions. Similarly, in XTRG, one needs to apply ρ(β) onto itself, i.e., A = B := ρ(β), in order to reach the density matrix C := ρ(2β) at half the temperature.
The direct SVD compression is straightforward [42] but computationally costly. Given the product of two MPOs, as depicted pictorially in Fig. 2(b), horizontal lines are fused pairwise with respect to the same horizontal bond position into a single "fat" index of dimension D a D b , and then truncated via SVD. Using symmetries, abelian or non-abelian symmetries alike, the fusion step includes a simple tensor product of two state spaces. For this it is important, that the arrows along the virtual bond state spaces (horizontal lines) are parallel and point all in the same direction [ Fig. 2(b)].
The computational cost of SVD compression scales as O D 3 a D 3 b . For the series expansions, say, when constructing H n ρ with D a = D H and D b = D, this is still relatively cheap, O D 3 , given that typically D H D. In contrast, the numerical cost becomes quickly prohibitive for XTRG since with D a = D b = D the cost scales as O D 6 .
The variational method can significantly reduce numerical cost, and therefore is mainly adopted in the present study. For this, we use a two-site update similar to standard DMRG procedures to allow adaptive adjustment of bond dimensions. This is particularly important when exploiting symmetries, abelian and non-abelian alike, to optimally adapt bond dimensions w.r.t. to each individual symmetry sector. This way, irrelevant symmetry sectors drop out automatically, whereas possibly new relevant symmetry sectors can emerge or get strengthened.
For the variational approach, we minimize the cost function (Frobenius norm squared), We then take the partial derivative with respect to the product of two adjacent local tensors C i C i+1 in the full MPO of C. This results in the linear system of equations, to be solved iteratively for i = 1, . . . , L − 1, with L the length of the MPO (while all simulations in this work are based on real numbers, for the simplicity of the derivation, nevertheless, we assume complex numbers). Both sides of Eq. (D2) can be expressed as fully contracted tensor networks, except for the missing tensors C * i and C * i+1 (i.e., with "punched holes"), as shown in Fig. A.4(a). Therefore both sides of Eq. (D2) represent tensors of rank-6. For a canonicalized MPO C, where all lines in C are directed towards the orthogonality center, here at sites (i, i + 1), the left-hand side is simply C i C i+1 . The right-hand side defines the generalized overlap tensor C E [cf. Fig. A.4(a)]. Equation (D2) therefore directly states the solution for (C i C i+1 ) to the local optimization problem w.r.t. to sites (i, i + 1).
To obtain C E , one needs to iteratively update the left/right environment tensors V L /V R of the three-MPO product until the structure in the right-hand-side of Fig. A.4(a) reached. As shown in Fig. A.4(c), to update the V L /V R tensors, we iteratively contract the local tensors A i , B i and C i of MPO A, B and C with it. The procedures of updating V L /V R constitute the most time-consuming ones, which scale as O D 4 (assuming D a = D b = D c =: D), when we square ρ and compress it in XTRG. Nevertheless, this is still computationally much cheaper as compared to the SVD compression of O D 6 cost as briefly discussed above.
Given the optimized product C i C i+1 = C E , we need to split C E into the actual product shape C i C i+1 by performing SVD, C E = U Λ V † as shown schematically in Fig. A.4 (d).
To gauge the MPO C in a canonical form, C i is updated with U , and C i+1 with ΛV † in a left-to-right sweep, while in a right-to-left sweep, the matrix Λ is contracted with U and thus associated with C i , instead. One typically only needs a few full sweeps (left to right and vice versa) to converge the cost function and obtain the optimal MPO for C. In this work, typically at most 4 sweeps were sufficient.

Compression of MPO sum
Summation of MPOs is an essential technique, e.g., in the series expansion ρ(β) = N n=0 n! H n . Here we generalize standard procedures for the addition of MPS to MPO. For this, we also resort to a 2-site variational approach as illustrated in Fig. A.4(b). To find a optimal MPO for C = A + B, where the generalization to more than two vectors is straightforward, we minimize the cost function, Again we take the partial derivative with respect to the product of two adjacent local tensors C i C i+1 of the MPO C, resulting in the linear system of equations, to be solved iteratively for i = 1, . . . , L − 1. Again, using a canonicalized MPO for C, Eq. (D4) simply reduces to C i C i+1 = C A + C B , with C A etc. generalized overlap matrices to be computed iteratively [cf. Fig. A.4(b)]. The remainder of the algorithm proceeds exactly the same as the compression of MPO products above.
Appendix E: Series expansion tensor network simulations with linear versus exponential β grid In this appendix, we compare SETTN calculations with three schemes of selecting expansion point set: (M) Maclaurin scheme which expands ρ(β) around β = 0, i.e., there is only one expansion point in the set; (L) point-wise Taylor expansion around a linear β set, β = nτ 0 with n an integer; (X) exponential set β = 2 n τ 0 . The results are summarized in Fig. A.5.
In Fig. A.5(a), we compare the accuracy of the above grid setups (M,L,X) in the calculation of the free energy for an L = 16 Heisenberg chain. Compared to (M), both (L) and (X) are clearly superior, as they gain four orders of magnitude in accuracy for the largest β.
We also compare these three grid setups within SETTN to the more challenging system of a 4×4 SLH. Here due to the significantly larger truncation errors across all approaches, the gain of (L,X) over (M) is significantly less pronounced, as seen in Fig. A.5(b). To reduce the truncation error, the number of virtual bond states would have to be increased significantly from the D * = 200 or even 400 in the present calculation. As for the numerical efficiency, although having very good accuracy, the computational overhead of the linear scheme (L) is significant in any case. In contrast, the logarithmic scheme (X) strongly reduces the computational cost [see caption to  For completeness, we provide analytical expression of the partition function in a 1D XY model. The Hamiltonian is given by Exploiting Jordan-Wigner transformation, the Hamiltonian can be mapped onto a plain fermionic tight-binding chain, which for open boundary condition (OBC), is diagonalized by the 1-particle eigenstates, where c k ] are fermionic annihilation (creation) operators at site i ("momentum" k), respectively [9,11].
The partition function is fully determined by the dispersion k , with free energy F = −T ln Z.
Appendix G: Entanglement measurements of XXZ model on the cylinder geometry The entanglement landscape, along with its top and side views, for the 2D XXZ model on a square-lattice cylinder (L = 10, W = 5, ∆ = 5) is analyzed in Fig. A.6. Similar to the analysis of the larger system in Fig. 13, we can again observe prominent peaks which can be help to pinpoint the critical temperature T c [vertical dashed lines in in the very center of the system deviates slightly from T c = 0.56 in that it approaches from above, i.e. T S c > T c . This is different from Fig. 13, in the main text where T S c approaches T c from the low-temperature side. The underlying reason is the different boundary conditions, i.e., cylindrical BC here vs. fully open BC in Fig. 13.
Interestingly, we find that T S c determined from cylinder geometry again serves as a good estimate for T c . As shown in Fig. A.6(c), we can see that the profile of entanglement curves overall (i.e., selecting the maximal entanglement over bonds for any given temperature), dubbed entanglement envelope, shows a peak in close proximity to T c = 0.56.