Detecting a $Z_2$ topologically ordered phase from unbiased infinite projected entangled-pair state simulations

We present an approach to identify topological order based on unbiased infinite projected entangled-pair states simulations, i.e., where we do not impose a virtual symmetry on the tensors during the optimization of the tensor network ansatz. As an example we consider the ground state of the toric code model in a magnetic field exhibiting $Z_2$ topological order. The optimization is done by an efficient energy minimization approach based on a summation of tensor environments to compute the gradient. We show that the optimized tensors, when brought into the right gauge, are approximately $Z_2$ symmetric, and they can be fully symmetrized a posteriori to generate a stable topologically ordered state, yielding the correct topological entanglement entropy and modular S and U matrices. To compute the latter we develop a variant of the corner-transfer matrix method, which is computationally more efficient than previous approaches based on the tensor renormalization group.

Figure 1

The square lattice for the toric code Hamiltonian considered in this paper. The black dots on the edges represent spin-$\frac{1}{2}$ particles. The spins involved in the action of a single $A_v$ or $B_p$ operation is shown in blue and red, respectively.

Figure 2

The ground state $|{\mathrm{\Psi }}_o\rangle$ cannot be turned into $|{\mathrm{\Psi }}_e\rangle$ under the action of the Hamiltonian, resulting in two distinct degenerate ground states. The two ground states can be distinguished by the nonlocal operator $W^{\left(z\right)}\left(C\right)={\prod }_{i\in C}{\sigma }_i^z$ along the line $C$ which measures the parity of the number of loops cutting the line $C$.

Figure 3

(a) Mapping the two-dimensional lattice to an iPEPS representation with a checkerboard pattern of $A$ (orange) and $B$ (blue) tensors. These tensors are related by a ${90}^{\circ }$ rotation. (b), (c) Tensors used for the construction of the exact ground state of the toric code model.

Figure 4

Tensor network diagrams describing the (symmetric) CTM algorithm used in this paper. (a) Infinite tensor network consisting of the checkerboard pattern of tensors $a$ and $b$, constructed by combining the iPEPS tensors $A$ with $A^{†}$ and $B$ with $B^{†}$, respectively, as shown in (c). The infinite tensor network surrounding a bulk tensor $b$ is effectively represented by an environment consisting of four corner tensors $C_{\times }$/$C_{\diamond }$ and four edge tensors $T$ as depicted in (b). The thin and thick lines have a dimension $D^2$ and $\chi$, respectively. (d) The initial boundary tensors are constructed from the iPEPS bulk tensors by projecting the legs in both the bra and ket layer onto the first element, depicted by the black dots. (e) Definition of the two different corner tensors and orientation of the $T$ tensor. (f), (g) Growth and renormalization step in the CTM algorithm (see text).

Figure 5

The tensor network diagrams to evaluate the toric code Hamiltonian using the CTM environment tensors. The plaquette contribution is labeled as $\widehat{H_1}$ and the vertex contribution $\widehat{H_2}$.

Figure 6

Tensor network diagrams to represent the gradient of the energy with respect to tensor $A^{†}$ (only the contributions from ${\widehat{H}}_1$ are shown). (a) The derivative of the bulk and boundary tensors with respect to $\partial A^{†}$ are represented by the green tensors, e.g., $a_{\otimes }={\partial }_{A^{†}}a$. The $\otimes$ symbols in the tensor shapes mark the different locations in the lattice at which the derivative is taken. In (b) the total derivative of $\langle {\widehat{H}}_1\rangle$ is represented in terms of a summation over the different tensors introduced in (a). The first term can be directly evaluated and its contribution can be added to the total gradient $\mathcal{G}$. The other contributions are rewritten in terms of the environments ${\mathrm{\Gamma }}_{C_{\diamond }}^{\left(N\right)},{\mathrm{\Gamma }}_T^{\left(N\right)}$ and $C_{\diamond ,\otimes }^{\left(N\right)},T_{\otimes }^{\left(N\right)}$ which are evaluated recursively as shown in Fig. 7.

Figure 7

Recursive procedure to sum up all contributions to the gradient (only the contributions from ${\widehat{H}}_1$ are shown). (a) The corner gradient tensor $C_{\diamond ,\otimes }^{\left(i\right)}$ is expanded in its four contributions. The first term can be directly evaluated and its contribution added to the total gradient $\mathcal{G}$. The second term, which is taken with a prefactor 2 due to the presence of two $T$ tensors which are equivalent by symmetry, becomes a contribution ${\mathrm{\Gamma }}_T^{{(i-1)}^{\prime }}$ to the edge environment for the next iteration, and the last term generates the corner environment tensors for the next iteration, ${\mathrm{\Gamma }}_{C_{\diamond }}^{(i-1)}$. Similarly, in (b) the edge tensor ${\mathrm{\Gamma }}_T^{\left(i\right)}$ is expanded into two terms, where the first term can be directly evaluated and added to the total gradient $\mathcal{G}$, and the second term yields a contribution ${\mathrm{\Gamma }}_T^{{(i-1)}^{\prime \prime }}$ to the edge environment for the next iteration. (Note that in the second term the orientation of the boundary legs flips.) After each iteration the different $(i-1)\mathrm{th}$ edge environment contributions obtained in the corner and the edge update step (from both ${\widehat{H}}_1$ and ${\widehat{H}}_2$) are added to obtain the edge environment ${\mathrm{\Gamma }}_T^{(i-1)}$ for the next iteration.

Figure 8

(a) A tensor with a virtual symmetry, i.e., which is invariant under the simultaneous action of $U_g$ on the virtual level. This naturally leads to the “pulling-through” condition shown in (b).

Figure 9

The PEPS tensor is brought into an approximate $Z_2$ symmetric form by performing a gauge change using a unitary matrix $q$ which minimizes ${\delta }_{Z_2}$.

Figure 10

Diagrams for the computation of the second Rényi entropy between two halves of an infinite cylinder of width $L$. (a) Representation of ${\sigma }_L$ and ${\sigma }_R$ of the left/right half of the infinite cylinder in terms of the edge tensors $T$ obtained from the CTM method. (b, c) Tensor networks representing the traces in Eq. (11). Instead of contracting a large network we diagonalize the double row of $T$ tensors shown in (d) to obtain the eigenvalue matrices $\mathrm{\Lambda }$ and $\mathcal{N}$, so that the traces in (b) can be computed as $\mathrm{Tr}\left[{\mathrm{\Lambda }}^{L/2}\right]$ and $\mathrm{Tr}\left[{\mathcal{N}}^{L/2}\right]$, respectively.

Figure 11

Variant of the CTM algorithm where the parity on the bra and ket level on each boundary can be controlled. The red lines are of dimension 2, with one state in the even (e) and odd (o) sector, respectively, and the green triangles denote a projection onto this space. A black dot corresponds to a projection onto the even-parity sector. (a) Initialization of the boundary tensors, where the boundaries of the corner tensor are projected onto the even-even sector, and the red parity legs of the $T$ tensor are kept open. (b) Coarse-graining step, in which the corner tensors are kept in the even-even boundary sector, and the red parity legs of the $T$ remain open. (c) Contraction of the network on the left yields the double-layer tensor $\mathbb{T}$ representing an infinite plane (or infinite torus when connecting the legs in a periodic way), where the parity in each layer on each boundary can be controlled via the red legs.

Figure 12

(a) Second Rényi entropy of the toric code model as a function of width $L$ of an infinite cylinder, obtained from the optimized tensors ($O$), combined with a gauge transformation to bring the tensors into an approximate $Z_2$ symmetric form ($O+GT$), and symmetrized to recover the $Z_2$ virtual symmetry ($O+GT+S$). The TEE $\gamma$ is obtained from the intersection of a linear fit (dashed lines) with the $y$ axis. At short distances the correct value of TEE is obtained after performing the gauge transformation even without symmetrizing the tensors. (b) The TEE $\gamma$ obtained from linear fits to the data between $L$ and $L+100$ on large cylinders. At large distances, only the symmetrized tensors yield the correct value for $\gamma$.

Figure 13

Error in the modular matrices with respect to the exact results, $\mathrm{\Delta }S=\left|\right|S-S_{\text{exact}}\left|\right|$ and $\mathrm{\Delta }U=\left|\right|U-U_{\text{exact}}\left|\right|$, as a function of CTM iteration $N$, corresponding to a linear system size $L=2N+2$, here for $D=2$.

Figure 14

(a) The TEE $\gamma$ and the dominant eigenvalues of the transfer matrix $\mathcal{N}$ across the phase transition as a function of the magnetic field $h_z$, obtained for $D=2$ ($\chi =30$). The dashed line shows the Monte Carlo result for the location of the phase transition [56]. (b) Correlation length $\xi$ as a function of $h_z$ across the phase transition for different values of the boundary dimension $\chi$, exhibiting a peak at the $D=2$ critical point. (c), (d) Same as in (a) and (b) for bond dimension $D=3$ [$\chi =60$ in (a)].