Thermal pure matrix product state in two dimensions: tracking thermal equilibrium from paramagnet down to the Kitaev honeycomb spin liquid state

We present the first successful application of the matrix product state (MPS) representing a thermal quantum pure state (TPQ) in equilibrium in two spatial dimensions over almost the entire temperature range. We use the Kitaev honeycomb model as a prominent example hosting a quantum spin liquid (QSL) ground state to target the two specific-heat peaks previously solved nearly exactly using the free Majorana fermionic description. Starting from the high-temperature random state, our TPQ-MPS framework on a cylinder precisely reproduces these peaks, showing that the quantum many-body description based on spins can still capture the emergent itinerant Majorana fermions in a ${\mathbb Z}_2$ gauge field. The truncation process efficiently discards the high-energy states, eventually reaching the long-range entangled topological state approaching the exact ground state for a given finite size cluster. An advantage of TPQ-MPS over exact diagonalization or purification-based methods is its lowered numerical cost coming from a reduced effective Hilbert space even at finite temperature.


Introduction
Characterizing a thermal quantum state, a quantum many-body state at finite temperature is an ongoing fundamental challenge in condensed matter physics and beyond, since it is often a matter of quantum and classical correlations studied in statistical and quantum information physics [1].Such a state has an intriguing aspect in that its representation is largely left facultative [2]; the Gibbs state is a mixture of an exponential number of states given by the density matrix ρ β of small purity P ∼ e −Θ(N ) , i.e., vanishing exponentially with the system size N .The thermal pure quantum (TPQ) state, on the other hand, is a single pure state of purity P = 1.In addition, there exist numerous thermal mixed quantum (TMQ) states with a purity between Gibbs and TPQ (see Fig. 1(a)).Canonical typicality guarantees that all these choices equivalently yield the same thermal equilibrium properties of the subsystem [3,4], and are macroscopically in the "same" thermal state.Since Gibbs, TPQ, and TMQ states rely on different design concepts, even when applying the "same" tensor network representation, its structure, convergence, or the amount of numerical resources required likely depend on which type of thermal state is chosen.
An important development concerning the Gibbs state is the matrix product density operator (MPDO), which provides a direct tensor network representation of the density matrix operator, ρ β [5,6].Another standard form of the Gibbs state is the purified state analog to thermofield double, consisting of the size-N system and the same numbers of ancilla degrees of freedom each suspended to a local site [7].Ancilla serve as an entanglement bath and tracing out the ancilla corresponds to taking the Gibbs ensemble.These doubled states also conform to a matrix product operator (MPO) approach, 1 whose schematic illustrations are shown in Fig. 1(a).Here, the entanglement entropy is meaningless as a measure to characterize the Gibbs state.Instead, the thermal area law of mutual information between subsystems determines the bond dimension χ of MPO's [10][11][12].The numerical drawback of MPDO or purification is the increase of the Hilbert space dimension due to the doubled degrees of freedom.Still, MPDO has been developed further recently using the XTRG algorithm [13], which realizes an exponential cooling down of the system by iteratively multiplying the matrix ρ β × ρ β = ρ 2β , allowing to reach very low temperatures rapidly.XTRG has successfully been applied to two dimensions [14,15] including our target, the Kitaev honeycomb model [16].
The TPQ state, in comparison, consisting only of physical degrees of freedom, is pure by construction, and does not need the doubling of the local Hilbert space.In MPDO and its analogues, the doubling or the ancilla play the role of an ensemble average-or the classical mixture of states-which provide the volume-law thermal entropy.The lack of doubling implies that the pure TPQ state needs to store the same amount of entropy internally as a volume-law entanglement entropy [17][18][19].For such purpose, the tensor-network-based representation bounded by the area law entanglement are thought to naturally be out of reach.Yet, the authors have recently exploited the specific form of matrix product state (MPS) practically recovering the volume law entanglement; only two ancilla/auxiliaries are attached to both edges of the one-dimensional (1D) MPS train, yet they have turned out to be sufficient to keep the nearly uniform distribution of entanglement entropy density throughout the system 2 which 1 The difference between MPDO and purification is that the MPDO is not necessarily positive definite after truncation, whereas purification using a canonical form is positive definite.However, purification generally requires larger χ than MPDO [8], and there are some examples [9] that the purification MPO shows a divergence of χ at low temperatures, which may indicate that the thermal area law may not safely apply. 2 If we take a bipartition of the TPQ-MPS system into left and right, each attached to the auxiliary, the entanglement entropy does not depend on the size of the left/right part, unlike the usual MPS that follows the size-dependent Page curve.This translational invariance of the entanglement entropy allows entanglement entropy between the center-n sites and the rest (with N − n sites and two auxiliaries) to follow the n-linear volume law (see Ref. [19]).  is essential for the volume law entanglement.We call this construction the TPQ-MPS [19].
The TPQ state itself has a numerically long history [20][21][22][23] far before the formulative seminal works [24,25].They mostly rely on a full Hilbert space representation using Lanczos-based methods that limit the system size to typically N ≲ 30 − 40.The TPQ-MPS largely shrinks the representation space and increases N by factors by efficiently choosing its constituent states to those representing the target temperature limited by the bond dimension of the MPS.We review a measure of the quality of a TPQ-MPS, which has been developed in Ref. [2], in Appendix A.
The present work advances a few steps in developing a TPQ-MPS for two dimensions (2D), particularly for a quantum mechanically nontrivial quantum spin liquid state with long-range entanglement.Encoding the substantial amount of entanglement expected for QSL within an MPS or a tensor-network is generically a challenging task, although reported in the case of ground state [26][27][28].Our result is the first to track the state by an MPS in the nearly pure form from the high-temperature random state down to the QSL with substantial entanglement between limited selection of basis states.
We finally refer to some TMQ-state-based approaches; the minimally entangled typical thermal state (METTS) [29,30] mixes (takes an equal weight average of) a series of MPS generated from the Markov process.The quantum Monte Carlo designs a local product state basis to suppress the sign problem [31,32], which are recently highlighted in combination with the iPEPS.

Construction of the TPQ-MPS state
We consider the standard imaginary-time evolution in generating the TPQ state at inverse temperature β = 1/T given as where H is the Hamiltonian of the system of interest, and the initial state |Ψ 0 〉 representing an 'infinite-T ' state is chosen as random, satisfying |Ψ 0 〉〈Ψ 0 | ∝ I, where • • • is the random average and I is the unit matrix.We now specify the construction of TPQ-MPS utilized here.The 1D tensor train of size-N and bond dimension χ is prepared with auxiliary degrees of freedom added to both ends to provide an entanglement bath (Fig. 1(b)).Here, instead of the χ × χ form proposed in Ref. [19], each auxiliary consists of N aux sites with the same local Hilbert space d as the physical sites of the system, i.e. d = 2 for S = 1 2 spins, resulting in rank-3 tensors of the form The number of auxiliary sites dictates the maximum bond dimension at the edge of the physical system as χ aux = d N aux and, hence, the maximum amount of entanglement between the auxiliary and the system. 3We emphasize that the auxiliary sites are not coupled to the physical system by any physical exchange, and therefore only the identity is applied to them during the imaginary time-evolution.
We extend TPQ-MPS to two spatial dimensions by wrapping the lattice on a cylinder with a finite circumference and wind the 1D MPS structure around, enumerating all the sites linearly (see Fig. 1(c)).Cylinder tensor networks are fairly standard techniques nowadays, involving various variants in the way of wrapping the lattice and subsequent enumeration schemes.The precise way of wrapping the lattice can have physical implications; The system, although gapless in the two-dimensional limit, maybe gapped if the gapless nodes are not on allowed momenta lines in the Brillouin zone [27].There are choices of particular cylindrical geometry known to capture the gapless state of the KH model [28,33], but are not used here.The choice of such cylinder is important for the ground state but not for the temperature we can reach in the present study.The enumeration scheme should, ideally, not alter the physical properties.However, in reality, it can influence the spatial distribution of correlations and entanglement in particular for relatively small bond dimensions [14].We employ a helical enumeration scheme with 8 × 3 × 2 (YC8×3×2, which has circumference L circ = 6 and is illustrated in Fig. 1(c)) and 8 × 4 × 2 sites (YC8×3×2, L circ = 8) conforming to YC3-1 and YC4-1, respectively, using the convention in Ref. [34].Both schemes treat the xand z-bond on equal footing, i.e. they are nearest neighbors in the 1D MPS structure, while the y-bonds turn into an exchange with range 2L circ − 1 sites.This choice results in the smallest χ M PO of the time-evolution unitary, while also reducing the number of nearest-neighbor bonds cut by a bipartition which, at sufficiently low T , enters the amount of entanglement entropy encoded in the TPQ-MPS.
The long-range interactions within the effective 1D model make the time-evolving block decimation scheme [35][36][37] in Eq. ( 1) infeasible.Instead, we rely on an MPO formulation of the time-evolution operator [38]. 4Specifically, we discretise U(β/2) = [U (dτ)] N with small imaginary time steps dτ and represent U(dτ) as MPO [38]. 5After each MPO-MPS product, the MPS is compressed using a variational scheme [41] reducing χ.We use an upper bound for the maximum χ = {362, 512, 724, 1024} to limit the computational resources needed.If the bound is not reached, small Schmidt values λ i are discarded provided either of the two criteria are met: (I) discard all λ i ≤ s trunc or (II) discard all λ i sufficing i λ 2 i < (s sum trunc ) 2 beginning from the smallest λ i .
Further technical details are given as follows; The initial random TPQ-MPS state |Ψ 0 〉 is prepared by applying a sequence of random two-site unitary matrices to a Néel-like product state in the z-basis, i.e. | • • • ↑↓ • • • 〉, in a TEBD-like way.We prepare N samples = 100 independent random initial states using 25 TEBD-iterations and cap the bond dimension at χ ini = 32.See Appendix A and Ref. [2] for further details regarding the random initial state.The imaginarytime step is chosen as dτ = 0.1 and smaller for β ≤ 0.8.Truncation thresholds are set to s trunc = 10 −6 and s sum trunc = 10 −5 unless stated otherwise.Measurements are not independent concerning β, since they are taken in a sequence of fixed β during each independent run of imaginary-time evolution, while the N samples averages are taken from the set of independent runs.The TenPy library [42] is used for all MPS-related numerical calculations.

Application to the Kitaev honeycomb model
We employ TPQ-MPS to the Kitaev honeycomb (KH) model defined as [16] where σ γ i are Pauli operators σ x , σ y , and σ z acting on sites i.The three sets of parallel bonds on the honeycomb lattice are labeled as γ = {x, y, z} (see Fig. 1(c)).The Kitaev interaction K γ couples a neighboring pair of spins 〈i, j〉 γ along the γ-bond by an Ising-like exchange The KH model features a gapless QSL ground state if K α ≤ K β + K γ is satisfied for all permutations of the bond labels {x, y, z}.Otherwise, a gapped QSL is found which adiabatically connects to the Toric Code [43].Here, we focus on the case of K x = K y = K z = 1  3 .The KH model features a double-peak structure in the specific heat, signalling crossovers and releasing an entropy of ∆S/N = 1 2 ln 2 each.The associated two energy scales are well known [44]: At the high-T peak, T H /K ≈ 0.5, nearest-neighbor spin-spin correlations develop and the fractionalization into itinerant and localized Majorana fermions occurs.The latter contributes to the formation of fluxes at each hexagonal plaquette given as W P = i∈P σ γ P (i) i , where γ P (i) = x, y, z is the label of bond connected to site i while not being part of the plaquette P. The fluxes give an extensive set of quantum numbers, w P = ±1, which are disordered at T ≲ T H /K. Below the low-T peak, T L /K ≈ 0.016, the fluctuation of fluxes is suppressed and we eventually find 〈W P 〉 → 1.They form the static 2 lattice-gauge field, fixing half of the Hilbert space per unit cell.A local Hilbert space dimension of 2 per site remains which is associated with itinerant Majorana fermions.Although the KH model at finite temperature YC8×4×2 exhibits a good quantitative agreement with MC+FFED down to T ∼ 0.05.We attribute the difference in the position of T L , in particular for YC8×3×2, to the finite circumference geometry used here; the ground state energies E GS of an infinitely long L circ = 6 cylinder (YC3-1, dash-dotted horizontal line) obtained by iDMRG has a ground state energy density lower than the bulk exact one [16] (solid horizontal line) by about ∼ 0.01K, yielding different crossover slopes in E near T L and a shifted peak.The XTRG result using a 6 × 4 × 2 cylinder is extracted from Fig. 4 in Ref. [15] and shown for comparison.The inset focuses on T ≤ 0.42 using a linear scale for the temperature.(c) The evolution of the plaquette flux average, 〈W p 〉.The flux-free state at T = 0 exhibits 〈W p 〉 = 1.
is not exactly solvable, once 2 bond variables constituting the 2 gauge field are treated as classical degrees of freedom, a combination of classical Monte Carlo method with free (Majorana) fermion exact diagonalization (MC+FFED) provides a nearly exact calculation in a relatively large cluster, as performed by Nasu, et.al [44].Whereas, its counterpart Eq.( 2) is a quantum many-body Hamiltonian which is generically difficult to solve at finite temperatures straightforwardly by an unbiased quantum many-body calculation.Therefore, the model provides a good platform and benchmark for our approach.We would like to emphasize that our approach, unlike MC+FFED, is not custom tailored to the KH model and can be applied to other quantum many-body Hamiltonian.
Our TPQ-MPS data in Fig. 2 exhibits a good qualitative agreement with the results obtained from MC+FFED [44] on a 12 × 12 × 2 cluster and XTRG using a 6 × 4 × 2 geometry [15]; The energy density6 E rapidly decreases near T H resulting in a crossover peak in the specific heat C.
A second step of energy reduction occurs near T L .The two-step behavior is already present for small χ = 362 with well converged behaviour down to T ∼ 0.2 including the high-T peak in the specific heat.Whereas for T ≲ 0.2, the finite-size and finite-χ effects inevitably influence the data; In Fig. 2(a) we display in two different lines the ground state energy obtained using iDMRG on an infinite cylinder with the same circumference L circ = 6, 8 and helical boundary condition YC3-1 and YC4-1, respectively.The circumference seriously affect the numerically achieved ground state energies and consequently the specific heat which can be summarized as follows: (I) The cylinder with L circ = 6 features an enhanced reduction in energy upon cooling down approaching the significantly lower ground state energy.The low-T peak in specific heat is of similar height to MC+FFED, while shifted to a two to three times higher temperature.(II) For L circ = 8 we obtain an evolution of the energy closer to MC+FFED, thus reducing the finite-size effect signicantly.For T L ≳ T ≳ 0.2, however, TPQ-MPS overestimates E compared to MC+FFED.Here, increasing χ gradually reduces E possibly approaching MC+FFED for sufficiently large χ.Near T ∼ T L and below, the effect of finite χ ceases and the energy eventually approaches both MC+FFED as well as the ground state energy.As a consequence of the overestimated energy density at intermediate T , we obtain an enhanced slope of E resulting in a higher peak in the specific heat.Again, increasing χ improves accuracy, reduces the height of the peak and results in a behaviour closer to MC+FFED.The peak position is very similar to MC+FFED at any χ.
The average of 2 fluxes in Fig. 2(c) nicely marks the two peaks by an onset of nonzero value (T H ) and the inflection point (T L ), finally approaching 〈W P 〉 → 1 at T → 0 systematically for various χ.
A recent XTRG calculation applied to the KH model reports the lower-T -peak at T L ∼ 0.023 with the peak-height of ∼ 0.3 using a 6 × 4 × 2 cylinder [15].While the circumference is similar to our YC8×4×2, the XTRG work uses a slightly shorter cylinder, does not use helical boundary condition, and employs a different winding scheme.The quantitative agreement of XTRG with MC+FFED and TPQ-MPS is very good above T ∼ 0.1 where finite-size effects become negligible.At lower temperature, however, deviations become apparent (see Fig. 2): Our geometry YC8×4×2 with helical boundary condition exhibits a ground state energy close to the thermodynamic limit, whereas XTRG uses a different winding scheme, which influences the location of T L .The comparison with our two geometries confirms that the size or shape of the cylinder shifts the T L peak.The height of the peak in XTRG is similar to that of MC+FFED and YC8×3×2.Both TPQ-MPS and XTRG give reasonable results for the given finite size system, but the choice of the cylinder can easily influence the quantitative agreement of the data against the bulk data at T ≤ 0.1.
In this context, we like to remark that in many frustrated spin models, the specific heat at T ≲ 0.1 naturally suffers large finite-size effect independent of the method employed.For example, in Kagome-lattice Heisenberg antiferromagnet, specific choices of clusters sometimes yield unphysical peaks or features not observed in other choices of cluster [25,45] possibly obscuring the physical behaviour.

How truncation affects the TPQ-MPS state
We now quantify the TPQ-MPS based on the error analysis during the run by focusing on two quantities: The first one is the sum of all discarded Schmidt values (λ i for i > i 0 (β j ) which fulfills aforementioned I or II in Sec. 2) accumulated over a single imaginary-time evolution, The second one is the product of the fidelities of the state |Ψ(β j )〉 = W II (dτ)|Ψ(β j−1 )〉 (see Ref. [36]) and | Ψ(β j )〉 just before and after truncation, respectively, for all truncations down to the temperature which evaluates how we deviate from the non-truncated wave function at β.The amount of truncated Schmidt values per unit of imaginary time is given as ∂ β Σ trunc .In Figure 3 we show the evolution of ∂ β Σ trunc , of the average bond dimension χ, and of F. Upon cooling down, ∂ β Σ trunc remains below 10 −7 until χ reaches the upper bound χ ∼ χ, which occurs near T ∼ 1. Larger χ or smaller system size generally lowers this threshold temperature.At these high temperatures, the evolution is very accurate reflected in a fidelity F ∼ 1. Upon lowering the temperature, ∂ β Σ trunc increases and then reaches a plateau at T L ≲ T ≲ T H with values ∂ β Σ trunc ∼ 10 −5 to 10 −4 depending on χ.Here, F starts to depart gradually from 1, which is more distinct for smaller χ.At T ≲ T L the error ∂ β Σ trunc reduces again and F starts to flatten out.In particular for YC8×3×2, a drop in χ is apparent, indicating the reduction in the size of the Hilbert space needed to effectively encode the low-temperature state.These observations suggest two effects of the truncation χ; For relatively small χ that is reached quickly, in particular at intermediate T L ≲ T ≲ T H , taking a larger χ lowers the energy towards the optimal value.This becomes evident upon inspection of E in Fig. 2(a), whose accuracy improves for larger χ approaching the MC+FFED data.
The second effect concerns the states at high energy.Let us expand the TPQ state constructed for the full Hilbert space for finite N .The system is split into a smaller part A (with dimension D A ) and a bigger part B, which is Schmidt decomposed as to the orthogonal basis sets {|n A 〉} and {|n B 〉}.The local part A is thermalized and its density operator is approximated by the Gibbs state in A as where {|n A 〉} is thought to be the energy eigenbasis of the subsystem's Hamiltonian H A .For its eigenvalues {E A n }, the Schmidt coefficient λ n is represented as e −β E A n /2 / Z A , and we find Note here that {|n B 〉} is left unknown.We finally truncate |Ψ β 〉 as D A → χ in Eq.( 7), discarding the basis states with small weight.Specifically, information of |n A 〉 belonging to higher E A n is lost.This explains the capability of TPQ-MPS to express qualitatively different quantum states from high to low temperatures; The truncation of the MPS efficiently compresses the information needed to represent the thermal state particular at low temperatures.
The above context of discarding high-temperature states-or high-energy states, respectively-explains a particular feature of TPQ-MPS: the variance of physical quantities among different initial states becomes smaller by more than one order for lower temperature [2,19].This tendency is opposite to the usual random sampling methods including standard TPQ or Monte Carlo methods, where the sampling error is by orders of magnitude larger in the lower temperature phase.We illustrate this point further by referring to the standard TPQ using the full Hilbert space for a limited system size N : It is shown in Ref. [25] that the variance increases at low temperatures by the order of ≤ e −S N /T , and when the entropy S N of size N is sufficiently large, the increase is moderately suppressed.This fact supports the application of TPQ methods to highly frustrated quantum magnets including the present KH model and Kagome or related lattice models [46][47][48][49][50][51].However, the idea of relying on the large entropy does not apply to TPQ-MPS: in the first proposal of TPQ-MPS in Ref. [19] some of the authors have shown that even for non-frustrated systems, the sample variance becomes smaller at lower temperature contrary to the prospect from TPQ.This is intuitively because MPS provides a good description of a quantum many body state at zero temperature.Our Eq. ( 7) and the related discussions support this observation irrespective of the choice of spatial dimensions, and suggest good applicability of the present 2D TPQ-MPS to non-frustrated models.

Conclusion
To summarize, the TPQ-MPS is applied to 2D by wrapping the MPS train into cylinders.The two peaks in the specific heat of the KH model signaling the fractionalization of spins into Majorana fermions and fixing the 2 gauge flux are both reproduced.While finite-size effects appear at T ≲ 0.1 as is common with other methods, finite-χ affects the MPS-TPQ only at intermediate temperatures T ∼ 0.1 and is less of a concern at very low temperatures T ≲ 0.01.This fact is in sharp contrast to other random sampling methods including the original TPQ method using the full Hilbert space.Here, the truncation process of TPQ-MPS efficiently discards the higher-temperature information explaining why it can track a nearly pure thermal state with its volume-law entanglement equivalent to the thermal entropy across a wide range of temperatures.This allows the state starting from random at high temperature (initial state) to gradually reach the qualitatively different long-range entangled topological ordered ground state.The application to non-frustrated model is expected to be promising because, unlike for the original TPQ, a high entropy density is not required at low temperature to attain a reasonable accuracy at a moderate numerical cost.where • • • is the random average.It is shown that the purity of the thermal state scales with δz 2 and the larger δz 2 means that the obtained state varies much with a sample.In fact, we showed that the number of samples needed, N sample , to obtain the same quality of Eq.(A.1) increases proportionally to δz 2 [2].
The results of δz 2 for the present calculation on the KH model are given in Fig. 4(a) for a set of data given in Fig. 2. The largest δz 2 at low-T ranges at 10 −1 − 10 1 , which is comparable to the value for the 1D Heisenberg model [2] using N sample = 100.Based on this comparison, we also adopt N sample = 100 for the KH model.The present calculation shows that the 2D TPQ-MPS is as capable as the 1D case despite the consensus that the calculations in 2D are much more difficult than in 1D.
The plateau of ∂ β Σ trunc observed in Fig. 3 agrees with the plateau of δz 2 , and as in Fig. 3, χ dependence appears at T ≲ 10 0 , showing that δz 2 is indeed a good measure to qualify the quantum state.We find a suppression of δz 2 by a log scale in terms of χ, indicating the high capability of TPQ-MPS to store the information required for a wide range of temperatures exhibiting different natures.
As mentioned above, TPQ-MPS stores substantial amount of entropy in the starting point of the imaginary time evolution, which is contrary to METTS.Therefore, the quality of the initial random state is important to have high purity and smaller δz 2 , respectively.When generating the initial random TPQ-MPS state, we use a TEBD-like algorithm with alternating application of random two-site unitary matrices.Although not relevant to our model, an advantage of this method is the possibility of utilizing charge or S z conservation, e.g. when studying a U(1) symmetric model.The iteration number, we use 25 iterations, is determined by ensuring a saturated entanglement entropy.Figure 4(b) shows the evolution of entanglement entropy of bipartition SE , averaged over the physical sites, as a function of iterations.While the bond dimension of the state generally doubles after each iterations, SE is not fully saturated after

Figure 1 :
Figure 1: (a) Schematic illustration of different thermal quantum states, classified by purity.The two well-known representations of the Gibbs state are sketched.Purification prepares a product state of singlets on a pair of system and ancilla sites, and operates ρ B ⊗ U aux where ρ B = e −β Ĥ/2 acts on the system and U aux on the ancilla.The resultant state is shown in the simplified tensor form.(b) Schematic illustration of TPQ-MPS with auxiliary degrees of freedom providing an entanglement bath, and the MPO-based imaginary time evolution.(c) Illustration of the honeycomb lattice composed of x-, y-, and zbonds.For the underlying 1D MPS structure, we use a cylindrical geometry with a helical enumeration scheme as is highlighted by orange dashed lines.Equivalent bonds across the boundary are marked with the same roman literal.Specifically shown is the YC8×3×2 geometry with a shifted (by one lattice vector) boundary condition.Orange dots connected by solid lines represent the auxiliary sites.

Figure 2 :
Figure 2: Temperature (T = β −1 )-dependent (a) energy density E and (b) specific heat density exhibiting the double peak obtained by TPQ-MPS on cluster YC8×3×2 and YC8×4×2 and several upper bounds for χ.Reference data (black dots) uses MC+FFED on 12 × 12 × 2 sites[44] and was provided by courtesy of J. Nasu.YC8×4×2 exhibits a good quantitative agreement with MC+FFED down to T ∼ 0.05.We attribute the difference in the position of T L , in particular for YC8×3×2, to the finite circumference geometry used here; the ground state energies E GS of an infinitely long L circ = 6 cylinder (YC3-1, dash-dotted horizontal line) obtained by iDMRG has a ground state energy density lower than the bulk exact one[16] (solid horizontal line) by about ∼ 0.01K, yielding different crossover slopes in E near T L and a shifted peak.The XTRG result using a 6 × 4 × 2 cylinder is extracted from Fig.4in Ref.[15] and shown for comparison.The inset focuses on T ≤ 0.42 using a linear scale for the temperature.(c) The evolution of the plaquette flux average, 〈W p 〉.The flux-free state at T = 0 exhibits 〈W p 〉 = 1.

Figure 3 :
Figure 3: Evolution of the truncation error and bond dimension of the KH model: (a) Accumulated truncation error per unit of imaginary-time ∂ β Σ trunc and bond dimension χ averaged over the system, and (b) fidelity F of the evolved MPS before and after the truncation accumulated for all β j < β.

Figure 4 :
Figure 4: (a) Normalized fluctuation partition function (NFPF) δz 2 as a function of T for the data calculated in Fig. 2 for the KH model with different χ.The number of samples required to obtain the same quality data scales linearly with δz 2 .(b) Evolution of the entanglement entropy of bipartition as a function of iterations.Each iteration consists of applying two-site random unitary matrices on even and odd bonds similar to the TEBD algorithm.SE is obtained by taking the average of S E over bipartitions between any of the physical sites.An additional sample averaging is done over 10 samples to obtain error bars indicating the variance.However, error bars are typically smaller than the marker size.