Variational approach to projected entangled pair states at finite temperature

The projected entangled pair state (PEPS) ansatz can represent a thermal state in a strongly correlated system. We introduce a variational algorithm to optimize this tensor network whose essential ingredient is an auxiliary tree tensor network (TTN). Since the full tensor environment is taken into account, with increasing bond dimension the PEPS-TTN ansatz provides the exact Gibbs state. Our presentation opens with a 1D version for a matrix product state (MPS-TTN) and then generalizes to PEPS-TTN in 2D. Benchmark results in the quantum Ising model are presented.

Figure 1

(a) The Trotter tensor $T_0$ in 1D. This spin operator with (red) spin indices depends on (black) bond indices connecting it with similar operators at the nearest-neighbor sites. (b) Evolution operator $U\left(\beta \right)$ is a product of $N$ time steps represented by horizontal rows of tensors. Each connecting line is a contracted index. (c) At every site the column of elementary $T_0$'s is compressed by an isometry $W$ into one tensor $T_n$ with a bond dimension $D<2^N$. (d) Evolution operator $U\left(\beta \right)$ is approximated by a matrix product operator of the compressed tensors $T_n$.

Figure 2

The huge compression in Fig. 1 is split into $n$ steps. (a) As a first step $k$ elementary Trotter tensors $T_0$ are compressed by isometry $W_1$ into tensor $T_1$ with a bond dimension $D\le 2^k$. (b) The first step is followed by iterative compressions $T_{m-1}\to T_m$ preserving the bond dimension $D$. (c) The final $T_n$ is the same as if the huge isometry $W$ were a tree tensor network (TTN) with $n$ layers of isometries: $W_1,\cdots ,W_n$. Here we show an example with $k=3$ and $n=3$.

Figure 3

(a) The partition function $Z=\langle \psi \left(\beta \right)\left|\psi \right(\beta )\rangle$ is a product of transfer matrices $t$. (b) The same as in (a) but with one of the tensors $T_n$ removed. This diagram is the tensor environment $E\left(n\right)$ for $T_n$. (c) A lower-level environment $E(m-1)$ for $T_{m-1}$ is obtained from $E\left(m\right)$. (d) Environment $E_W\left(m\right)$ for isometry $W_m$ is obtained from $E\left(m\right)$. (e) Environment $E_W\left(1\right)$ is obtained from $E\left(1\right)$. (f) Isometric environment is subject to singular value decomposition, $E_W\left(m\right)=U\lambda V^{†}$, with $D$ nonzero singular values $\lambda$, a $D^2\times D$ isometry $U$, and a $D\times D$ unitary $V$. The isometry is updated as $W_m=U{V}^{†}$ to maximize the figure of merit $Z=\mathrm{Tr}{E}_W\left(m\right)W_m^{†}$.

Figure 4

(a) Energy per site of the quantum Ising chain at the critical field $h=1$ in function of $\beta$ for different bond dimensions $D$. For larger $D$ the results remain accurate down to lower temperatures that are closer to the quantum critical point. (b) The correlation length in the tail of the ferromagnetic correlator (13). In both (a) and (b), the number of isometry layers $n=11$, the bond dimension $D=2^k=2,\cdots ,64$, and the elementary time step in the Suzuki-Trotter decomposition $d\beta =\beta /\left(2^{n-1}k\right)$.

Figure 5

(a) The elementary Trotter tensor $T_0$ on a square lattice. (b) Compression of two tensors $T_{m-1}$ into one $T_m$ with four isometries $W_m$. (c) Tensors $T_n$ are contracted into a PEPS network for the evolution $U\left(\beta \right)$. (d) A contraction of two tensors $T_n$ makes a transfer tensor $t$. (e) Contraction of transfer tensors is the partition function $Z$. (f) The partition function with one of the transfer tensors removed is a tensor environment $E_t$ for the removed tensor $t$. We enclose the ends of the free bonds by a transparent reddish sphere. Each free bond has dimension $D^2$. (g) When the environment/sphere $E_t$ is filled and contracted with one Trotter tensor $T_n$ we obtain tensor environment $E\left(n\right)$ for $T_n$. This environment is the partition function with one $T_n$ removed. (h) Environment $E\left(m\right)$ filled and contracted with two $T_{m-1}$'s and three $W_m$'s makes environment $E_W\left(m\right)$ for isometry $W_m$. (i) $E(m-1)$ is obtained from $E\left(m\right)$.

Figure 6

(a) Spontaneous magnetization $\langle Z\rangle$ in the ferromagnetic phase at the transverse field $h=\frac{2}{3}{h}_0$. Different colors correspond to bond dimensions $D=2,\cdots ,8$. Plots for $D=6,7,8$ collapse demonstrating convergence for $D\ge 6$. (b) The ferromagnetic correlator $C_R$ in Eq. (15) for $D=6$ at $h=\frac{2}{3}{h}_0$ and different $\beta$ in the ferromagnetic phase. Its correlation length diverges closer to criticality requiring a diverging $M$. For $M=68$ the longest length achieved is $\xi =290$ lattice sites.

Figure 7

(a) Spontaneous magnetization $\langle Z\rangle$ near the critical point. With increasing bond dimension $D$ the plots become “critical” enough to make a fit $\langle Z\rangle \propto {\left(\beta -{\beta }_c\right)}^{{\beta }^{\prime }}$ with a critical-point ${\beta }_c=0.5898$ and a critical exponent ${\beta }^{\prime }=0.198$. (b) The ferromagnetic correlation function (15) for $D=6$ at ${\beta }_c$. In the intermediate range $1<R<30$ the correlator is fitted well by the power law $C_R\propto R^{-\eta }$ with the critical exponent $\eta =0.265$.

Figure 8

Comparison between the present variational method in (a) and the direct imaginary time evolution in (b). The TTN ansatz is exponentially more compact than MPS since the number of different isometries is only logarithmic in the number of time steps. On top of this, in the present variational method the isometries in TTN are optimized by repeated up- and down-sweeps to provide only the most accurate state at the final $\beta$, while in the time evolution all intermediate states between 0 and the final $\beta$ need to be accurate compromising the accuracy of the final one.

Figure 9

Spontaneous magnetization $\langle Z\rangle$ in function of $\beta$ obtained with different algorithms. The direct imaginary time evolution produces a plot with a discontinuity near the second-order transition due to insufficient $M$. Adding a tiny symmetry breaking bias $\delta ={10}^{-6}$ smooths the transition making the necessary $M$ finite. Both evolution curves come from Ref. [15]. They are compared with the PEPS-TTN curve from Fig. 7. Here all the results for $D=6$.

Figure 10

On the left, a planar version of the network in Fig. 5 representing the partition function. Here each (black) bond represents two bond indices in Fig. 5. Its dimension is $D^2$. This infinite contraction cannot be done exactly, hence it is approximated by the finite network on the right. The finite corner matrices $C$ and top tensors $T$ effectively represent corresponding infinite sectors of the network on the left separated by the dashed blue lines. Their (red) environmental bonds have dimension $M$. The environmental tensors $C$ and $T$ should be such that to the transfer tensor $t$ in the center its effective environment on the right appears the same as its exact environment on the left as much as possible. They are obtained by iterating until convergence of the corner matrix renormalization in Fig. 11.

Figure 11

The corner $C$ and top $T$ are obtained by repeating until convergence of a renormalization procedure. The procedure has four steps. In the first step, the tensors $C,T,t$ are contracted to an enlarged corner $C^{\prime \prime }$. In the second step, the symmetric $M{D}^2\times M{D}^2$ matrix $C^{\prime \prime }$ is diagonalized and its $M$ eigenvectors with the largest eigenvalues define an isometry $Z$. The diagonalization that scales like $M^3{D}^6$ is the leading cost of this variant of corner matrix renormalization. In the third step, $Z$ is used to renormalize/truncate the indices of $C^{\prime \prime }$ back to the original dimension $M$ giving a new (diagonal) corner $C^{\prime }$. In the fourth step, the same $Z$ renormalizes the contraction of $T$ with $t$ to a new $T^{\prime }$. The four-step procedure is repeated until convergence of the $M$ leading eigenvalues of $C^{\prime \prime }$.

Figure 12

(a) A finite tensor network approximating the infinite environment $E_t$ in Fig. 5. Notice that the corner matrix $C$ can be singular-value-decomposed, $C=U\lambda V^{†}$, and then absorbed into the neighboring top tensors $T$. With $X=\sqrt{\lambda }{U}^{†}TU\sqrt{\lambda }$ and $Y=\sqrt{\lambda }{V}^{†}TV\sqrt{\lambda }$ (here the matrix products are understood in the environmental bond indices) we obtain an equivalent network in panel (b). This network is a matrix product state (MPS) of four tensors $X,Y,X,Y$ with a bond dimension $M$, hence in principle a finite $M$ could suffice to represent the exact $E_t$ even at the critical point.