SciPost Phys. 15, 206 (2023)

Thermal pure matrix product state in two dimensions: Tracking
thermal equilibrium from paramagnet down to the Kitaev
honeycomb spin liquid state
Matthias Gohlke1 , Atsushi Iwaki2 and Chisa Hotta2
1 Theory of Quantum Matter Unit, Okinawa Institute of Science andTechnology Graduate University, Onna-son, Okinawa 904-0495, Japan
2 Department of Basic Science, The University of Tokyo, Meguro-ku, Tokyo 153-8902, Japanand Komaba Institute for Science, The University of Tokyo,
Meguro-ku, Tokyo 153-8902, Japan

Abstract
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 specificheat peaks previously solved nearly exactly using the free Majorana fermionic description. Starting from the high-temperature random state, our TPQ-MPS framework on acylinder precisely reproduces these peaks, showing that the quantum many-body description based on spins can still capture the emergent itinerant Majorana fermions in a
Z2 gauge field. The truncation process efficiently discards the high-energy states, eventually reaching the long-range entangled topological state approaching the exact groundstate for a given finite size cluster. An advantage of TPQ-MPS over exact diagonalizationor purification-based methods is its lowered numerical cost coming from a reduced effective Hilbert space even at finite temperature.
Copyright M. Gohlke et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 03-08-2023
Accepted 10-11-2023
Check forupdatesPublished 24-11-2023
doi:10.21468/SciPostPhys.15.5.206

Contents
1 Introduction

2

2 Construction of the TPQ-MPS state

4

3 Application to the Kitaev honeycomb model

5

4 How truncation affects the TPQ-MPS state

7

5 Conclusion

9

A Random sampling average

10

References

121

SciPost Phys. 15, 206 (2023)

1

Introduction

Characterizing a thermal quantum state, a quantum many-body state at finite temperature isan 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 thedensity matrix ρβ of small purity P ∼ e−Θ(N ) , i.e., vanishing exponentially with the systemsize N . The thermal pure quantum (TPQ) state, on the other hand, is a single pure state ofpurity P = 1. In addition, there exist numerous thermal mixed quantum (TMQ) states with apurity between Gibbs and TPQ (see Fig. 1(a)). Canonical typicality guarantees that all thesechoices 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 relyon 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 tothermofield 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 tracingout the ancilla corresponds to taking the Gibbs ensemble. These doubled states also conformto 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–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 twodimensions [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–19]. For such purpose, the tensor-network-based representation
bounded by the area law entanglement are thought to naturally be out of reach. Yet, theauthors 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 edgesof the one-dimensional (1D) MPS train, yet they have turned out to be sufficient to keep thenearly uniform distribution of entanglement entropy density throughout the system2 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 thesize-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]).

2

SciPost Phys. 15, 206 (2023)

(a)

(a)
MPDO

...

...

...

QMC, random sampling

purification(b) aux χ
L
L
1
1

1
1

1

2

3

W

W

W

W

W

N −2N −1

W

W

W

W

W

χR auxR
N

W
W

1

1

1

1

|ψ(β = 0)iU(β/2) = WτN(c) YC8×3×2 geometry
I

auxL

5
3

1

2

4

VII

II

6

N

χR

y

···
7

χL

auxR
I

x

VI

z

VII

Figure 1: (a) Schematic illustration of different thermal quantum states, classified bypurity. 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, andpp
operates ρB ⊗ Uaux where ρB = e−β Ĥ/2 acts on the system and Uaux on the ancilla.
The resultant state is shown in the simplified tensor form. (b) Schematic illustrationof TPQ-MPS with auxiliary degrees of freedom providing an entanglement bath, andthe MPO-based imaginary time evolution. (c) Illustration of the honeycomb latticecomposed of x-, y-, and z- bonds. 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 sameroman 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.
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–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 therepresentation space and increases N by factors by efficiently choosing its constituent statesto 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
3

SciPost Phys. 15, 206 (2023)

MPS or a tensor-network is generically a challenging task, although reported in the case ofground state [26–28]. Our result is the first to track the state by an MPS in the nearly pureform 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 typicalthermal 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 statebasis to suppress the sign problem [31, 32], which are recently highlighted in combination
with the iPEPS.

2

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β|Ψβ 〉 = e− 2 H |Ψ0 〉 = U (β/2) |Ψ0 〉 ,(1)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 sizeN 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 inRef. [19], each auxiliary consists of Naux sites with the same local Hilbert space d as the physical sites of the system, i.e. d = 2 for S = 12 spins, resulting in rank-3 tensors of the form
χi−1 × χi × d. The number of auxiliary sites dictates the maximum bond dimension at the edgeof the physical system as χaux = d Naux and, hence, the maximum amount of entanglement
between the auxiliary and the system.3 We emphasize that the auxiliary sites are not coupled
to the physical system by any physical exchange, and therefore only the identity is applied tothem during the imaginary time-evolution.
We extend TPQ-MPS to two spatial dimensions by wrapping the lattice on a cylinder with afinite 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, althoughgapless 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 choiceof such cylinder is important for the ground state but not for the temperature we can reach inthe 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 inparticular for relatively small bond dimensions [14]. We employ a helical enumeration schemewith 8 × 3 × 2 (YC8×3×2, which has circumference Lcirc = 6 and is illustrated in Fig. 1(c)) and8 × 4 × 2 sites (YC8×3×2, Lcirc = 8) conforming to YC3-1 and YC4-1, respectively, using theconvention in Ref. [34]. Both schemes treat the x- and z-bond on equal footing, i.e. they arenearest neighbors in the 1D MPS structure, while the y-bonds turn into an exchange with range
2Lcirc − 1 sites. This choice results in the smallest χ M PO of the time-evolution unitary, while
3
N
Using Naux spins of dimension d is equivalent to preparing a single degree of freedom with daux
. However, theformer has practical advantage regarding the ease of implementation and the physical intuition about the degrees
of freedom included in the bath.

4

SciPost Phys. 15, 206 (2023)

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–37] in Eq. (1) infeasible. Instead, we rely on an MPO formulation of
the time-evolution operator [38].4 Specifically, we discretise U (β/2) = [U (dτ)]N with small
imaginary time steps dτ and represent U (dτ) as MPO [38].5 After each MPO-MPS product,
the MPS is compressed using a variational scheme [41] reducing χ. We use an upper bound forthe maximum χ = {362, 512, 724, 1024} to limit the computational resources needed. If thebound is not reached, small Schmidt values λi are discarded provided
the two criteria
P either ofsum
are met: (I) discard all λi ≤ strunc or (II) discard all λi sufficing i λ2i < (strunc)2 beginningfrom the smallest λi .
Further technical details are given as follows; The initial random TPQ-MPS state |Ψ0 〉 isprepared by applying a sequence of random two-site unitary matrices to a Néel-like productstate in the z-basis, i.e. | · · · ↑↓ · · · 〉, in a TEBD-like way. We prepare Nsamples = 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 tosumstrunc = 10−6 and strunc= 10−5 unless stated otherwise. Measurements are not independent
concerning β, since they are taken in a sequence of fixed β during each independent run ofimaginary- time evolution, while the Nsamples averages are taken from the set of independentruns. The TenPy library [42] is used for all MPS-related numerical calculations.

3

Application to the Kitaev honeycomb model

We employ TPQ-MPS to the Kitaev honeycomb (KH) model defined as [16]
X
X
X
y y
σ i σ j + Kz
σzi σzj ,H = Kxσix σ xj + K y〈i, j〉 x〈i, j〉 y(2)〈i, j〉zγ

where σi are Pauli operators σ x , σ y , and σz acting on sites i. The three sets of parallel bondson 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 σi σ j .
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 adiabaticallyconnects to the Toric Code [43]. Here, we focus on the case of K x = K y = Kz = 31 .
The KH model features a double-peak structure in the specific heat, signalling crossovers
and releasing an entropy of ∆S/N = 21 ln 2 each. The associated two energy scales are wellknown [44]: At the high-T peak, TH /K ≈ 0.5, nearest-neighbor spin-spin correlations develop
and the fractionalization into itinerant and localized Majorana fermions occurs. The latter conQγ (i)tributes to the formation of fluxes at each hexagonal plaquette given as WP = i∈P σi P ,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, wP = ±1, which are disorderedat T ≲ TH /K. Below the low-T peak, TL /K ≈ 0.016, the fluctuation of fluxes is suppressed
and we eventually find 〈WP 〉 → 1. They form the static Z2 lattice-gauge
field, fixing half ofp
the Hilbert space per unit cell. A local Hilbert space dimension of 2 per site remains whichis associated with itinerant Majorana fermions. Although the KH model at finite temperature
We note that time-dependent variational principle (TDVP) [39, 40] can be utilized as well.
The MPO representation of the imaginary time evolution is given as W II (dτ) ≡ U (dτ), following Ref. [38].
Splitting dτ = τ1 + τ2 with sufficiently chosen complex τ1 and τ2 such that U (dτ) ≈ W II (τ2 )W II (τ1 ) reduces theerror in dτ by one order.
4

5

5

SciPost Phys. 15, 206 (2023)

E

−0.24

0.0

(a)

−0.25
EGS at T = 0

−0.26
10−2

exact
iDMRG YC3-1
iDMRG YC4-1

10−1
−2
(b)10

10−1

1000.4

−0.1
−0.2

101

0.4
C

0.2

0.2
0.0
0.0

TL

10−2

0
10−1 YC8×3×2, χ10= 512

YC8×3×2, χ = 1024
YC8×4×2, χ = 362
YC8×4×2, χ = 512

0.5
0.0

0.2

0.4

TH

hWP i

0.0
1.0

1
10
YC8×4×2, χ = 724

YC8×4×2, χ = 1024
MC+FFED 12×12×2
XTRG 6×4×2

(c)

10−2

10−1

100
T = 1/β101

Figure 2: Temperature (T = β −1 )-dependent (a) energy density E and (b) specificheat 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 TL , in particular for YC8×3×2, to the finite circumference geometry used here; the ground state energies EGS of an infinitely
long Lcirc = 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 horizontalline) by about ∼ 0.01K, yielding different crossover slopes in E near TL and a shiftedpeak. 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 forthe temperature. (c) The evolution of the plaquette flux average, 〈Wp 〉. The flux-freestate at T = 0 exhibits 〈Wp 〉 = 1.
is not exactly solvable, once Z2 bond variables constituting the Z2 gauge field are treated asclassical 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 aquantum 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 ourapproach, unlike MC+FFED, is not custom tailored to the KH model and can be applied toother 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]; Theenergy density6 E rapidly decreases near TH resulting in a crossover peak in the specific heat C.
6
We are computing the energy density neglecting the left Nl = Lcirc and right Nr = Lcirc sites of the physicalsystem to obtain a better estimate of the energy density in the bulk.

6

SciPost Phys. 15, 206 (2023)

A second step of energy reduction occurs near TL . The two-step behavior is already present forsmall χ = 362 with well converged behaviour down to T ∼ 0.2 including the high-T peak inthe 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 Lcirc = 6, 8 and helical boundary
condition YC3-1 and YC4-1, respectively. The circumference seriously affect the numericallyachieved ground state energies and consequently the specific heat which can be summarized asfollows: (I) The cylinder with Lcirc = 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 Lcirc = 8 we obtain an evolution of the energy closer to MC+FFED, thus reducing thefinite-size effect signicantly. For TL ≳ T ≳ 0.2, however, TPQ-MPS overestimates E comparedto MC+FFED. Here, increasing χ gradually reduces E possibly approaching MC+FFED for
sufficiently large χ. Near T ∼ TL and below, the effect of finite χ ceases and the energy
eventually approaches both MC+FFED as well as the ground state energy. As a consequenceof 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 Z2 fluxes in Fig. 2(c) nicely marks the two peaks by an onset of nonzero
value (TH ) and the inflection point (TL ), finally approaching 〈WP 〉 → 1 at T → 0 systematicallyfor various χ.
A recent XTRG calculation applied to the KH model reports the lower-T -peak at TL ∼ 0.023
with the peak-height of ∼ 0.3 using a 6 × 4 × 2 cylinder [15]. While the circumference issimilar to our YC8×4×2, the XTRG work uses a slightly shorter cylinder, does not use helicalboundary condition, and employs a different winding scheme. The quantitative agreementof 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 TL . The comparison with our two geometries confirms that the size or shape ofthe cylinder shifts the TL 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 thedata against the bulk data at T ≤ 0.1.
In this context, we like to remark that in many frustrated spin models, the specific heatat 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.

4

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 > i0 (β j ) which
fulfills aforementioned I or II in Sec. 2) accumulated over a single imaginary-time evolution,
X X
Σtrunc (β) =λ2i (β j ) .
(3)
β j <β i>i0 (β j )7

SciPost Phys. 15, 206 (2023)

YC8×3×2, χ = 512
YC8×3×2, χ = 1024
YC8×4×2, χ = 362
YC8×4×2, χ = 512
YC8×4×2, χ = 724
YC8×4×2, χ = 1024

(a)

10−5
10−7
10−9
−2
1000 10
800
600
400
200
0
10−2

10−1

100

−1

0

101 (b)

1.00
F

χ̄

∂β Σtrunc

10−3(c) 0.99

1010
T = 1/β

10

1

10

−2

10

−1

0

10
T = 1/β101Figure 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 < β.
The second one is the product of the fidelities of the state |Ψ(β j )〉 = W II (dτ)|Ψ(β j−1 )〉 (seeRef.[36]) and |Ψ̃(β j )〉 just before and after truncation, respectively, for all truncations down
to the temperature β −1 ,Y(|〈Ψ̃(β j )|Ψ(β j )〉|2 ) ,(4)F (β) =
β j <βwhich evaluates how we deviate from the non-truncated wave function at β. The amountof truncated Schmidt values per unit of imaginary time is given as ∂β Σtrunc . In Figure 3 weshow 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. Uponlowering the temperature, ∂β Σtrunc increases and then reaches a plateau at TL ≲ T ≲ TH withvalues ∂β Σtrunc ∼ 10−5 to 10−4 depending on χ. Here, F starts to depart gradually from 1,
which is more distinct for smaller χ. At T ≲ TL 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 χ thatis reached quickly, in particular at intermediate TL ≲ T ≲ TH , taking a larger χ lowers theenergy 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 (withdimension DA) and a bigger part B, which is Schmidt decomposed as|Ψβ 〉 =DAX

λn |nA〉|nB 〉 ,

(5)

n=1

to the orthogonal basis sets {|nA〉} and {|nB 〉}. The local part A is thermalized and its density
operator is approximated by the Gibbs state in A as
ρA =

DA
X

λ2n |nA〉〈nA| ≃n=18e−β HA,ZA(6)

SciPost Phys. 15, 206 (2023)

where {|nA〉} is thought to be the energy eigenbasis of the subsystem’s Hamiltonian
HA. ForpA−β EnA/2
its eigenvalues {En }, the Schmidt coefficient λn is represented as e
/ ZA, and we findDA −β E A/2Xne|Ψβ 〉 ≃|nA〉|nB 〉 .
p
ZA
n=1

(7)

Note here that {|nB 〉} is left unknown. We finally truncate |Ψβ 〉 as DA → χ in Eq.(7), discarding
the basis states with small weight. Specifically, information of |nA〉 belonging to higher EnA
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 theinformation 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 thevariance increases at low temperatures by the order of ≤ e−SN /T , and when the entropy SNof size N is sufficiently large, the increase is moderately suppressed. This fact supports theapplication of TPQ methods to highly frustrated quantum magnets including the present KH
model and Kagome or related lattice models [46–51]. However, the idea of relying on thelarge entropy does not apply to TPQ-MPS: in the first proposal of TPQ-MPS in Ref. [19] someof 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.

5

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 Z2 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 rangeof temperatures. This allows the state starting from random at high temperature (initial state)
to gradually reach the qualitatively different long-range entangled topological ordered groundstate. The application to non-frustrated model is expected to be promising because, unlike forthe original TPQ, a high entropy density is not required at low temperature to attain a reasonable accuracy at a moderate numerical cost.

9

SciPost Phys. 15, 206 (2023)

Acknowledgments
We thank J. Nasu for providing us with reference data. We acknowledge the use of computational resources of the supercomputer Fugaku provided by the RIKEN AICS through the HPCI
System Research Project (Project ID: hp210321) and of the Scientific Computing section of
the Research Support Division at the Okinawa Institute of Science and Technology Graduate
University (OIST). M.G. acknowledges support by the Theory of Quantum Matter Unit at OIST.
Funding information This work was supported by a Grant-in-Aid for Transformative Research Areas “The Natural Laws of Extreme Universe— A New Paradigm for Spacetimeand Matter from Quantum Information” (No. 21H05191) and JSPS KAKENHI (Grants No.
JP21K03440 and JP22K14008). A.I. was supported by JSPS Research Fellowship (Grant No.
21J21992).

A

Random sampling average

TPQ-MPS is a random sampling method using the MPS representation of the quantum manybody wave function. Since the quality of the MPS state relies on its bond dimension χ practically accessible in the computation, a smaller χ would require a larger number of independentruns to be averaged over. This number is generally by orders of magnitude smaller than with
METTS when applying them to the same system.
Let us first highlight the difference from METTS, which is a moderately mixed quantum
thermal state. METTS starts from a (classical) product state with χ = 1, and χ grows upon
imaginary time evolution. Accordingly, the entanglement stored is moderate and does not suffice the amount of entropy of the thermal state: This is compensated by taking an average oversamples and constructing a mixed state. In METTS, a Markov chain scheme is employed to increase the efficiency of the sampling, while the sample number typically amounts to order-102 .
Their χ can be kept smaller by increasing the sample average. TPQ-MPS does not rely muchon the sample average by storing the entanglement as maximally as possible in a single pure
TPQ-MPS form, which usually requires large χ. However, by relying on the auxiliary degrees
of freedom on the edges, the entanglement entropy distribution becomes nearly flat (without
auxiliaries, χ at the edge becomes very small), and this allows us to store the entanglement
maximally for a given χ. There exists a trade-off between the size of χ and of the sample
average, but each method has a specific control parameter to keep them in the practicallyreasonable range. A more detailed comparison with METTS is given in Ref. [2].
One may anticipate that the present two-dimensional KH model might require more samples than a 1D system [19]. As illustrated in Fig. 1 in the main text, the higher the purity is,
fewer samples Nsamp are needed to safely reproduce the thermal quantum state. In TPQ-MPS,
the Nsample -independent runs are performed starting from the independent initial random MPS,
Nsampleyielding a set of unnormalized states over different β j for each, {|Ψ (l) (β j )〉}l=1
average of physical quantities O is taken as. The random

PNsample〈Ψ (l) (β j )|O|Ψ (l) (β j )〉l=1
〈O〉 = PN.
sample(l) (β )|Ψ (l) (β )〉
〈Ψ
j
j
l=1

(A.1)

Here, the summations over samples are taken independently betweenP
the numerator and deM
nominator, since the partition function is given by the denominator Z = l=1 〈Ψ (l) (β j )|Ψ (l) (β j )〉(the reasoning for why the average taken by the normalized |Ψ (l) (β j )〉 does not provide thecorrect sampling average is analytically shown in Ref. [2]).
10SciPost Phys. 15, 206 (2023)

3

(a)

10−1

10−4
10−5
10−6

2

YC8×3×2, χ = 512
YC8×3×2, χ = 1024
YC8×4×2, χ = 362
YC8×4×2, χ = 512
YC8×4×2, χ = 724
YC8×4×2, χ = 1024

10−3

10

−2

10

S̄E

δz 2

10−2

1

(b)

−1

10

0

10

1

0

0

T = 1/β

5

10

15
20
iterations

25

30

Figure 4: (a) Normalized fluctuation of 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 andodd bonds similar to the TEBD algorithm. S¯E 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.
The number Nsample required can be measured using a quantity called normalized fluctuation of partition function (NFPF),
Var(〈Ψ(β j )|Ψ(β j )〉)δz 2 = 2 ,〈Ψ(β j )|Ψ(β j )〉(A.2)

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, Nsample , 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 aset of data given in Fig. 2. The largest δz 2 at low-T ranges at 10−1 − 101 , which is comparable
to the value for the 1D Heisenberg model [2] using Nsample = 100. Based on this comparison,
we also adopt Nsample = 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 aremuch 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 ≲ 100 , 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 thehigh 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 applicationof random two-site unitary matrices. Although not relevant to our model, an advantage ofthis 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 asaturated entanglement entropy. Figure 4(b) shows the evolution of entanglement entropy ofbipartition S¯E , averaged over the physical sites, as a function of iterations. While the bonddimension of the state generally doubles after each iterations, S¯E is not fully saturated after
11

SciPost Phys. 15, 206 (2023)

N = 5 iterations, which is required to reach χRMPS = 32. Instead, saturation is reached ataround N ∼ 10 iterations. Our choice of N = 25 iterations is well within the saturated regime.

References[1] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information, Cambridge University Press, Cambridge, UK, ISBN 9781107002173 (2000),
doi:10.1017/CBO9780511976667.
[2] A. Iwaki and C. Hotta, Purity of thermal mixed quantum states, Phys. Rev. B 106, 094409(2022), doi:10.1103/PhysRevB.106.094409.
[3] S. Popescu, A. J. Short and A. Winter, Entanglement and the foundations of statisticalmechanics, Nat. Phys. 2, 754 (2006), doi:10.1038/nphys444.
[4] S. Goldstein, J. L. Lebowitz, R. Tumulka and N. Zanghì, Canonical typicality, Phys. Rev.
Lett. 96, 050403 (2006), doi:10.1103/PhysRevLett.96.050403.
[5] F. Verstraete, J. J. García-Ripoll and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004),
doi:10.1103/PhysRevLett.93.207204.
[6] M. Zwolak and G. Vidal, Mixed-state dynamics in one-dimensional quantum lattice systems:
A time-dependent superoperator renormalization algorithm, Phys. Rev. Lett. 93, 207205(2004), doi:10.1103/PhysRevLett.93.207205.
[7] A. E. Feiguin and S. R. White, Finite-temperature density matrix renormalization using an enlarged Hilbert space, Phys. Rev. B 72, 220401 (2005),
doi:10.1103/PhysRevB.72.220401.
[8] G. De las Cuevas, N. Schuch, D. Pérez-García and J. I. Cirac, Purifications of multipartite states: Limitations and constructive methods, New J. Phys. 15, 123021 (2013),
doi:10.1088/1367-2630/15/12/123021.
[9] S. Goto, R. Kaneko and I. Danshita, Matrix product state approach for a quantum systemat finite temperatures using random phases and Trotter gates, Phys. Rev. B 104, 045133(2021), doi:10.1103/PhysRevB.104.045133.
[10] M. M. Wolf, F. Verstraete, M. B. Hastings and J. I. Cirac, Area laws in quantum
systems: Mutual information and correlations, Phys. Rev. Lett. 100, 070502 (2008),
doi:10.1103/PhysRevLett.100.070502.
[11] T. Barthel, One-dimensional quantum systems at finite temperatures can be simulated efficiently on classical computers, (arXiv preprint) doi:10.48550/arXiv.1708.09349.
[12] T. Kuwahara, Á. M. Alhambra and A. Anshu, Improved thermal area law and quasilinear time algorithm for quantum Gibbs states, Phys. Rev. X 11, 011047 (2021),
doi:10.1103/PhysRevX.11.011047.
[13] B.-B. Chen, L. Chen, Z. Chen, W. Li and A. Weichselbaum, Exponential thermal tensor network approach for quantum lattice models, Phys. Rev. X 8, 031082 (2018),
doi:10.1103/PhysRevX.8.031082.

12

SciPost Phys. 15, 206 (2023)[14] H. Li, B.-B. Chen, Z. Chen, J. von Delft, A. Weichselbaum and W. Li, Thermal tensor
renormalization group simulations of square-lattice quantum spin models, Phys. Rev. B
100, 045110 (2019), doi:10.1103/PhysRevB.100.045110.
[15] H. Li, D.-W. Qu, H.-K. Zhang, Y.-Z. Jia, S.-S. Gong, Y. Qi and W. Li, Universal
thermodynamics in the Kitaev fractional liquid, Phys. Rev. Res. 2, 043015 (2020),
doi:10.1103/PhysRevResearch.2.043015.
[16] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006),
doi:10.1016/j.aop.2005.10.005.
[17] J. R. Garrison and T. Grover, Does a single eigenstate encode the full Hamiltonian?, Phys.
Rev. X 8, 021026 (2018), doi:10.1103/PhysRevX.8.021026.
[18] Y. O. Nakagawa, M. Watanabe, H. Fujita and S. Sugiura, Universality in volumelaw entanglement of scrambled pure quantum states, Nat. Commun. 9, 1635 (2018),
doi:10.1038/s41467-018-03883-9.
[19] A. Iwaki, A. Shimizu and C. Hotta, Thermal pure quantum matrix productstates recovering a volume law entanglement, Phys. Rev. Res. 3, L022015 (2021),
doi:10.1103/PhysRevResearch.3.L022015.
[20] M. Imada and M. Takahashi, Quantum transfer Monte Carlo method for finite temperature
properties and quantum molecular dynamics method for dynamical correlation functions,J. Phys. Soc. Jpn. 55, 3354 (1986), doi:10.1143/JPSJ.55.3354.
[21] J. Jaklič and P. Prelovšek, Lanczos method for the calculation of finite-temperature quantities in correlated systems, Phys. Rev. B 49, 5065 (1994), doi:10.1103/PhysRevB.49.5065.
[22] A. Hams and H. De Raedt, Fast algorithm for finding the eigenvalue distribution of verylarge matrices, Phys. Rev. E 62, 4365 (2000), doi:10.1103/PhysRevE.62.4365.
[23] T. Iitaka and T. Ebisuzaki, Random phase vector for calculating the trace of a large matrix,
Phys. Rev. E 69, 057701 (2004), doi:10.1103/PhysRevE.69.057701.
[24] S. Sugiura and A. Shimizu, Thermal pure quantum states at finite temperature, Phys. Rev.
Lett. 108, 240401 (2012), doi:10.1103/PhysRevLett.108.240401.
[25] S. Sugiura and A. Shimizu, Canonical thermal pure quantum state, Phys. Rev. Lett. 111,
010401 (2013), doi:10.1103/PhysRevLett.111.010401.
[26] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand and T. Xiang,
Gapless spin-liquid ground state in the s = 1/2 Kagome antiferromagnet, Phys. Rev. Lett.
118, 137202 (2017), doi:10.1103/PhysRevLett.118.137202.
[27] Y.-C. He, M. P. Zaletel, M. Oshikawa and F. Pollmann, Signatures of Dirac cones ina DMRG study of the Kagome Heisenberg model, Phys. Rev. X 7, 031020 (2017),
doi:10.1103/PhysRevX.7.031020.
[28] M. Gohlke, R. Verresen, R. Moessner and F. Pollmann, Dynamics of the Kitaev-Heisenberg
model, Phys. Rev. Lett. 119, 157203 (2017), doi:10.1103/PhysRevLett.119.157203.
[29] S. R. White, Minimally entangled typical quantum states at finite temperature, Phys. Rev.
Lett. 102, 190601 (2009), doi:10.1103/PhysRevLett.102.190601.
[30] E. M. Stoudenmire and S. R. White, Minimally entangled typical thermal state algorithms,
New J. Phys. 12, 055026 (2010), doi:10.1088/1367-2630/12/5/055026.
13

SciPost Phys. 15, 206 (2023)[31] J. L. Jimenez et al. A quantum magnetic analogue to the critical point of water, Nature
592, 370 (2021), doi:10.1038/s41586-021-03411-8.
[32] A. Wietek, P. Corboz, S. Wessel, B. Normand, F. Mila and A. Honecker, Thermodynamic
properties of the Shastry-Sutherland model throughout the dimer-product phase, Phys. Rev.
Res. 1, 033038 (2019), doi:10.1103/PhysRevResearch.1.033038.
[33] M. Gohlke, R. Moessner and F. Pollmann, Dynamical and topological properties of
the Kitaev model in a [111] magnetic field, Phys. Rev. B 98, 014418 (2018),
doi:10.1103/PhysRevB.98.014418.
[34] S. Yan, D. A. Huse and S. R. White, Spin-liquid ground state of the s = 1/2 Kagome Heisenberg antiferromagnet, Science 332, 1173 (2011), doi:10.1126/science.1201080.
[35] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev.
Lett. 93, 040502 (2004), doi:10.1103/PhysRevLett.93.040502.
[36] S. R. White and A. E. Feiguin, Real-time evolution using the density matrix renormalizationgroup, Phys. Rev. Lett. 93, 076401 (2004), doi:10.1103/PhysRevLett.93.076401.
[37] A. J. Daley, C. Kollath, U. Schollwöck and G. Vidal, Time-dependent density-matrix
renormalization-group using adaptive effective Hilbert spaces, J. Stat. Mech.: Theor. Exp.
P04005 (2004), doi:10.1088/1742-5468/2004/04/p04005.
[38] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore and F. Pollmann, Time-evolving
a matrix product state with long-ranged interactions, Phys. Rev. B 91, 165112 (2015),
doi:10.1103/PhysRevB.91.165112.
[39] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde and F. Verstraete, Timedependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011),
doi:10.1103/PhysRevLett.107.070601.
[40] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken and F. Verstraete, Unifying time
evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016),
doi:10.1103/PhysRevB.94.165116.
[41] U. Schollwöck, The density-matrix renormalization group in the age of matrix productstates, Ann. Phys. 326, 96 (2011), doi:10.1016/j.aop.2010.09.012.
[42] J. Hauschild and F. Pollmann, Efficient numerical simulations with tensor networks: Tensor network Python (TeNPy), SciPost Phys. Lect. Notes 5 (2018),
doi:10.21468/SciPostPhysLectNotes.5.
[43] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. 303, 2 (2003),
doi:10.1016/S0003-4916(02)00018-0.
[44] J. Nasu, M. Udagawa and Y. Motome, Thermal fractionalization of quantum spins in a Kitaev model: Temperature-linear specific heat and coherent transport of Majorana fermions,Phys. Rev. B 92, 115122 (2015), doi:10.1103/PhysRevB.92.115122.
[45] T. Shimokawa and H. Kawamura, Finite-temperature crossover phenomenon in the s = 1/2
antiferromagnetic Heisenberg model on the Kagome lattice, J. Phys. Soc. Jpn. 85, 113702(2016), doi:10.7566/JPSJ.85.113702.
[46] Y. Yamaji, T. Suzuki, T. Yamada, S.-i. Suga, N. Kawashima and M. Imada, Clues and criteria
for designing a Kitaev spin liquid revealed by thermal and spin excitations of the honeycomb
iridate Na2 IrO3 , Phys. Rev. B 93, 174425 (2016), doi:10.1103/PhysRevB.93.174425.
14

SciPost Phys. 15, 206 (2023)[47] H. Endo, C. Hotta and A. Shimizu, From linear to nonlinear responses of thermal pure quantum states, Phys. Rev. Lett. 121, 220601 (2018), doi:10.1103/PhysRevLett.121.220601.
[48] C. Hotta and K. Asano, Magnetic susceptibility of quantum spin systems calculated by sine
square deformation: One-dimensional, square lattice, and Kagome lattice Heisenberg antiferromagnets, Phys. Rev. B 98, 140405 (2018), doi:10.1103/PhysRevB.98.140405.
[49] I. Rousochatzakis, S. Kourtis, J. Knolle, R. Moessner and N. B. Perkins, Quantum spin
liquid at finite temperature: Proximate dynamics and persistent typicality, Phys. Rev. B
100, 045117 (2019), doi:10.1103/PhysRevB.100.045117.
[50] J. Richter, O. Derzhko and J. Schnack, Thermodynamics of the spin-halfsquare Kagome lattice antiferromagnet, Phys. Rev. B 105, 144427 (2022),
doi:10.1103/PhysRevB.105.144427.
[51] J. Richter and J. Schnack, Magnetism of the s = 12 J1 − J2 square-Kagome lattice antiferromagnet, Phys. Rev. B 107, 245115 (2023), doi:10.1103/PhysRevB.107.245115.

15