Exponential Thermal Tensor Network Approach for Quantum Lattice Models

We speed up thermal simulations of quantum many-body systems in both one- (1D) and two-dimensional (2D) models in an exponential way by iteratively projecting the thermal density matrix $\widehat{\rho }=e^{-\beta \widehat{H}}$ onto itself. We refer to this scheme of doubling $\beta$ in each step of the imaginary time evolution as the exponential tensor renormalization group (XTRG). This approach is in stark contrast to conventional Trotter-Suzuki-type methods which evolve $\widehat{\rho }$ on a linear quasicontinuous grid in inverse temperature $\beta \equiv 1/T$. As an aside, the large steps in XTRG allow one to swiftly jump across finite-temperature phase transitions, i.e., without the need to resolve each singularly expensive phase-transition point right away, e.g., when interested in low-energy behavior. A fine temperature resolution can be obtained, nevertheless, by using interleaved temperature grids. In general, XTRG can reach low temperatures exponentially fast and, thus, not only saves computational time but also merits better accuracy due to significantly fewer truncation steps. For similar reasons, we also find that the series expansion thermal tensor network approach benefits in both efficiency and precision, from the logarithmic temperature scale setup. We work in an (effective) 1D setting exploiting matrix product operators (MPOs), which allows us to fully and uniquely implement non-Abelian and Abelian symmetries to greatly enhance numerical performance. We use our XTRG machinery to explore the thermal properties of Heisenberg models on 1D chains and 2D square and triangular lattices down to low temperatures approaching ground-state properties. The entanglement properties, as well as the renormalization-group flow of entanglement spectra in MPOs, are discussed, where logarithmic entropies (approximately $\ln \beta$) are shown in both spin chains and square-lattice models with gapless towers of states. We also reveal that XTRG can be employed to accurately simulate the Heisenberg $XXZ$ model on the square lattice which undergoes a thermal phase transition. We determine its critical temperature based on thermal physical observables, as well as entanglement measures. Overall, we demonstrate that XTRG provides an elegant, versatile, and highly competitive approach to explore thermal properties, including finite-temperature thermal phase transitions as well as the different ordering tendencies at various temperature scales for frustrated systems.

Popular Summary

Strong interactions among atomic particles can lead to materials, such as some high-temperature superconductors, with unusual electronic and magnetic properties. To explore these materials and their potential applications, researchers rely on computer simulations. Reliable numerical access to thermodynamic properties such as magnetization and heat capacity allows researchers to make direct comparisons to experiments. However, such simulations are exceptionally difficult in the presence of strongly correlated quantum many-body behavior. We advocate a simple attractive way to exponentially speed up thermal simulations of quantum many-body systems in both one- and two-dimensional materials.

The thermal equilibrium of a quantum system is fully described by a thermal density matrix, from which all thermodynamic properties can be derived. Starting from close to infinite temperature, our technique repeatedly multiplies the thermal density matrix onto itself, which lowers the temperature at each iteration by a factor of 2. This approach allows us to swiftly jump across finite-temperature phase transitions. In general, we can reach low temperatures exponentially fast, and thus not only save computational time but also achieve better accuracy because of significantly fewer truncation steps.

Our approach is a simple yet accurate approach to thermal simulations for one- and two-dimensional quantum lattice models that removes several obstacles inherent to other methods. We believe that this opens up promising applications in exploring experimentally relevant problems in correlated quantum matter at low temperatures.

Figure 1

(a) Finite-size spectra $E_s$ with low-energy level spacing $\delta E\sim 1/L$. By requiring $T>\delta E$ for thermal averaging, this suggests $L\gtrsim \beta$, and therefore a thermal correlation length $\xi \sim \beta$. (b) A large or infinite system has an effective thermal cutoff $\xi \sim \beta$ in system length when measuring local observables. Therefore, provided $\xi \lesssim L$, finite systems can be used, as a very good approximation, to simulate thermal properties in the thermodynamic limit.

Figure 2

(a) Linear versus logarithmic temperature scale employed in thermal simulations. (b) A single step in XTRG evolution by projecting MPO ${\rho }_n$ (at $\beta =2^n\tau$) to itself. Following common notation, tensor networks are graphically depicted by blocks (i.e., tensors) connected by lines which are to be contracted. Here, vertical lines indicate physical state spaces, whereas horizontal lines indicate virtual or bond state spaces. The exploitation of symmetries, quite generally, mandates directed lines; hence, each line carries an arrow.

Figure 3

(a) SETTN initialization of ${\rho }_0$ (with $\beta ={\tau }_0$) using Maclaurin expansion with coefficients ${\omega }_k\equiv (-{\tau }_0{)}^k/k!$ [see Eq. (7a)]. (b) Pointwise SETTN algorithm exploiting the logarithmic $\beta$ scale, here with coefficients ${\omega }_k\equiv ({\tau }_n-{\tau }_{n+1}{)}^k/k!$ [see Eq. (7b)].

Figure 4

(a) Relative errors of free energy of an $L=18$ spin-$1/2$ Heisenberg chain (PBC), calculated by TTN algorithms including LTRG, XTRG, and SETTN relative to ED, with SU(2) symmetry implemented, throughout. Data for different methods share the same color for the same $D^{*}$. Here, $\rho ({\tau }_0=0.01)$ is initialized using SETTN. (b) Bipartite entanglement entropy $S_E$ in the middle of the MPO, which increases monotonically. The time costs of the LTRG calculations are 8 (10) times as long as SETTN and 180 (113) times those of XTRG, in $D^{*}=100$ (200) calculations, respectively.

Figure 5

(a) Relative errors of free energy in an $L=50$ spin-$1/2$ $XY$ chain with OBCs, calculated by TTN algorithms including LTRG, XTRG, and SETTN, with U(1) symmetry encoded, relative to the analytical solution. Similar representation as in Fig. 4, otherwise. Here, $\rho ({\tau }_0=0.01)$ is again also initialized using SETTN. (b) Entanglement entropy in the middle of the system. Overall, time costs of LTRG are 7 (7.4) times that of the SETTN run and 363 (183) times that of XTRG, for bond dimension $D=100$ (200), respectively.

Figure 6

Specific heat of a Heisenberg chain of length $L=300$, with up to $D^{*}=250$ multiplets retained. For low temperatures $T\ll 1$, a universal linear behavior versus T is observed, i.e., $c_V=\alpha T^{\eta }$ with fitted exponent $\eta =0.996$ and slope $\alpha \simeq 2/3$, in the regime $T\le 0.025$. Exploiting the fact that $\alpha =[(\pi c)/(3v)]$ ($v=\pi /2$ for a spin-$1/2$ Heisenberg chain), we extract the central charge $c\simeq 1$. For large $T$, the specific heat shows a universal $1/T^2$ temperature dependence (the fit shown is performed for $T>15$). The inset shows the same data on a linear vertical scale.

Figure 7

Free energy $f$, internal energy $u$, and specific heat $c_V$ of the isotropic Heisenberg square lattice on an $L=4$, $W=4$ [upper panels (a)–(c)] and $L=16$, $W=5$ [lower panels (d)–(f)] square lattice. The XTRG results are in very good agreement with those of ED for a $4\times 4$ lattice and quantum Wang-Landau (QWL) (expansion order 1200, sweep number 30, down to $T/J=0.1$) for the $16\times 5$ lattice. With $D^{*}$ the number of retained bond multiplets in the SU(2) MPO and $D$ the corresponding number of individual U(1) states, data is shown for $D^{*}$ ($D$) equal to $D^{*}=400$ (1462), $D^{*}=500$ (1832), and $D^{*}=600$ (2400), resulting in a maximal ratio $D/D^{*}\simeq 4$. The maximal truncation errors, i.e., at largest $\beta$, are $\delta \rho \sim 1.4\times {10}^{-5}$ for a $4\times 4$ lattice at $D^{*}=400$ and $\delta \rho \sim 1.9\times {10}^{-4}$ for a $16\times 5$ lattice at $D^{*}=600$. These truncation errors directly scale with the relative error in the partition function and, thus, the free energy. Extrapolating $1/D^{*}\to 0$ by a quadratic polynomial based on $D^{*}\sim 300–600$ data results in perfect agreement (green dashed line) with QWL data in (d). The low-temperature regions in (c),(f) are shaded where we have a limited accuracy of $c_V$. The insets in (c),(f) fit the high-$T$ specific heat with $T^{\mu }$, resulting in $\mu =-2.01$ in both cases (based on the last eight points, i.e., $T>30$).

Figure 8

Thermodynamics of the triangular lattice Heisenberg model defined on two $W=4$ geometries, OS with OBCs on both directions and $Y$ cylinder (YC) with various lengths $L$, i.e., PBCs along the vertical direction [see the inset in (d) for our specific choice of the YC geometry, as well as the TLH lattice layout in Fig. 9 below]. The presented quantities include (a) specific heat $c_V$, (b) static structure factor $S(q)$ at $q=K$ and $q=M$ at the boundary of the first Brillouin zone (see the inset), (c) uniform magnetic susceptibility $\chi$, and (d) thermal entanglement versus temperatures. The inset in (a) depicts the temperature $T_l$ of the peak in $c_V$ at a lower temperature versus inverse length $1/L$, with the dashed line a guide for the eye. Note that the high-temperature scale in $c_V$ stays at $T_h\sim 0.5$. The HTSE and Padé data is taken from Refs. [58, 59] and BDMC from Ref. [60].

Figure 9

(a) Finite-temperature phase diagram of TLH on a YC4 lattice, which consists of high-$T$ “gas,” intermediate-$T$ liquid, and low-$T$ solid states. (b)–(e) show the static structure factor at four representative temperature points (left to right, from high temperature to low), and (f)–(i) exhibit corresponding bond textures at finite $T$, and the values are the bond energies to whose absolute value the thickness of nearest-neighbor bonds are proportional. In (i), we also show the bond energy distribution at $T=0$ calculated by DMRG, where the agreement between XTRG results (left half) at low-$T$ and DMRG data (right half) is apparent, and one can recognize the strip solid structure quite distinctly. Dashed lines are guide for the eye, representing in (e) the line connecting two inequivalent $M$ points (labeled as $M_1$ and $M_2$) with enhanced intensity. In (f), we label the $X$ and $Y$ directions, as well as the way we wrap the lattice into a $Y$ cylinder.

Figure 10

(a) Bipartite MPO entanglement entropy $S_E$ across the center of the system of a spin-$1/2$ Heisenberg chain (length $L=100$, OBC). (b) Entanglement spectra obtained in the center of the system versus a wide range of temperatures presented as an RG flow in energy scales, where lines from the same symmetry sectors, i.e., $S=0,1,2,\dots$, are plotted in the same color. Vertical markers depict different temperature regimes (see text). The inset enlarges the region around $\beta =100$, and the labels $d^{*}$ ($d$) indicate the degeneracy in the RG fixed point ES in terms of individual multiplets (or states), respectively.

Figure 11

Bipartite MPO entanglement $S_E$ across the center of the system in an OBC Heisenberg chain of length $L=200$. At low temperatures ($\beta \ge 10$), $S_E$ diverges logarithmically. The line in the main panel for large $\beta$ (small $T$) is a fitting to $S_E=c/3\ln \beta +\text{const}$, with central charge $c\simeq 0.999$ determined by fitting the range $\beta >12$. The left inset shows the same data on a log-log plot, which emphasizes a power law of entanglement versus $\beta$ for $\beta \lesssim 1$. The purple dashed lines represent the fit in Eq. (13) with $\alpha =0.05$, while the yellow dashed lines represent ${\beta }^2$ as a guide to the eye. The right inset shows the slope of the log-log data in the left inset, i.e., the power-law exponent $\gamma$ of the entanglement versus $\beta$.

Figure 12

Bipartite MPO entanglement across the center of the system in the SLH, up to width $W=6$. The system length is fixed to $L=10$ for width $W=2$, 3, 4, 5 and to $L=12$ for $W=6$. For $W\le 4$, the entanglement entropy in the center of the system is well converged by retaining $D^{*}=500$ bond multiplets. For $W=5$ and 6, we show data with $D^{*}(D)=300$ (1000) and 550 (2000) (color matched lines) as well as extrapolated data (symbols) based on a quadratic polynomial extrapolation in $1/D^{*}\to 0$ that includes additional data between $D^{*}=300$ and $D^{*}=550$. The inset presents the entropy per leg ($S_E/W$) on a log-log scale, where the power-law increase for $\beta \ll 1$ again clearly follows the large-temperature fit in Eq. (13) with approximate exponent $\gamma \simeq 2$, using $\alpha =0.05$. The dashed line shows $S_E\sim \beta$ as a guide to the eye for reference.

Figure 13

Finite-temperature phase transition in the $XXZ$ Heisenberg model at $\mathrm{\Delta }=5$, keeping $D=600$ states in U(1) XTRG. (a) Entanglement landscape versus temperatures $T/\mathrm{\Delta }$ and bond indices on a $L=16$, $W=8$ open square lattice which exhibits pronounced peaks around the exact transition temperature $T_c=0.56$ [97] [dashed black lines in (b)–(d)]. The dash-dotted line in (b) depicts the maximal entropy point for each bond, which converges to ${T_c}^S\simeq 0.552$ in the central region of the system. (c) Block entanglements $S_E$ at bond $[(N+W)/2]$ cutting across the center of the system (where $N=LW$ and $[.]$ means the integer part) and also the specific heat $c_V$ with a comparison to QMC data. The data is shown for various widths $W=5$, 6, 7, 8 with a fixed aspect ratio, using $L=2W$. Both $S_E$ and $c_V$ curves show peaks near $T_c$ (vertical dashed line). The horizontal dashed line represents constant $\ln 2$, manifesting the global ${\mathbb{Z}}_2$ symmetry which can be spontaneously broken at low temperatures if applying a small pinning field, otherwise. The peak positions ${T_c}^S$ of the $S_E$ and $c_V$ are collected and shown versus $1/W$ in the inset. The value of ${T_c}^S$ from $S_E$ suffers smaller finite-size corrections than that obtained from the specific heat. (d) Binder ratio $U_4$ for various system sizes $L=2W$. XTRG data is plotted with open circles and looper QMC data with crosses with perfect agreement between the two data sets. The cross point determines ${T_c}^X$. By enlarging the data in the lower left inset, we obtain the crossing points ${T_c}^X(W)$ for pairs of consecutive system sizes $W$, as plotted versus $1/W$ in the upper right inset, again with good agreement between the XTRG and the QMC data. From this we obtain ${T_c}^X\simeq 0.554$ for our largest systems, which still trends towards $T_c$ (horizontal dashed line) with increasing system size and, hence, agrees with $T_c$ to within an error bar of 1%.

Figure 14

(a) The local identity and spin operator for the SU(2) spin-symmetric Heisenberg model, which transform according to $S\equiv q=0$ and $q=1$, respectively. All relevant local operators, including the identity operator, can be combined into the rank-3 tensor $X$ [(b)]. (c) The super-MPS described by the rank-3 tensor $A$ can be contracted with the operator index in $X$ to form the rank-4 local tensor $T$ in the SU(2)-invariant MPO [(d)]. (e) Typical sequence that occurs in the construction of the MPO via the super-MPS that describes a specific individual interaction term ${\widehat{S}}_i·{\widehat{S}}_j^{†}$. This demonstrates how the quantum number of the IROP $S$ simply stretches like a string (red line) along the virtual bonds in between the two sites $i$ and $j$ where the spin operators act. The values ${\alpha }_1,{\alpha }_2,\dots$ on given $A$-tensors are reserved for this very specific interaction term, where, in general, the indices ${\alpha }_i>2$ are not all the same. Note that the arrows on the red line are reversed with respect to site $j$, which thus needs to incorporate the Clebsch-Gordan coefficients $(1,1;0)$ that combine the $S=1$ multiplet of the spinor $S$ with its dual (also $S=1$) into a scalar. It is this CGC on the $A$ tensor of site $j$ within the super-MPS that takes properly care of the dagger in the scalar product ${\widehat{S}}_i·{\widehat{S}}_j^{†}$. Similarly, also the coupling strength $J$ is encoded with the $A$ tensor in the super-MPS.

Figure 15

Comparison between the first-order Trotter-Suzuki and the series-expansion initialization schemes on an $L=12$ Heisenberg chain. Two sets of data are shown from calculations starting with ${\tau }_0=0.1$ (symbols) and ${\tau }_0=0.01$ (lines). After the initialization, the MPO $e^{-{\tau }_{0}H}$ is fed into an XTRG evolution, where the number of retained bond states is set to $D^{*}=150$. The same color represents the same type of initialization; i.e., blue represents Trotter initialization, while green, red, and yellow represent initialization by SETTN with bond dimensions $D_0^{*}=2$, 3, 4 ($D_0=4$, 5, 8), respectively.

Figure 16

(a) Relative error of free energy for the Heisenberg model on a $4\times 4$ open square lattice ($D_H^{*}=6$), with various initial $\tau$ and $D_0^{*}$ values in the XTRG scheme. Here, the $\rho (\beta ={\tau }_0)$ is initialized via SETTN [cf. Eq. (7)]. (b) Similar to (a) but using the lowest-order, i.e., linear, expansion of $\rho ({\tau }_0)=\mathbb{I}-{\tau }_{0}H$ for the initialization of $\rho ({\tau }_0)$ instead of SETTN [only the green curve at ${\tau }_0={10}^{-2}$ uses SETTN and is copied from (a) for direct reference]. Here, we use $D^{*}=D_0^{*}$ for all temperatures down to $\beta ={10}^{-2}$ (black dashed line), where we switch to $D^{*}=200$. Purple dashed lines indicate linear behavior, i.e., $|\delta F/F|\propto \beta$.

Figure 17

(a) The environment tensor $C_E$ in two-site variational optimization of compressing the product of two MPOs into a single MPO, i.e., $C≔A*B$ [cf. Eq. (d2)]. Arrows on the horizontal lines demonstrate orthonormalization or, equivalently, canonicalization of the MPO. The environment $C_E$ can be evaluated by contracting the cluster which consists of two environment tensors $V_L$ and $V_R$, as well as relevant MPO tensors. (b) The environment tensor $C_E$ in two-site variational scheme of rewriting the sum of two MPOs into a single MPO, i.e., $C≔A+B$ [cf. Eq. (d4)]. An asterisk inside a box indicates a “daggered” tensor. (c) Update the left environment tensor $V_L$ with three local tensor $A_i$, $B_i$, and $C_i$ of MPOs $A$, $B$, and $C$, respectively. The computational costs of the three substeps (from left to right) are $\mathcal{O}(D_a^2{D}_b{D}_c)$, $\mathcal{O}(D_a{D}_b^2{D}_c)$, and $\mathcal{O}(D_a{D}_b{D}_c^2)$, respectively. $V_R$ can be updated in a similar way. (d) Update of the local tensor $T_i^{*}$ and $T_{i+1}^{*}$ of MPO $C$, by performing SVD on either environment tensor $C_E$ or $(C_A+C_B)$.

Figure 18

Relative errors of the free energy fully within SETTN in (a) an $L=16$ Heisenberg chain and (b) a $4\times 4$ Heisenberg model on the square lattice (OBC). The data is computed along various $\beta$ grids, including Maclaurin ($M$) expansion around $\beta =0$ and Taylor expansion in linear ($L$) and exponential ($X$) $\beta$ scales. In (a), we use $D^{*}=200$, resulting in relative CPU run times $M:L:X=2.34:5.25:1$. In (b), we use $D^{*}=200$, 400. For $D^{*}=200$, this results in relative CPU run times $M:L:X=2.1:4.95:1$ and, for $D^{*}=400$, $M:L:X=1.4:4.2:1$.

Figure 19

(a) Block entanglement landscapes versus temperature $T/\mathrm{\Delta }$ and bond index for the $XXZ$ model on a square lattice (using $\mathrm{\Delta }=5$, $L=10$, $W=5$ keeping up to $D=500$ states, i.e., similar analysis as in Fig. 13 but for a smaller system size). (b) Top view and (c) side view versus $T/\mathrm{\Delta }$ of the same data as in (a). The black dashed line represents the exact value $T_c=0.56$, and the dash-dotted line is tracking the “ridge” of entanglement landscape (i.e., maxima versus $T$ on each bond) in (a); the latter coincides with $T_c$ in the central bonds of the system. Incidentally, the envelope of entanglement curves also shows a peak around $T_c$ (with the resolution $\delta T/\mathrm{\Delta }\simeq 0.02$).