Energy Scales and Exponential Speedup in Thermal Tensor Network Simulations

# Energy Scales and Exponential Speedup in Thermal Tensor Network Simulations

Bin-Bin Chen Department of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and Physics (Ministry of Education), Beihang University, Beijing 100191, China    Lei Chen Department of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and Physics (Ministry of Education), Beihang University, Beijing 100191, China    Ziyu Chen Department of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and Physics (Ministry of Education), Beihang University, Beijing 100191, China    Wei Li Department of Physics, Key Laboratory of Micro-Nano Measurement-Manipulation and Physics (Ministry of Education), Beihang University, Beijing 100191, China International Research Institute of Multidisciplinary Science, Beihang University, Beijing 100191, China    Andreas Weichselbaum Physics Department, Arnold Sommerfeld Center for Theoretical Physics, and Center for NanoScience, Ludwig-Maximilians-Universität, Theresienstrasse 37, 80333 Munich, Germany
July 4, 2019
###### Abstract

In this work, we exponentially speed up thermal simulations of quantum many-body systems in both one- (1D) and two-dimensional (2D) models using matrix product operators (MPOs). Instead of evolving the density operator linearly in inverse temperature as conventional Trotter-Suzuki methods do, we cool down the system by iteratively projecting the MPO representation of to itself, i.e., doubling in the process of imaginary time evolution. This exponential tensor renormalization group (XTRG) scheme, compared to linear evolution schemes, reaches low temperatures much faster, and thus not only saves computational time but also merits better accuracy due to significantly fewer truncation steps. For similar reasons, we also find that the series expansion thermal tensor network (SETTN) approach benefits in both efficiency and precision, from the logarithmic temperature scale setup. For both thermal algorithms, XTRG as well as SETTN, we fully implement non-Abelian and Abelian symmetries to greatly enhance their numerical performance. We employ these cutting-edge techniques for finite temperature simulations to explore low-temperature thermal states of both 1D and 2D Heisenberg models. The entanglement properties, as well as the renormalization group flow of entanglement spectra in MPOs, are discussed, where logarithmic entropies () are shown in both spin chains and square lattice models with gapless tower of states.

###### pacs:
05.10.Cc, 05.30.-d, 75.10.Jm

## 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) Xiang et al. (2017); Lake et al. (2013) and spin liquid materials in two dimensions (2D) Balents (2010); Nasu et al. (2015); Shen et al. (2016); Do et al. (2017); Yamashita et al. (2017); Kelly et al. (2016) 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. For example, thermal states near 1D quantum critical point are revealed to show universal entropy in the partition function due to emergent conformal symmetry Tu (2017); Tang et al. (2017); Chen et al. (2017a) in the low-energy regime.

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 with 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 Eisert et al. (2010) (as well as operator space entanglement Žnidarič et al. (2008)). However, it was estimated that the required MPO bond dimension has an upper bound scaling as Hastings (2006) which still seems to pose a severe barrier towards obtaining low- properties.

Nevertheless, on the other hand, various renormalization group algorithms have been proposed to accurately compute thermodynamics in the (1+1)D problems, practically even down to very low temperatures. These methods include the transfer matrix renormalization group (TMRG) Bursill et al. (1996); Wang and Xiang (1997); Xiang (1998) and finite-temperature DMRG Feiguin and White (2005) which are based either on traditional density matrix renormalization group (DMRG), or tensor network algorithms such as the linearized tensor renormalization group (LTRG) Li et al. (2011); Ran et al. (2012); Dong et al. (2017) and variational projected entangled pair operator (PEPO) method in (2+1)D Czarnik et al. (2012); Czarnik and Dziarmaga (2015); Czarnik et al. (2016). Besides, a combination of finite-temperature DMRG and Monte Carlo samplings called minimally entangled typical thermal states (METTS) was proposed White (2009); Stoudenmire and White (2010) and recently generalized to (2+1)D Bruognolo et al. (2017). The success of these algorithm in (1+1)D, and partly in (2+1)D, strongly suggests that does not necessarily scale exponentially with .

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. Prosen and Pižorn, 2007; Žnidarič et al., 2008. It has been revealed that this MPO entanglement saturates for gapped systems and scales logarithmically (as ) for quantum critical chains Prosen and Pižorn (2007); Žnidarič et al. (2008), which ideally characterizes the computational complexity for simulating a thermal CFT in classical computers. Very recently, Barthel Barthel (2017) deploys the time-dependent DMRG methods (adapted to imaginary time), as well as CFT arguments, to show that the related Renyi entropy generally scales as .

Intuitively, this scaling can be understood simply by considering finite-size spectra with many-body low-energy level spacing . In order to sample a thermal average, it must hold that at the very least , or equivalently (in other words, finite introduces an effective cutoff to the system size). Now given that low-energy states are beyond area law due to logarithmic corrections, one has the Renyi entropy Bazavov et al. (2017) for individual low-energy states using OBC. By choosing with fixed constant , and by going from individual low-energy pure states to a thermal state with weights , i.e., , the overall entropy acquires another factor of leading to the earlier result above. Consequently, the numerical cost of simulating a thermal state on a sufficiently large system with OBC scales similar vs. as a ground state calculation vs. using PBC.

The above intuitive argument fits the holographic picture in terms of thermal multiscale entanglement renormalization ansatz (MERA) Evenbly and Vidal (2015), 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 Swingle (2012), i.e., is determined by an effective system size at low temperatures.

Furthermore, the argument of translating the scaling of the entropy in to that of is completely consistent with the notion within CFT that and are equivalent directions connected via a modular transformation. This has direct consequences for conformal TTN framework in (1+1) dimensions, i.e., with one spatial (horizontal) and one imaginary time (vertical) axis. The horizontal transfer matrix across different temperatures has the ground state of as its dominant eigenvector which thus contains logarithmic entanglement. For thermodynamics of an infinite-size quantum chain (), we therefore expect that the vertical transfer matrix (across different real-space sites) also has a dominant eigenvector containing logarithmic entropy, i.e., . 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 Barthel (2017). This suggests that the bond dimension 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. 1 by doubling , one can therefore design an exponentially faster cooling procedure, in contrast to current standard simulation techniques which linearly evolve the full density matrix Bursill et al. (1996); Wang and Xiang (1997); Xiang (1998); Feiguin and White (2005); Li et al. (2011); Ran et al. (2012); Dong et al. (2017) or the typical sampling states White (2009); Stoudenmire and White (2010); Bruognolo et al. (2017) in imaginary time. We note that essentially a similar, even though much more involved strategy was pursued in Refs. Czarnik and Dziarmaga, 2015; Czarnik et al., 2016. Their approach was based on a dimensional reduction via a nested contraction of linear Trotter gates, followed by a variational optimization of a tree tensor network. In contrast, our approach does not rely on Trotter gates, and hence is straightforwardly applicable to arbitrary Hamiltonians (1D and 2D). Overall, it represents 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) Wilson (1975); Bulla et al. (2008); Weichselbaum (2012a). There also, with every new iteration the energy scale is reduced by a factor, typically . 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 SU(2) symmetric spin-half Heisenberg models

 H=J∑⟨i,j⟩→Si⋅→Sj (1)

both, in 1D Heisenberg chains of length , as well as in the 2D square lattice for systems of width and length , and thus with a total of sites, using open (OBC) as well as periodic boundary conditions (PBC). We only include nearest neighbor couplings as indicated by the sum . For the purpose of benchmarking, we also consider the U(1) symmetric spin-half XY-chain with , i.e., for OBC,

 H=J∑⟨i,j⟩(SxiSxj+SyiSyj) (2)

as this can be mapped to a fermionic tight-binding chain. Symmetries, whether non-abelian or abelian, are fully exploited, throughout. We also set as the unit of energy, unless specified otherwise. Furthermore we use units .

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 Chen et al. (2017b), named pointwise Taylor expansion algorithm, exploiting 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.

## Ii Symmetric Thermal Tensor Networks in Logarithmic Temperature Scale

### ii.1 Symmetric Matrix Product Operator

Symmetries appear very naturally in TTNs. By construction, the thermal state itself inherits all symmetries of the Hamiltonian , which then also needs to be preserved on the algorithmic level. For example, symmetries are also preserved for Trotter-Suzuki type TTNs Li et al. (2011); Ran et al. (2012), where the local tensors (storing Boltzmann weights) are symmetric. Similarly, in the series-expansion TTNs, it is also clear that arbitrary ’s have exactly the same symmetry as (any unitary symmetry transformation that leaves intact also leaves intact), and so does the resulting tensor network representation of .

The explicit implementation of non-Abelian symmetries has been regarded as a standard technique in ground state DMRG simulations () McCulloch and Gulácsi (2002), which has many important applications including exploring quantum spin liquids in frustrated quantum magnets Depenbrock et al. (2012); Gong et al. (2014), and is also shown to be useful in METTS type thermal simulations Bruognolo et al. (2015); Binder and Barthel (2017). 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 Weichselbaum (2012b), 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 . As this involves reduced matrix elements in the spirit of Wigner-Eckart theorem Weichselbaum (2012b), we refer to this as the reduced MPO in contrast to the full MPO in the absence of or only abelian symmetry representation. By switching from a state-based to a multiplet-based description, we can reduce the overall bond dimension from states to multiplets. The reduced MPO representation of can be constructed by automata Crosswhite and Bacon (2008); Fröwis et al. (2010); Pirvu et al. (2010) method for 1D Hamiltonians and MPO sum-and-compress scheme Hubig et al. (2017) for more complicated 2D lattice models. For a Heisenberg chain with nearest-neighbor interactions, a full MPO requires bond states. As this corresponds to two singlets and one triplet, i.e., where specifies multiplets of dimension each, the reduced SU(2) invariant MPO only involves multiplets. For the Heisenberg model on a 2D square lattice, we map a system of width to a 1D snake-like chain with “long-range” interactions (up to distance ). Then e.g., using OBC, the full MPO requires bond states, while the reduced SU(2) invariant MPO has a significantly more compact representation with only multiplets (). More details on the symmetric MPO representation of total Hamiltonian can be found in App. A.

The computational cost in a tensor network algorithm typically scales with some power aside other factors concerning number of sites etc., where for the MPO structure of this paper we encounter . By exploiting non-abelian symmetries, the computational cost can be effectively reduced to which leads to a gain in numerical efficiency by . For a single SU(2) symmetry, it roughly holds, on average, 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 since after all, with the overall scale factor absorbed into the definition of .

### ii.2 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 with some constant 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 system at finite temperature, with high efficiency and accuracy. We start by preparing an MPO of the (unnormalized) thermal state 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

 ρ0≡ρ(τ)→ρ1≡ρ(τ)⋅ρ(τ)=ρ(2τ). (3)

Feeding the last MPO iteratively into the next step, with and therefore , we obtain,

 ρn≡ρ(τn)=ρn−1⋅ρn−1 . (4)

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 with the largest local energy scale here given by the Heisenberg coupling strength, we can use an efficient series expansion scheme [cf. Sec. II.3]. For as small as even simplest lowest-order linear expansion of can suffice, which extremely simplifies initialization even for longer-ranged Hamiltonians e.g., that are not suitable for Trotter like decomposition 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 is provided in App. B.

Given an MPO representation for , we can compute the (unnormalized) thermal state at temperature via , i.e., by contracting with its conjugate. This guarantees positivity of the thermal state even in the presence of truncation of the MPO for . Furthermore, we can also compute the partition function at via , and thus gain another factor of two to reach lower temperatures [see Fig. 1(b)]. The latter can be simply obtained by computing the Frobenius norm squared of . Not incidentally, many of the features above are directly related with common procedures within the setup of a purified thermal state. Bursill et al. (1996); Wang and Xiang (1997); Li et al. (2011); Xiang (1998); Feiguin and White (2005); Dong et al. (2017); Czarnik and Dziarmaga (2015); Czarnik et al. (2016); Schwarz et al. (2017).

In case the grid of inverse temperature values is too sparse, intermediate values can be easily obtained by shifting the initial value of

 τ0→τ0⋅2zwith z∈[0,1) , (5a) a procedure that is entirely analogous to z-shifts within the NRG. In order to obtain a uniform logarithmic grid over nz shifts, one may simply choose zi=inzwith i=0,…,nz−1% . (5b)

Different “-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 for various .

### ii.3 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 Chen et al. (2017b). By exploiting the series expansion of density matrix in Eq. (6a), SETTN is essentially free of discretization error, making it distinct from previous Trotter-Suzuki type RG methods including TMRG Wang and Xiang (1997); Xiang (1998), finite-temperature DMRG Feiguin and White (2005), LTRG Li et al. (2011); Ran et al. (2012); Dong et al. (2017), and METTS White (2009), etc. The efficient MPO representations of 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 long-range 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 for small , a series expansion yields [cf. Fig. 2(a)]

 ρ0≡ρ(τ0)≃Nc∑k=0(−τ0)kk!Hk . (6a) The required cutoff order of the expansion is Nc∼Nτ0 , i.e., proportional to the total number of sites N. In practice, Nc 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 Nc≲10. The high-temperature Maclaurin expansion in Eq. (6a) can be employed not only in the intialization stage, but also for simulating low-temperature thermal states, as shown in Ref. Chen et al., 2017b. Despite its competitive performance, this method still leaves room for further improvement. Since Eq. (6a) expands ρ around the infinitely high temperature, i.e., β=0, the power series in Hn involves large Nc∝Nβ for large system size N and low temperature 1/β. The precision of SETTN is limited by the truncation in Hn [see Sec. II.4], which generally increases as n increases Chen et al. (2017b). In this sense, a pointwise Taylor expansion can help reduce the expansion order Nc and improve the accuracy, i.e., ρ(τn+1)≃(Nc∑k=0(τn−τn+1)kk!Hk) e−τnH≡ρn . (6b)

Equation (6b) expands the density operator around an arbitrary but fixed . For generality, the initialization in Eq. (6a) may be viewed as iteration having and .

Now given the density operator obtained by initialization via Eq. (6a), with can be obtained via Taylor expansion around . In particular, this also hold for 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 . Given , the MPO representation of can be expanded as in Eq. (6b). Each term in the summation there can be obtained iteratively by projecting onto and compressing the product. For the overall sum then [cf. Fig. 2(b)] we also employ variational optimization (see App. C) to finally arrive at the MPO representation of .

By repeating this procedure, we also can follow the XTRG protocol to cool down the temperature along the inverse temperature grid . In contrast to the plain doubling scheme in XTRG, however, in case of SETTN the step size 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 with only modest cost.

Note that for this improved SETTN, as we will refer to it, using an exponentially increasing series does not acquire exponential acceleration as XTRG does, since in the case of SETTN one still needs to perform projection and compression operations times. Nevertheless, from the point of view of SETTN, the exponential series is computationally preferable to, say, a linear series as expansion points, since the former can reduce expansion overhead and thus save computational time, in practice, without losing any accuracy (see App. D for a detailed comparison). Overall, SETTN can serve as useful complementary approach to XTRG methods.

### ii.4 MPO compression and numerical cost

In SETTN we start with a reduced SU(2) invariant MPO for . Then we iteratively apply the projections onto to obtain , with Eq. (6a) represented by . 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 can affect the accuracy for large . Hence we truncate by number of multiplets, throughout. For this, we introduce the control parameters which stand for the maximum number of multiplets to be kept in the -th iterative term when computing . For simplicity, we set this parameter constant, i.e., , which also stands for the bond dimension of the target state . Furthermore, we choose constant but, for the sake of the analysis, may use a different value for for the initialization in Eq. (6a) if specified.

For an extremely small (say, as small as to ), the initialization of can be simplified to lowest-order, i.e., linear expansion . Having in Eq. (6a), the result shares the same bond-dimension as itself. In constrast, when expanding around finite as in Eq. (6b), the bond dimension in typically needs to grow significantly, and therefore is fixed to some specified .

The compression of the SETTN projections above can be achieved either by a singular value decomposition (SVD) technique quite similar to that in Ref. Chen et al. (2017b), apart from the fact that the MPOs here have SU(2) symmetry, or by a variational optimization which can greatly improve numerical efficiency (see App. C for more details on related MPO compression techniques). Within SETTN, the cost of either compression scheme scales like . We tested both and found comparable numerical accuracy. Finally, we variationally add up the MPOs for with coefficients as in Eqs. (6) to obtain (cf. App. C.2).

In contrast, XTRG projects onto itself in Eq. (4). So both MPOs involved have large bond dimension . Here a direct SVD compression is numerically costly, , to be compared to for SETTN. For the XTRG iteration in Eq. (4) we therefore constrain ourselves to a variational compression which scales like [see App. C.1].

More explicitly, the numerical cost of SETTN for an entire run up to inverse temperature scales as assuming for large with in Eqs. (6), whereas XTRG scales as . For practical simulations as in Figs. 3 and 4, 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 ( for CFTs, say, spin-1/2 Heisenberg chain, see Ref. Barthel, 2017), we can thus estimate the relative run time of SETTN over XTRG as . Therefore for large. Similarly, also Trotter-Suzuki type linear thermal RG methods, like the finite-temperature DMRG Feiguin and White (2005) and LTRG Li et al. (2011); Dong et al. (2017), with scaling , are much slower by a factor 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 Bruognolo et al. (2017). The numerical (time) cost of the latter scheme scales like , where the additional factor stems from the number of required swap gates which is proportional to the width . For 2D, however, typically such that the relative cost of XTRG scales like . For a typical 2D simulation parameter setting this is of [e.g., with , and (correspondingly ) in SU(2) simulations, or in U(1) calculations one obtains ]. Nevertheless, XTRG is still expected to be clearly advantageous over Trotter gates due to the far fewer truncation steps involved. Besides, XTRG can be effectively parallelized based on -shifts [cf. Eq. (5)].

## Iii Benchmark calculations: 1D and 2D Heisenberg models at finite temperature

An equilibrium thermal state is described by the partition function . Typical interesting thermodynamic quantities, which constitute important tasks for the TTN algorithms to compute, including

 f≡FN=−1NβlnZ(β) free energy (7a) u≡EN=∂(βf)∂β=1NTr[H⋅ρ(β)]Z(β) (internal) energy (7b) cV=∂u∂T=−β2∂u∂β specific heat (7c)

etc. The computation of free energy and energy density , are straightforward, where 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 , i.e.,

 u = 1β∂(βf)∂lnβ (8a) cV = −β∂u∂lnβ , (8b)

instead. This is also more stable numerically for small temperature since, the quotient of numerical differences is divided by for the specific heat in Eq. (8b) and not by as in Eq. (7c), and is multiplied by for the internal energy in Eq. (8a). This formula is used to compute the specific heat in Figs. 5, 6. In order to reduce numerical differential errors, independent calculations with slightly different initial values are run in parallel, e.g., using in Eq. (5b) which produces interleaved data points with , i.e., .

Finally, we note that the LTRG approach adopted in this work, e.g., for the data in Figs. 3 and 4 below, has been streamlined with the remainder of the TTN procedures used in this work. It differs from the original LTRG algorithm in Refs. Li et al. (2011); Dong et al. (2017) in that it successively projects the MPO for 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.

### iii.1 Heisenberg chain

Firstly, we benchmark XTRG results with conventional linear evolution in Fig. 3, where an -site spin-1/2 Heisenberg chain (PBC) [cf. Eq. (1)] is calculated up to . From Fig. 3(a) it is clear that the accuracy in the low- regime gets continuously improved as increases in XTRG. Starting from a fixed , XTRG reaches a precision as good as also at the lowest temperatures for as compared to 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 (), but are clearly less accurate in the low- region where truncation errors dominate. These remarkable results suggest that as XTRG targets the low- 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. 3(b), which offers a quantitative estimate of computational complexity. The entanglement data suggest that truncation errors start to develop for and the simulation errors stop increasing due to the convergence of entanglement for . This is also clearly reflected in the overall error in physical quantities such as the free energy in Fig. 3(a).

Besides the Heisenberg chain, we also benchmark XTRG, LTRG and SETTN for an XY chain [cf. Eq. (2)] with size , where analytical solutions are available (App. E). As shown in Fig. 4, 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 . 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 , utilizing the XTRG algorithm with bond dimension up to . As shown in Fig. 5, in the low- region, the specific heat of the system shows a universal linear relation versus temperature, as indicated by the polynomial fitting (purple dashed line) with . In addition, the fitted slope is also in perfect agreement with the well-known value from CFT prediction Affleck (1986), from which we extract the central charge . On the other hand, in the high- regime, the specific heat is also universal, and decays as (a polynomial fit in the log-log scale, depicted by the yellow dashed line yields an exponent ). This exponent can be confirmed by a high-temperature expansion up to the second order, which approximates the energy as

 E=Tr(e−τHH)Z(τ)≃1Z0[TrH−Tr(H2)τ+(TrH)2Z0τ+O(τ2)]

where , with the high- limit .

### iii.2 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 App. A. 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 () in Fig. 6(a-c) and QMC for larger systems () in Fig. 6(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. 6(a-c), we show the free energy, energy and specific heat results of a 44 Heisenberg square lattice. Very nice agreement between XTRG and ED data is observed in all three plots. As seen in the inset of Fig. 6(a), the relative accuracy is quite high, i.e., () for at low temperatures (). The energy density shown in Fig. 6(b) was obtained by taking derivatives of interleaved XTRG free energy data [cf. Eq. (8a)]. The error in the energy density is small even down to , as seen in the zoom in the low- region in the inset of Fig. 6(b). XTRG data only differ from the ED results in the fourth digit for , and is only bounded by numerical differentiation error for . In Fig. 6(c), we show our results of the specific heat, which was calculated by taking derivatives of energy data as in Eq. (8b). Inset plots on a log-log scale, from which we again observe an algebraic behavior () at high temperatures. In the low- region, it shows a very rapid (exponential) decay versus the temperature.

For a 165 Heisenberg square lattice which is far beyond the scope of ED calculations, we compare our XTRG results to those of quantum Wang-Landau (QWL) simulations Bauer et al. (2011) in Figs. 6(d). We run the calculation down to . For the smallest temperature for which with have well-converged QWL reference data at comparable numerical cost, the error in the data for the free energy is . Since the truncation errors is generally larger on lattice, we extrapolate the free energy in Fig. 6(d) to and observe a perfect agreementwith QWL data, with the error further reduced by about an order of magnitude. Besides the free energy, in Figs. 6(e-f) we also show that the energy density as well as specific heat all have very good accuracy. In the high- region, again shows a relation, as shown in the inset of Fig. 6(f) and as already discussed with Fig. 5.

In Figs. 6(b, e), for comparison we also included METTS data exploiting symmetry only Bruognolo et al. (2017). 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 plot in Fig. 6(e), note that is currently the typical lowest temperatures that 2D METTS simulations can reach Bruognolo et al. (2017) at comparable computational resource as 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.

## Iv Entanglement in Thermal Tensor Networks

### iv.1 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” Zwolak and Vidal (2004)

 |Ψ(β)⟩≡1√Z(β)|e−β2H⟩≡|~Ψ(β)⟩ (9)

which vectorizes the MPO for . 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 , whereas . Note that this definition, is a specific (and most natural) choice of purification Verstraete and Cirac (2004); Feiguin and White (2005), which in some other context is also called the thermofield double (TFD) state , where is the eigen-energy of eigenstate and its duplicate in the auxiliary state space Takahashi and Umezawa (1996); Maldacena (2003); Kallin et al. (2011); Schwarz et al. (2017). 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. (9). Specifically, the entanglement spectrum is given by the eigenspectrum of where is the ‘super’-density-matrix of the purified thermal state as in Eq. (9).

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 Terhal et al. (2002). 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 have attracted recent interest Nguyen et al. (2017); Takayanagi and Umemoto (2017); Barthel (2017).

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 . 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. 7, 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. 7: (i) a low entanglement region , (ii) an intermediate region where entanglement rises quickly, and (iii) the saturation region for 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. 7. 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) low-energy excitations, yet also suggest efficient computational complexity of thermal simulations for the specific model systems considered.

### iv.2 Universal entanglement behavior in (1+1)D conformal thermal states

In Fig. 8 we plot the entanglement entropy in a spin-1/2 Heisenberg chain. By simulating twice the length () as in Fig. 7, we can observe a logarithmic divergence in the low temperature region () The logarithmic entropy was already observed in the past and related to the computational complexity of finite temperature simulations Prosen and Pižorn (2007); Žnidarič et al. (2008). More recently, this was further analyzed by numerical simulations on Renyi entropy and also conformal field theory (CFT) analysis Barthel (2017). Notably, the finite-temperature entanglement calculation in Fig. 8 provides a convenient and accurate way to extract the central charge of CFT: without going into ground states calculation at Holzhey et al. (1994); P. and J. (2009), one can fit the MPO entanglement at finite temperatures. In Fig. 8, we observe that, by fitting data, the estimate of central charge is already very accurate (). For this we fit the data to the CFT prediction Barthel (2017). 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. 8 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. 8, the MPO entanglement shows a power-law behavior for . The slope on the log-log plot at very high temperatures is analyzed in the right inset in Fig. 8, which suggests a power-law exponent for . The growth in the entanglement , however, slows down strongly for , i.e., , where drops significantly below 1 as seen in the right inset of Fig. 8. The entanglement behaviors at high temperatures can be understood from a lowest, i.e., first-order expansion of density operator, . The singular value spectrum of this MPO is given by the vector , with another numerical vector. The resulting normalized “density matrix” of the supervector in Eq. (9) has eigenvalues with lowest-order thermal contributions . With the von Neumann entropy and

 γ(τ) ≡ dlnSdlnτ, (10a) one obtains that γ0 ≡ limτ→0+γ(τ)=2. (10b)

For simplicity, one may consider a single value for the vector , resulting in the two normalized weights

 (r1,r2)=11+α2τ2(1,α2τ2) (11a) with von Neumann entropy S(τ;α)=−2∑i=1rilnri. (11b)

A subsequent one-parameter fitting of with respect to to the actual entanglement entropy for nicely reproduces the high- entanglement data, as shown in the right inset in Fig. 8 for . The slope decreases monotonously, starting from and undergoing a sharp decrease around . Note, however, that the convergence towards the power-law exponent of for small is extremely slow, as also clearly supported by the simple asymptotic analysis above. For example, for , one only has .

Nevertheless, it follows from the generality of the above asymptotic argument, that the exponent 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 only holds for where , 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 .

### iv.3 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 Prosen and Pižorn (2007); Žnidarič et al. (2008), rather than exponentially as originally estimated Hastings (2006). 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. 9 we plot the entanglement property versus inverse temperature , for the 2D Heisenberg model on the square lattices of various system sizes, ranging from width to . As shown in the inset of Fig. 9, in the high temperature regime one still recovers the universal power-law , c.f. Eq. (10). Specifically, the data is very well fitted by the function in Eq. (11) with exactly the same parameter as in Fig. 8 for the 1D quantum chain.

For the inset of Fig. 9, since the entanglement entropy satisfies area law, i.e., , we divide the entropy by the width . For high temperatures , 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, , deviations from the strict scaling collapse of the area law can be observed.

In the main panel of Fig. 9, the entropy changes gradually into a logarithmic divergence versus as 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 is separated roughly by equal vertical offsets when incrementing the width for ). Since the calculations of the entropy in the system center are not fully converged for the wider systems, we also extrapolate the width and systems in . This actually further reinforces the regime of logarithmic increase of vs. for .

A similar additive logarithmic scaling of the entanglement entropy () in 2D Heisenberg model has been also found numerically via QMC calculations in the ground state Kallin et al. (2011); Song et al. (2011); Humeniuk and Roscilde (2012); Kulchytskyy et al. (2015); Luitz et al. (2015); Laflorencie et al. (2015). The coefficient of logarithmic correction was argued to be universal Metlitski and Grover (2011) (proportional to the number of Goldstone modes in the system). In the 2D Heisenberg model here, according to studies of ED and DMRG studies of the tower of states (ToS) in the energy and the entanglement spectra Alba et al. (2013); Kolley et al. (2013), respectively, the relevant low-energy ToS has characteristic level spacing that scales as with 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., .

These ToS excitations are responsible for the logarithmic entanglement at , and possibly also relate to the scaling of the entropy observed in the present study. In Fig. 9, we restrict the length to be as small as , which suggests that the magnon excitations are gapped out at temperatures as low as – 1/10. Therefore the relevant energy scale are likely only ToS modes with energy level spacing , which is smaller than the temperature in the regime where logarithmic scaling is observed.

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

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 high-temperature behaviors. We have also investigated the temperature dependence of entanglement properties in the MPO, and observed logarithmic entropies with constants and at low-temperatures not only in gapless quantum chains, but also in the square lattice Heisenberg model at fixed system size due to gapless ToS modes.

XTRG can be straightforwardly generalized to challenging 2D systems, such as frustrated antiferromagnets like the kagome and the square lattice Heisenberg models, as well as interacting fermionic models. The present MPO algorithms may be improved in several directions, including combining them with METTS samplings at low , or linked-cluster expansion to reduce finite-size effects, etc., which certainly deserves further exploration.

###### Acknowledgements.
Acknowledgments.— The authors thank Benedikt Bruognolo for nicely providing the METTS data and for very constructive discussions. BBC and WL also would like to acknowledge Han Li and Yun-Jing Liu for helpful discussions on related topics. This work was supported by the National Natural Science Foundation of China (Grant No. 11504014, 11474015, 61227902, 11774018) and the Beijing Key Discipline Foundation of Condensed Matter Physics. BBC, LC, and WL thank the hospitality of Arnold Sommerfeld Center for Theoretical Physics, University of Munich, where this work was finished. A.W. acknowledges support from the German Research Foundation (DFG) WE4819/2-1 and WE4819/3-1.

## Appendix A 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 Weichselbaum (2012b) that can deal with abelian and non-abelian symmetries such as SU() or the symplectic symmetry Sp() on a generic footing.

In order to emphasize the generality of the argument, we will frequently use the notation for a label of a generic irreducible representation (irep) of a given symmetry. Here, specifically, it may either stand for the spin-projection or spin label in the case of a U(1) [an SU(2)] symmetry, respectively. For this reason, we also refer to 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 of some given irep is defined by the unique representation that allows to form a scalar, i.e., with Clebsch-Gordan coefficients (CGCs) . These CGCs when properly normalized define a unitary matrix . This will be referred to as -symbol by analogy e.g., to 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., , abelian U(1) symmetries are not, since one has such that properly adds up to the scalar representation . For self-dual symmetries, such as SU(2), however, one must be careful in that , when written simply reduced to a unitary matrix of rank-2, it becomes indistinguishable from which, however, may differ by a sign, e.g., for half-integer spins in the case of SU(2). Importantly, -symbols allow to revert arrows in lines in a tensor network by inserting .

### a.1 Automata approach

Firstly, we briefly recapitulate the automata approach for constructing MPOs of the Hamiltonian Fröwis et al. (2010); Pirvu et al. (2010); Crosswhite and Bacon (2008). Consider, for example, the quantum Ising chain , 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 ). We need to compute and store the matrix elements of the spin operators {,}. Together with the identity operator, , these form a basis of local operators that enter a rank-4 tensor , where by rank we refer to the number of indices (or legs in a graphical depiction) of a given tensor. The tensor 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 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 -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 contains matrix elements, where is the dimension of the local state space , and is the bond dimension of the virtual state space . Every matrix element of 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 {,,}.

The virtual bond state space is given by a start state [, or equivalently, ], a stop state [, or equivalently, ], followed by which assigns an index position to every one of the interaction terms in the Hamiltonian that stretches across a given bond in the MPO (strictly speaking, 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 . For example, for the Ising model above, the -th bond in the system in between sites and carries the single interaction term , hence and .

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 (). This is carried through the MPO (therefore ) until an interaction term in the Hamiltonian occurs, say at site , which brings the automaton into the state (therefore . Having only nearest neighbor terms, the subsequent tensor e.g., at site immediately brings down the automaton to the end state (therefore ). By having completed the interaction term, the automaton stays in that state (hence ). 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 , is described by , which directly brings the automaton from the start into the end state. By translational invariance, there is nothing special about site , though. Therefore all of the matrix elements of the tensor 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.

### a.2 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 , let us continue with the model system of interest in this work, the (anisotropic) Heisenberg model

 ^H=∑⟨i,j⟩^Sxi^Sxj+^Syi^Syj=12(^S+i^S−j+^S−i^S+j)+J′^Szi^Szj , (12)

where we temporarily introduce hats on top of operators, in order to differentiate them from symmetry labels (e.g., vs. ). The model in Eq. (12) is U(1) symmetric as it preserves . In the isotropic case, , it becomes SU(2) spin symmetric. Then the spin operators need to be grouped into a spinor,

 ^S≡⎛⎜ ⎜ ⎜ ⎜⎝−1√2^S+^Sz+1√2^S−⎞⎟ ⎟ ⎟ ⎟⎠ (13)

such that Eq. (12) can be rewritten in SU(2) invariant form

 ^H=∑⟨i,j⟩^S†i⋅^Sj . (14)

Note that the relative weights and signs in Eq. (13) are important for consistency with standard conventions on SU(2) spin multiplets. In particular, the operators in the spinor in Eq. (13) exactly represent, top to bottom, the states of an spin multiplet (e.g., see Ref. Weichselbaum (2012b)). In contrast, the Hermitian set of operators does not. However, using Eq. (13), the dagger on one of the spinors in Eq. (14) 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 [cf. Fig. A.1(a)]. This can be written as the irreducible operator (IROP) where the composite index naturally specifies entire state spaces Weichselbaum (2012b), or here an operator space: the index differentiates between local IROPs that transform according to the same irreducible representation . By definition of an index, , we therefore also introduces an arbitrary but fixed order to the local operators. The label , finally, fully differentiates the operators within a given spinor Weichselbaum (2012b). For example, within SU(2), simply stands for .

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 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 operator . In contrast, the operators carry , hence switch between symmetry sectors, and therefore are not considered scalar operators.

All local operators eventually can be combined into a single rank-3 tensor [cf. Figs. A.1(a,b)]. The third index then represents the state space Weichselbaum (2012b) of the “supervectors” . For efficiency, the set of local operators should be orthogonal in the sense

 tr[(X[nq;qz])†X[n′q′;q′z]]∝δnn′δq,q′δqz,q′z , (15)

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 of the MPO is given, note that the intermediate supervector index that connects the super-MPS with the local operators [cf. Fig. A.1(c)] may also be generated by the reverse operation of splitting off the local state space from the tensor via SVD. Then by construction, the operators in 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 Fig. A.1(a-b)]. By finally contracting the super-MPS (rank-3 tensors) with the local operators along the intermediate index, this leads to the final rank-4 tensors that constitutes the MPO [cf. Figs. A.1(c,d)]. Note that the intermediate index also specifies an arbitrary but fixed order of the local operators (‘supervectors”).

Now the structure of an interaction term as in Eq. (14) is generic: a non-scalar irreducible operator acting on site must be paired up, i.e., contracted on the spinor index into a scalar term of the Hamiltonian with another operator acting on site that transforms according to exactly the same irreducible representation (typically ; 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. (12) or Eq. (14). In partiuclar, the super-MPS can be built completely analogous to the automata approach above, while paying simple attention to symmetry sectors. The start () and the stop () state on the virtual bonds transform like scalars (i.e., have ), whereas the bond states 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. (12), the set of local operators is given by for the U(1) symmetric setup, and by for the SU(2) symmetric setup. In either case, the set of local operators is orthogonal as in Eq. (15). 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 . Specifically, with a unitary transformation that corresponds to the Clebsch-Gordan coefficients which combine a spin with its dual (again ) into a singlet [cf. -symbol earlier], the spin-spin interaction can be written as . Therefore the action of the dagger on one of the spin operators can be transferred via the unitary -symbol into the definition of the super-MPS itself, proper sign-convention on 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 increases [see also Fig. A.1(e)]. The period of the translational invariance of the MPO also increases from to the width of the system and hence requires at least different -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. (15)]. If the bond-dimension 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 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 singular value decomposition (SVD) or variational techniques. The latter is unavoidable for MPO products or sums in any case, as will be discussed next.

## Appendix B Initialization in the XTRG algorithm

In this section, we compare three different initializations of in the XTRG algorithm. The quality of our initial for small is measured by estimating the relative error of the free energy , starting from exponentially small (i.e., the first data point after initialization at ) down to intermediate temperatures .

### b.1 Series expansion vs. Trotter-Suzuki initialization

Firstly, we compare the Trotter-Suzuki initialization with series expansion [Eq. (6a) in the main tex] followed by XTRG. Trotter-Suzuki decomposition breaks into product of local evolution gates, which for nearest-neighbor spin-1/2 chains within first-order can be represented as an MPO with bond dimension (, comprised of . In Fig. A.2 we plot the relative errors of the free energy after Trotter-Suzuki initialization at two values of and . At , the respective errors are and , respectively.

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 (; green data). The resulting errors are much lower and for and 0.01, respectively.

By increasing the initial bond dimension to and slightly above, i.e., (), as seen in Fig. A.2, SETTN initialization as in Eq. (6a) 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 , the SETTN initialization with is still an order of magnitude more accurate as compared to the Trotter-initialized data, and orders of magnitude more accurate if is only marginally increased to relative to the of the subsequent XTRG.

### b.2 Linear initialization ρ(τ0)≃I−τ0H

By sweeping over several orders of energy scales, the XTRG algorithm permits a simple linear initialization