Thermal Pure Quantum Matrix Product States Recovering a Volume Law Entanglement

We propose a way to construct thermal pure quantum matrix product state (TPQ-MPS) that can simulate finite temperature quantum many body systems with a minimal numerical cost comparable to the matrix product algorithm for the ground state. The MPS was originally designed for the wave function with area-law entanglement. However, by attaching the auxiliary sites to the edges of the random matrix product state, we find that the degree of entanglement is automatically tuned so as to recover the volume law of the entanglement entropy that characterizes the TPQ state. The finite temperature physical quantities of the transverse Ising and the spin-1/2 Heisenberg chains evaluated by a TPQ-MPS show excellent agreement even for bond dimension ∼ 10-20 with those of the exact results.

Introduction. Finding a good description of typical wavefunctions of quantum many body states at finite temperature has been a challenge in condensed matter theory. Traditional quantum mechanical representation of equilibrium states rely on the density operators of Gibbs' ensembles, which are the classical mixtures of an exponentially large numbers of pure quantum states.
This construction does not allow us to simulate sufficiently large quantum systems of physical interest, since the available numerical devices, e.g. stochastic quantum Monte Carlo (QMC) methods or some size-free methods [1][2][3][4] are limited.
The key conceptual development against the ensemble physics is the typicality [5][6][7][8][9][10]; there exists a single thermal pure quantum (TPQ) state that solely represents the thermal equilibrium [11][12][13][14][15][16][17]. It then happens that for the description of any of the equilibrium quantum states, one can choose an arbitrary degrees of classical mixture, from the purity-1/(e Θ(N ) ) ensemble to purity-1 TPQ state. Let us consider a size-N system with a Hamiltonian H, and its subsystem-A of size-n. Since the equilbrium state at inverse temperature β can be described equivalently by the TPQ state |ψ β and by the Gibbs state, the local observables O A also fulfill β ψ|O A |ψ β = tr(O A e −βH )/tr(e −βH ). ( This directly indicates the equality of the local density matrix and the local canonical ensemble, ρ n (|ψ β ) = trĀ(e −βH )/tr(e −βH ). Accordingly, the von Neumann entropy S n = −Tr(ρ n ln ρ n ) obtained by the reduced density operator ρ n of the subsystem-A becomes the entanglement entropy(EE) of a TPQ state, and is related to thermodynamic entropy density s th as [18] S n /n = s th , (1 ≪ n ≪ N ).
The EE of the TPQ state thus needs to fulfill the volume law, and indeed, in the similar context, the Page curve of the second Renyi entropy in a finite open boundary system is observed in thegexacth TPQ-state [19].
Practically, however, constructing suchgexact" TPQ state, which we call full-TPQ state [11,15,16], requires a cost only slightly smaller than the conventional finite temperature diagonalization methods [20][21][22], and is available only up to ∼ 2N of that of the Gibbs state. As for the ground state, the approximate forms of the pure wave functions are established by the density matrix renormalization(DMRG) or matrix product states (MPS) [23][24][25]. However, their application to the TPQ state has never been tested because the MPS description has an area law entanglement by construction [26,27], and is apparently unsuitable for the finite temperature case where the entanglement blows up massively.
In this Letter, we show that the TPQ-MPS state is realized by attaching appropriate auxiliary sites at both edges of the system, which work as an entanglement-bath and make the system highly entangled. By successively operating the Hamiltonian to the random matrix product state (RMPS) with the auxiliaries, the MPS is annealed down to lower temperature where we find that the system recovers the volume law entanglement when measured from the system center toward the very edges. We demonstrate that the TPQ-MPS wave functions give accurate evaluation of the physical quantities in typical quantum spin models without taking ensemble average. The computational cost is significantly reduced to that of the MPS of the ground state, and the accessible system size increases to N 100.
So far, the MPS methods at finite temperature like minimally entangled typical thermal states (METTS) [28,29] and matrix product purification(MPP) [30][31][32], both with intermediate degrees of purity, have relied on some sort of "ensemble" averages to compensate for the low entanglement properties of MPS. The METTS starts from the minimally entangled product state which requires a large number of sampling average of 100. The MPP adopts extra N -ancillary system which are traced out to obtain the mixed state of the system [33]; This expression is mathematically redundant since they replace a mixed state in an irreducible representation in the minimum Hilbert space of dimension d N by a reducible representation in a space of dimension d 2N , and does not save our numerical cost by itself. In our TPQ-MPS construction, the redundancy is limited to d N +χ /d N = d χ , which does not grow with increasing N [34].
Random initial state.
We consider a onedimensional(1D) lattice system consisting of N sites with open boundary condition (OBC) where each site hosts d-dimensional degrees of freedom, and two auxiliaries attached at both edges having χ-dimensions. The RMPS of such a system shown in Fig.1(a) is given as [35] where the A [m]im is the dχ 2 matrix on the m-th site with i m = 1, · · · , d, and is explicitly given using the dχ × dχ random unitary matrix U as A [m]im αβ = U (i,α)(1,β) or U (1,α)(i,β) which fulfill the left or right canonical form, respectively. Here |aux L/R α L/R is the right/left auxiliary state with α L/R = 1-χ. This RMPS reproduces the physical quantities at T = ∞ with the variance of order χ −2 , which can be shown analytically as follows; taking A on the l.h.s/r.h.s. of a one-site operatorÔ i as left/right canonical form, we have ψ RM |ψ RM = χ. By taking account of the formula for the random average,

and its variance as
Typicality of the RMPS is studied and confirmed numerically in the similar context [35,36], followed by several proposals to stochastically construct microcanonical and canonical ensembles of RMPS [37][38][39]. In our work, the RMPS is constructed not by using Eq.(3) but by preparing a tensor Γ [m]im αmβm (∈ C) of bond dimension χ, whose elements follow the Gaussian distribution (see Fig.1(a)). It can be straightforwardly shown that after transforming Eq.(3) into the canonical form by the successive Schmidt-decomposition [40,41], the obtained matrices also form an equivalent RMPS.
whereĥ is the Hamiltonian divided by N , and l is a parameter larger than the maximum eigenenergy, which is necessary to generate sharp microcanonical energy distribution(see Ref. [15].) Here, |k is the TPQ state at a temperature k B T (k) = N (l − u k )/2k [42] with energy u k = k|ĥ|k / k|k . In our algorithm, (l −ĥ) is represented by a matrix product operator (MPO) of bonddimension D that depends onĥ, and applying this MPO at each step to |k multiplies the matrix dimension to Dχ. Before truncating the dimension of the enlarged matrix down to χ, we transform the MPS to its canonical form including auxiliary sites in order to minimizes the truncation error [43]. The process is repeated until k B T (k) reaches low enough temperature β −1 max , and the effective bond dimension, χ eff i (i = 0, N ), which is the number of finite eigenvalues of the Schmidt decomposition λ i on the i-th bond, can change automatically within 1 ≤ χ eff i ≤ χ. We consider an operatorÂ that can be described by a low-order polynomial of local observables. At each step, suchÂ is evaluated and stored, Â = k|Â|k / k|k is the physical quantities at k B T (k) . Instead of directly adopting this form, one can generate the physical quantities for arbitrary temperatures β k|Â|k + 1 . The quality of the TPQ-MPS is tested by the comparison of A β,N with the exact or nearly exact solution by the counterparts like the exact diagonalization (ED), QMC or quantum transfer matrix (QTM) [44,45] methods [46].
The required number of random averages of mTPQ runs N ran is as small as in the full-TPQ that uses the full Hilbert space, e.g. typically less than 5 for N 32. This can be seen from a benchmark results in Fig. 1(b); When performing 10-mTPQ runs with and without auxiliaries, the variance of the former turned out to be small by orders of magnitude, particularly at higher temperature. The variance of TPQ-MPS degreases at lower T , contrarily to the full-TPQ. The higher performance despite its being an approximate method is possibly because the basis in the MPS construction is optimally biased at low temperature to those favoring the low energy states. The difference from the original TPQ also lies that the range of mTPQ temperature k B T (k) depends much on l. Usually, starting from smaller l will accelerate the convergence and we reach the same temperature with smaller k max . However, because of finite χ, the starting temperature at k ∼ 1 for l ≪ Θ(h) is kept to k B T Θ(h) (h being the typical local energy scale of the Hamiltonian) in which case the canonical summation becomes inaccurate. Whereas, l Θ(10h) sacrifices the low temperature information, and one needs to set proper l depending on the models. Volume law and the auxiliaries. Let us first examine the basic entanglement properties of the TPQ-MPS.
Here, we benchmark the 1D transverse Ising chain, i , with σ z = ±1, whose MPO has D = 3, and take J = g = 1. Figures 1(c) and 1(d) show the EE, S i = −Tr(ρ i lnρ i ) whereρ i is the density matrix when dividing the system into left and right parts at the i-th bond. It is known that S i of the TPQ state follows a Page curve [47], which increases linearly from the edges that reflect the volume law and saturates toward the maximum at the center. The standard MPS without auxiliaries in Fig.1(d) indeed follows a page-like curve at the edge up to the bond dimension where S i becomes flat because of the upper bound by the bond dimension χ. Here, the lower temperature/larger k requires the smaller χ. Whereas in TPQ-MPS, the auxiliaries introduces to the edge a large S 0 that depends only on the effective bond dimension at the edge, generating a nearly flat S i throughout the system. However, if we divide the system by setting the subsystem at the center, namely cutting the two bonds at equal distances from the center, we find physically meaningful entanglement properties. Figure 1(e) is the corresponding EE, S c n = −Tr(ρ n lnρ n ), as a function of the size n of the subsystem. One finds that S c n increases linearly in n until it saturates to the upper bound. At k 600 where the system reaches the temperature k B T (k) 0.33, S c n continues to increase up to the very edges of the system.

Entropy.
These results indicate that the TPQ-MPS wave function has acquired qualitatively different degrees of freedom in its state space, just by attaching the two auxiliaries. While the upper bound of the EE, 2 ln χ, is only twice as large as the case without auxiliaries, the EE continues to go up until that bound since the MPS does not feel the edge. This is in sharp contrast to the standard MPS where the entanglement is suppressed toward both edges( Fig. 1(e), bold dots). The present scheme allows the evaluation of the thermodynamic entropy density s th using Eq.(2); We perform a linear fit of a series of S c n for N = 64 and for various χ, and compare its slope with s th obtained by integrating the specific heat of the QMC calculation. As shown in Fig.1(f), the slope of S c n agrees with s th , asymptotically approaching the QMC data with increasing χ. Benchmarks. Figure 2(a) shows the energy density e = E/N of the transverse Ising model for N = 16 which is in good agreement with the QMC results with open boundary condition (OBC), where we put together the periodic boundary (PBC) ones of the same size. The N = 16 and 64 cases are also compared in the inset, demonstrating that the finite size effect is much larger than the difference between the TPQ-MPS and QMC results. Already at χ 10, the error is converged as we see in Fig.2(b) which is smaller than the variance of the QMC results over 100 independent runs each with 200000 averages, regardless of size-N and χ. The specific heat C/N also gives excellent agreement with the QMC results of the same size (see Fig.2(c)). Figure 2(d) gives the spatial distribution of the bond-and site-energies (J-and g-terms of the Hamiltonian) for several temperatures at N = 64. They perfectly follow those of OBC obtained by QMC for k B T = 0.5 even at the very edges that show downturn/upturn. The variance of mTPQ runs are kept small enough; 10-20 averages for N = 16 and 5 for N = 64, 96 in the figures.
Next, we test our method with the Heisenberg chain described byĤ = N −1 i=1 Jŝ iŝi+1 , with s z i = ±1/2. The MPO for this Hamiltonian has D = 5 that would increase the truncation error. In fact, it is known for the MPS ground state that the model requires much larger χ compared to the product-type ground state of the transverse Ising model. Figure 3(a) shows E/N and C/N for several system sizes in comparison with the N = 16 ED with OBC and N = ∞ exact solution (QTM) [45], which shows good agreement for the same order of χ as the transverse Ising model. We also plot in Fig. 3(b) the susceptibility χ s to see how much a well-known logarithmic singularity at the lowest temperature can be traced by the TPQ-MPS. The drop of χ s at N = 16 is almost perfectly reproduced already for χ = 20. Also, the larger size results are in reasonable agreement with QTM.
Conclusion. We realized the TPQ-MPS by attaching the edge-auxiliaries of dimension χ to the MPS, showing that the EE of the subsystem cut out from the center continues to show a volume law up to the very edges of the system, particularly at low temperature, by setting the bond dimension χ to the realistically small values. Physical properties are evaluated almost free of random sampling. Such construction enables the application to higher dimensions, possibly less costly than the aforedeveloped tensor network approaches [48][49][50][51].