Thermal tensor renormalization group simulations of square-lattice quantum spin models

In this work, we benchmark the well-controlled and numerically accurate exponential thermal tensor renormalization group (XTRG) in the simulation of interacting spin models in two dimensions. Finite temperature introduces a finite thermal correlation length $\xi$, such that for system sizes $L\gg \xi$ finite-size calculations actually simulate the thermodynamic limit. In this paper, we focus on the square lattice Heisenberg antiferromagnet (SLH) and quantum Ising models (QIM) on open and cylindrical geometries up to width $W=10$. We explore various one-dimensional mapping paths in the matrix product operator (MPO) representation, whose performance is clearly shown to be geometry dependent. We benchmark against quantum Monte Carlo (QMC) data, yet also the series-expansion thermal tensor network results. Thermal properties including the internal energy, specific heat, and spin structure factors, etc. are computed with high precision, obtaining excellent agreement with QMC results. XTRG also allows us to reach remarkably low temperatures. For SLH, we obtain an energy per site $u_g^{*}\simeq -0.6694\left(4\right)$ and a spontaneous magnetization $m_S^{*}\simeq 0.30\left(1\right)$ already consistent with the ground-state properties, which is obtained from extrapolated low-$T$ thermal data on $W\le 8$ cylinders and $W\le 10$ open strips, respectively. We extract an exponential divergence versus $T$ of the structure factor $S\left(M\right)$, as well as the correlation length $\xi$, at the ordering wave vector $M=(\pi ,\pi )$, which represents the renormalized classical behavior and can be observed over a narrow but appreciable temperature window, by analyzing the finite-size data by XTRG simulations. For the QIM with a finite-temperature phase transition, we employ several thermal quantities, including the specific heat, Binder ratio, as well as the MPO entanglement to determine the critical temperature $T_c$.

Figure 1

Various MPO paths utilized in XTRG simulations that map the 2D lattices into quasi-1D system with long-range interactions, including (a) the snakelike, (b) zigzag, (c) diagonal, and (d) slash paths. The line width visualizes the low-temperature bond entanglement $S_E$ along the MPO obtained on the $\mathrm{OS}6\times 6$ and $\mathrm{OS}6\times 12$ lattices (at $T\simeq 0.06$), where we used a width of $w=(4S_E-11)\mathrm{pts}$, yet enforcing $w\ge 1$ for visibility.

Figure 2

Relative errors of free energy $f$ vs $1/D^{*}$, $D^{*}$ the number of multiplets kept, for the $\mathrm{YC}4\times 4$ SLH at $T\simeq 0.06$, and the ED data are taken as the exact reference. Four mapping paths are compared, including the snakelike, slash, diagonal, and zigzag. The inset shows the maximum of $S_E$ over all MPO bonds vs $T$. One can observe that $S_E$ values coincide in the essentially equivalent snakelike and zigzag paths, and are significantly smaller compared to the diagonal and slash paths.

Figure 3

The free energy $f$ of SLH on the (a) $\mathrm{OS}6\times 6$ and (b) $\mathrm{OS}6\times 12$ lattices at $T\simeq 0.06$, obtained for the four different MPO paths in Fig. 1 by retaining $D^{*}=100$ to 500 multiplets in all cases. (c)–(f) show comparisons of the free energy $f$ vs $T$ for all paths using $D^{*}=500$ for OS [(c) and (d)] and YC [(e) and (f)]. Here $f_{\min }\left(T\right)$ represents the minimal value amongst all four paths at any given $T$.

Figure 4

Internal energy $u$ and specific heat $c_V$ for [(a) and (b)] $\mathrm{OS}6\times 6$, and [(c) and (d)] $\mathrm{YC}6\times 12$ SLH systems. The insets zoom into the low-$T$ data, where the SETTN and XTRG data for various $D^{*}$ are shown to agree excellently with QMC. In the legends, $u^a$ and $u^b$ refer to the two schemes in Eqs. (7a) and (7b), respectively. The specific heat $c_V$ in the lower panels is obtained from $u$ using Eq. (8). The vertical dashed line represents the temperature scale $T_S\sim 0.6$ in SLH.

Figure 5

(a) Internal energy $u$ on the $\mathrm{OS}L\times L$ lattices up to $L=9$ calculated by XTRG keeping up to $D^{*}=1000$ multiplets. The data are obtained in three different ways (total, center, and “torus”) as described in the main text, with total shown in (a), and finite size scaling of the extrapolated data vs $1/L^2$ shown in (b). In order to reduce the finite $D^{*}$ effects, we extrapolate the internal energy $u$ to $1/D^{*}=0$, as seen in the inset of (a). In (b), we collect the low-temperature ($T\simeq 0.1$) data extrapolated in (a) $1/D^{*}\to 0$ and analyze it here vs $1/L^2\to 0$. For center, the error bars are estimated from $1/D\to 0$ extrapolations (see text), while for “torus,” we extrapolate the four largest system sizes (i.e., data in gray shaded area were excluded), to the thermodynamic limit, via a second-order polynomial fitting vs $1/L^2$. The horizontal dashed line represents the ground-state energy $u_g\simeq -0.6694$ from QMC [52]. For comparison, (c)–(h) analyses the internal energy $u$ of SLH on $\mathrm{YC}W\times L$ cylinders of widths $W=6,8$ and lengths $L=8,10,12$. Exemplary extrapolations of $u_{\mathrm{center}}$ vs $1/D^{*}\to 0$ are shown in (c) and (f) for $T\simeq 0.11$ and 0.06, respectively, where $u_{\mathrm{center}}$ is evaluated via a weighted average around the center as illustrated in the inset of (f) (see main text for more details). The results at $1/D^{*}\to 0$ are collected vs $1/WL$ in (d), (e), (g), and (h) (green stars). There they are also compared to similarly extrapolated data for $u_{\mathrm{tot}}$ (black squares), as well as to $u_{\mathrm{subtr}}$ (blue horizontal line) obtained by subtracting the length $L=8$ from the $L^{\prime }=12$ cylinder. With $u_{\mathrm{tot}}$ also extrapolated to $1/WL\to 0$, we find good agreement across our data towards the thermodynamic limit.

Figure 6

(a) SLH specific heat $c_V$ on the $\mathrm{OS}L\times L$ lattices on a log-log scale to emphasize the algebraic behavior at high and low temperature, obtaining, $c_V\sim T^{-2}$ and $\sim T^2$, respectively (see dashed lines as guide to the eye). The inset zooms into the peak at intermediate temperatures around $T_S\sim 0.6$. (b)–(e) show the static spin structure factors $S\left(q\right)$ on $\mathrm{OS}9\times 9$ at $T\simeq 0.20$, 0.55, 2.21, and 10.51, respectively. The grid lines demarcate the Brillouin zone, where the white dots in (e) indicate specific high-symmetry points therein.

Figure 7

The static spin structure factors $S\left(q\right)$ on the $\mathrm{OS}9\times 9$ lattice at (a) $T\simeq 12.5$ and (c) $\simeq 0.2$, calculated by XTRG. Fittings using the antiferromagnetic Ising (AFI) model and independent dimer approximation (IDA) (see main text) are presented at the right half of each panel [(b) and (d)], which enjoy excellent agreement with XTRG data on the left half of each panel [(a) and (c), respectively].

Figure 8

The finite-size analysis of spontaneous magnetization $m\left(L\right)\equiv \sqrt{S\left(M\right)/L^2}$ of the SLH on the $\mathrm{OS}L\times L$ lattices up to $L=10$. In the inset, we collect $m^2\left(L\right)$ data at $T=0.1$ where convergence vs $T$ is reached, and extrapolate it to the thermodynamic limit $1/L\to 0$ using a parabolic fit. The data in the gray shaded area were excluded from this fit. From this, we estimate the value for the thermodynamic limit $m_S^{*}\simeq 0.30\left(1\right)$, in good agreement with the literature, $m_S\simeq 0.3070\left(2\right)$ [56].

Figure 9

(a) Static structure factor $S\left(M\right)$ vs $T$ (same data as in Fig. 8). Its derivative in (b) has a maximum which defines a specific temperature $T_f$ at which finite size effects become significant. The inset in (a) then analyzes $S_f/T_f^2$ vs $1/T_f$. The inset in (b) shows that $T_f$ extrapolates to 0 in the thermodynamic limit, via a cubic fit as shown based on the data with $L\ge 4$. (c) Correlation length $\xi$ vs $T$. The results for $S\left(M\right)$ and $\xi$ evaluated at $T_f$, are marked by black circles in (a) and (c), respectively. The data are collected and analyzed in the respective insets vs $1/T_f$ which, overall, shows good agreement with RC predictions (dashed line). For comparison, the insets in (a) and (c) also plot the data for $S\left(M\right)$ and $\xi$ vs $1/T$ for the largest system size.

Figure 10

MPO entanglement entropy $S_E$ on (a) cylinders $\mathrm{YC}4\times 8$, $\mathrm{YC}6\times 12$, and $\mathrm{YC}8\times 12$, and (b) $\mathrm{OS}L\times L$ with $L=4$ to 8. The tilted dashed lines in both (a) and (b) represent $S_E\propto aln\beta +\mathrm{const}.$, with slope $a\simeq 0.4$ (see dashed guides on top of the curves). The vertical dotted line represents the temperature scale $T_S\sim 0.6$ [e.g., see Fig. 6].

Figure 11

Specific heat of QIM for fixed $h_x=\frac{2}{3}{h}_c$ on $\mathrm{YC}W\times L$ of width up to $W=8$ and length $L=2W$ ($W=7$ curve not shown in the main panel for better readability). The XTRG results retaining up to $D=240$ states coincide with the QMC reference data. In the inset, we collect the peak position $T_c^{*}$ of $c_V$ curves calculated by QMC and XTRG, which also coincide, and extrapolate towards the exact critical temperature in the thermodynamic limit $x\equiv 1/W^2\to 0$ by a second order polynomial fit, having $T_c^{*}\left(x\right)\simeq -14.5x^2+2.8x+T_c^{*}\left(0\right)$, with an extrapolated value of $T_c^{*}\left(0\right)\simeq 0.4184$ (a polynomial fit in $1/W$ leads to a similar value).

Figure 12

Binder ratio curves of QIM for the (a) $h_x=\frac{2}{3}{h}_c$ and (b) $h_x\simeq 0.82h_c$ for $\mathrm{YC}W\times L$ systems with a fixed aspect ratio $L/W=2$. The data in (a) is compared to the QMC data on the same lattice, and in (b) the $D=200$ and 300 Binder ratio curves lie on top of each other, showing the convergence of results vs bond dimension $D$ in both cases. The left inset in (a) and the upper inset in (b) zoom in the region near the crossing points $T_c^{*}$, and the rest insets show subsequent second-order polynomial extrapolations of $T_c^{*}$ to $1/W^2\to 0$. In (a), the $W=7$ data, not shown in the main panel for better readability, are included for the right inset, and a second-order polynomial fitting as shown yields an extrapolated $T_c^{*}\left(0\right)\simeq 0.4212$, which is in excellent agreement with $T_c\simeq 0.4239$ in the thermodynamic limit [25]. In (b), we obtain an extrapolated $T_c^{*}\left(0\right)\simeq 0.3216$, again in good agreement with the QMC [70] and iPEPS results [28, 29].

Figure 13

(a) Entanglement landscape of QIM thermal states on the $\mathrm{OS}10\times 10$ lattice, for the same $h_x=\frac{2}{3}{h}_c$, vs bond indices and temperature. (b) Top view, showing that the peak temperatures at central bonds (away from boundaries) locate right at the critical temperature $T_c$ (horizontal dashed line). (c) $S_E$ vs $T$ on various $\mathrm{OS}L\times L$ lattices, cut at a central bond [indicated by the black arrow in (b)] with maximal $S_E$. The blue vertical dashed lines in both the main panel and right inset indicate $T_c\simeq 0.4239$ in the thermodynamic limit [25]. The left inset shows the peak temperature $T_c^S$ vs inverse system size $1/L^2$, which approach the true critical temperature $T_c$. $D=240$ bond states are kept in the calculations. The right inset shows the entanglement difference vs $T$ between consecutive system sizes i.e., $\mathrm{OS}L\times L$ and $\mathrm{OS}{L}^{\prime }\times L^{\prime }$ with $L^{\prime }=L+2$ bearing in mind even-odd effects. The gray vertical dashed line corresponds to the first-order extrapolated value at $1/L^2\to 0$ in the left inset.

Figure 14

(a) Comparison of XTRG and Trotter linear evolution (with swap gates) schemes in the benchmark calculations of an SLH model on $\mathrm{OS}4\times 4$. (b) Computational runtime of the Trotter approach relative to XTRG, which scales roughly as $1/D$ (with the dashed line a guide to the eye). (c) The truncation error analysis in the linear Trotter calculations: the swap gates contribute significantly larger (over two orders of magnitude) truncation errors than those of imaginary-time evolution gates.

Figure 15

Truncation error (solid lines; left axis set) and entanglement entropy $S_E$ (dashed lines; right axis set) at a low temperature ($T\simeq 0.06$). The data are obtained by retaining $D^{*}=500$ multiplets. We compared the four different mapping paths in Fig. 1, namely, the snakelike, zigzag, diagonal, and slash, on (a) $\mathrm{OS}6\times 6,$ (b) $\mathrm{OS}6\times 12,$ (c) $\mathrm{YC}6\times 6$, and (d) $\mathrm{YC}6\times 12$. The entanglement entropy is precisely also the data that was used for the width of the lines in Fig. 1 in a graphically more organized way.

Figure 16

The energy per site $u$ vs truncation error $\delta \rho$ on (a) YC6 and (b) YC8 geometries at $T\simeq 0.06$ [same as in Fig. 5 of the main text]. A clear linear $u$ vs $\delta \rho$ relation is observed and employed to perform the extrapolation. The thus extrapolated $u$ values are in a very good consistency to $1/D^{*}\to 0$ analysis (asterisks), with relative error $\sim 0.1\%$. Insets in (a) and (b) show $\delta \rho$ vs $1/D^{*}$ on a log-log scale, showing polynomial scaling.